Bayesian Wavelet-packet Historical Functional Linear Models

TL;DR

This paper introduces a novel Bayesian approach using wavelet-packet transformations for historical functional linear models, enabling improved inference on temporal exposure-outcome relationships.

Contribution

It proposes the first use of wavelet-packet transformations in Bayesian HFLMs and evaluates two Bayesian inference methods in this context.

Findings

Wavelet-packet HFLM performs well in simulations.

Bayesian inference procedures provide reliable uncertainty quantification.

Application reveals effects of particulate matter exposure on heart rate variability.

Abstract

Historical Functional Linear Models (HFLM) quantify associations between a functional predictor and functional outcome where the predictor is an exposure variable that occurs before, or at least concurrently with, the outcome. Current work on the HFLM is largely limited to frequentist estimation techniques that employ spline-based basis representations. In this work, we propose a novel use of the discrete wavelet-packet transformation, which has not previously been used in functional models, to estimate historical relationships in a fully Bayesian model. Since inference has not been an emphasis of the existing work on HFLMs, we also employ two established Bayesian inference procedures in this historical functional setting. We investigate the operating characteristics of our wavelet-packet HFLM, as well as the two inference procedures, in simulation and use the model to analyze data on the impact of lagged exposure to particulate matter finer than 2.5$\mu$g on heart rate variability in a cohort of journeyman boilermakers over the course of a day's shift.

Tables

Table 1. Table 1: RMISEs, average percent of energy preserved in X W P superscript X subscript 𝑊 𝑃 \textbf{X}^{W_{P}} (PE), and computation time for the simulation settings lagged ( L 𝐿 L ), cumulative ( C 𝐶 C ), time-specific ( T 𝑇 T ), and delayed time-specific ( D 𝐷 D ). Table values represent averages taken over 200 simulated datasets. RC denotes retained coefficients. Time is averaged over all simulated datasets, over all settings.

N	T	RC		RMISE					PE		Time (s)
				${𝜷}_L$	${𝜷}_C$	${𝜷}_T$	${𝜷}_D$
50	64	25%		0.015	0.012	0.095	0.045		73.67%		516.2
		50%		0.018	0.017	0.056	0.017		88.42%		609.6
	128	25%		0.010	0.008	0.040	0.008		74.91%		1035.7
200	64	25%		0.009	0.007	0.088	0.009		73.55%		515.6
		50%		0.010	0.008	0.045	0.008		88.34%		615.5
	128	25%		0.004	0.004	0.036	0.003		75.63%		1032.1

Table 2. Table 2: Point-wise and joint credible interval coverage probabilities for the simulation settings lagged ( L 𝐿 L ), cumulative ( C 𝐶 C ), time-specific ( T 𝑇 T ), and delayed time-specific ( D 𝐷 D ). All intervals are at the 95% level. Table values represent averages taken over 200 simulated datasets. RC denotes retained coefficients.

N	T	RC		Point-wise Interval					Joint Interval
				${𝜷}_L$	${𝜷}_C$	${𝜷}_T$	${𝜷}_D$		${𝜷}_L$	${𝜷}_C$	${𝜷}_T$	${𝜷}_D$
50	64	25%		0.811	0.870	0.524	0.565		0.993	0.997	0.770	0.778
		50%		0.923	0.942	0.826	0.934		0.999	0.999	0.984	0.999
	128	25%		0.855	0.885	0.767	0.940		1.000	1.000	0.978	1.000
200	64	25%		0.724	0.784	0.452	0.735		0.957	0.986	0.673	0.957
		50%		0.865	0.890	0.718	0.930		0.999	0.999	0.933	1.000
	128	25%		0.884	0.881	0.590	0.897		0.999	0.999	0.906	1.000

Equations

$$y_i(t)=\alpha (t)+{\int }_{v\in \mathcal{V}}{x}_i(v)\beta (v,t)dv+E_i(t),$$

$$y_i(t)=\alpha (t)+{\int }_{\{v\le t\}}{x}_i(v)\beta (v,t)dv+E_i(t).$$

$$\text{Y}=\text{X}\mathbf{\beta }+\text{E},$$

$$\boldsymbol{\beta}=\left(\begin{array}[]{cccc}\beta(v_{1},t_{1})&\beta(v_{1},t_{2})&\cdots&\beta(v_{1},t_{T})\\ 0&\beta(v_{2},t_{2})&\cdots&\beta(v_{2},t_{T})\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\beta(v_{V},t_{T})\end{array}\right)$$

$${\beta }_{s\ell ,jk}^{W_P}|\text{rest}\sim 1(\ell \le k){\gamma }_{s\ell ,jk}N({\mu }_{s\ell ,jk},{ϵ}_{s\ell ,jk})+(1-{\gamma }_{s\ell ,jk})d_0,$$

$${\gamma }_{s\ell ,jk}|\text{rest}\sim Bern({\alpha }_{s\ell ,jk}),$$

$$P\left({\sigma }_{jk}^2|\text{rest}\right)$$

$${\tau }_{s\ell ,j\cdot }|\text{rest}$$

$${\pi }_{s\ell ,j\cdot }|\text{rest}$$

$$q^{(m)}=\underset{(v,t)}{\max }\left|\frac{{\beta }^{(m)}(v,t)-\hat{\beta }(v,t)}{\widehat{\text{St.Dev}}\left\{\hat{\beta }(v,t)\right\}}\right|,\text{ s.t. }v\le t.$$

Full text

Keywords: Function-on-function regression; Functional data analysis; Wavelet regression; Bayesian Statistics; Functional Inference.

1 Introduction

Historical Functional Linear Models (HFLMs) are used to analyze the relationship between a functional “exposure” and a functional “outcome” where only exposures occurring in time before or concurrently with the outcome can affect the outcome. HFLMs are a special case of the Function-on-Function Regression (FFR) model which fits an unconstrained surface and is therefore inappropriate for modeling functional predictors that are lagged exposures. For example, suppose that for subject $i$, $x_{i}(v)$ represents levels of a pollutant sampled on a grid $v\in\mathcal{V}$ and $y_{i}(t)$ represents measurements of heart rate variability (HRV) sampled on a grid $t\in\mathcal{T}$. A general FFR model with no constraints takes the form

[TABLE]

where the surface $\beta(v,t)$ is the primary quantity of interest for estimation and $E_{i}(t)$ is typically assumed to be distributed as a Gaussian Process. For example, see Ivanescu et al. (2015), Meyer et al. (2015), Morris (2015), Scheipl, Staicu, and Greven (2015), Scheipl and Greven (2016), Kim et al. (2018), and references therein.

However, when applied to the HRV and pollutant example, a model like (1) would allow for HRV measurements at time $t$ to be associated with pollutants occurring both prior to time $t$, as well as after time $t$ despite the implausibility of such a relationship. The HFLM addresses this issue, constraining $\beta(v,t)$ to prevent such spurious associations by limiting the integration in (1) to the set of coefficients such that $\{v\in\mathcal{V},t\in\mathcal{T}:v\leq t\}$. The basic HFLM takes the form

[TABLE]

The observed data is usually discrete, so that $\beta(v,t)$ can be expressed as a matrix of coefficients. Thus the problem reduces to constraining the estimate of $\beta(v,t)$ to be non-zero for the upper triangle of the matrix.

Several authors explore ways of implementing the constraint in (2). Malfait and Ramsay (2003) propose the use tent-shaped basis functions with support over a two-dimensional region. They estimate the surface using a multivariate linear model approximation to a finite dimensional model. Thus after dimension reduction via the basis-space expansion, they use least squares to estimate $\beta(v,t)$. Harezlak et al. (2007) also use basis functions defined over a two-dimensional region but specify a large number of basis functions and penalize the fit. The authors consider both the LASSO and $L_{2}$-norm penalty on triangular basis functions, using restricted maximum likelihood (REML) for the latter. Both Malfait and Ramsay (2003) and Harezlak et al. (2007) allow for a pre-defined lag beyond which the effect of exposure is zero, further constraining the surface to a trapezoidal region defined by $\{v,\in\mathcal{V},t\in\mathcal{T}:t-\Delta\leq v\leq t\}$ for some pre-defined lag $\Delta$. Kim, Şentürk, and Li (2011) take the constraint further by proposing a Recent History Functional Linear Model where the surface is constrained to a trapezoidal region defined by $\{v,\in\mathcal{V},t\in\mathcal{T}:t-\Delta_{1}\leq v\leq t-\Delta_{2}\}$ for $0<\Delta_{1}<\Delta_{2}<T$. The authors estimate the constrained surface with a varying coefficient model representation using B-spline basis functions, although they suggest Fourier, truncated power, and Eigen basis functions can also be used.

More recently, Pomann et al. (2016) examine two HFLMs that allow for multiple functional predictors with estimation constrained to a fixed window similar to that used in Harezlak et al. (2007). The authors implement two approaches to estimation using semi-local smoothing, which performs point-wise estimation, and global smoothing, which smooths over $\mathcal{T}$ globally. The methods select smoothing parameters via cross-validation and REML, respectively, and use B-spline basis expansions to model the functional form. In each of these methods, the covariance of the error term, $E_{i}(t)$, is assumed independent or to have “working” independence.

To the best of our knowledge, the existing body of work on HFLMs is largely limited to spline-based methods which can over smooth signals and peaks in spiky and irregular data. Further, the current literature does not consider inferential procedures, focusing on estimation and model fit criterion instead. As such, the performance of the proposed methods with respect to uncertainty quantification is not clear. Wavelet-based functional regression models, such as the work of Morris and Carroll (2006) and Malloy et al. (2010), consider the function-on-scalar and scalar-on-function regression settings, respectively, in the Bayesian context. Meyer et al. (2015) extend Morris and Carroll (2006) to the FFR case using wavelets for the basis function of the outcome and wavelet-principal components for the expansion of the predictor. One advantage of the wavelet-based framework is that it does not require the assumption of independence in the data-space. However, their approach estimates an unconstrained surface and therefore is inappropriate for modeling lagged exposures. The wavelet domain Meyer et al. (2015) use, which results from a discrete wavelet transformation, lacks a convenient relationship between the wavelet coefficients and the time domain, thus making it inefficient for use in an HFLM. Furthermore, the wavelet-principal components basis function does not preserve the temporal relationship between exposure and outcome in the wavelet domain. Thus to implement an HFLM using wavelets, a different basis function is warranted. Wavelet-packets, which result from a discrete wavelet-packet transformation, are a variant of the wavelet basis function that have a convenient relationship to the time domain that we can exploit to implement a wavelet-packet based HFLM.

In this work, we propose a novel use of wavelet-packets that allows us to build a Bayesian wavelet-space HFLM. We formulate our model within the framework of Morris and Carroll (2006) and Meyer et al. (2015), thus our method does not require the assumption of independence in the data space. A benefit of the Bayesian context is that we can implement several Bayesian techniques for conducting multiplicity adjusted inference including joint credible intervals to quantify uncertainty and the Bayesian false discovery rate to identify exposure lags and exposure times that are associated with the outcome. We assess the operating characteristics of the methodology and inference procedures in simulation and present an application to data from a study of journeyman boilermakers exposed to varying levels of particulate matter during the course of the day. The data consists of five-minute assessments of HRV, as defined by standard deviation of the normal-to-normal intervals (SDNN) at each time-point, and particulate matter finer than 2.5 $\mu$m (PM2.5) which results from exposure to residual oil fly ash and cigarette smoke (Magari et al., 2001; Cavallari et al., 2008). Harezlak et al. (2007) present an analysis of part of this data and found both negative and positive time-specific associations in the morning that corresponded to the workers’ break times. Our analysis focuses on the morning hours of the work day where Harezlak et al. (2007) saw the largest effects. We also make available MATLAB code for the implementation of our method at https://github.com/markjmeyer/WPHFLM.

The remainder of the paper is organized as follows: Section 2 provides a brief introduction to the discrete wavelet-packet transformation. Section 3 details the formulation of the the Bayesian wavelet-packet HFLM along with a discussion of inferential procedures in this modeling framework. Sections 4 and 5 present the results of our simulation study and the application of our model to the Journeyman data, respectively. Finally, in Section 6, we provide a discussion of the methodology.

2 Discrete Wavelet-packet Transformation

We begin with a brief description of the partial discrete wavelet transformation (DWT). Consider a 1-dimensional function, $x(v)$, which we discretely observe as \textbf{x}=\left[\begin{array}[]{ccc}x_{1}&\cdots&x_{V}\end{array}\right]^{\prime}, where $V=2^{k}$ for some positive integer value $k$. For a given mother wavelet, a partial DWT to $J=3$ levels results in a set of wavelet coefficients that can be further separated into approximation coefficients, a, and detail coefficients, d. The pyramid algorithm for performing the partial DWT is graphically depicted in Figure 1. Here, $\mathcal{G}$ and $\mathcal{H}$ denote the low and high pass filters for the corresponding mother wavelet. The resulting decomposition consists of the detail coefficients from each level plus the final approximation coefficients. Post-transformation, the wavelet-space representation of x is then \textbf{w}=\left[\begin{array}[]{cccc}\textbf{a}_{3}&\textbf{d}_{3}&\textbf{d}_{2}&\textbf{d}_{1}\end{array}\right]^{\prime}. The DWT can be expressed as a matrix multiplication by an orthogonal matrix, $W$, thus it can be shown that $\textbf{w}=\textbf{x}W$.

From Figure 1, we see the DWT applies the filters only to the successive sets of approximation coefficients. For a given mother wavelet, the partial discrete wavelet-packet transformation (DWPT) begins in the same way, passing first the low and then high-pass filters over the original signal. However, after the first level, the DWPT applies $\mathcal{G}$ and $\mathcal{H}$ to all coefficients at each level, reversing the order when applied to detail coefficients. Figure 2 illustrates the partial DWPT decomposition for $J=3$ levels. Once again, $\mathcal{G}$ and $\mathcal{H}$ denote the low and high pass filters of the corresponding mother wavelet. The resulting representation of x is then the final set of coefficients from the last level, \textbf{w}_{P}=\left[\begin{array}[]{cccccccc}\textbf{a}_{3,0}&\textbf{d}_{3,1}&\textbf{a}_{3,2}&\textbf{d}_{3,3}&\textbf{a}_{3,4}&\textbf{d}_{3,5}&\textbf{a}_{3,6}&\textbf{d}_{3,7}\end{array}\right]^{\prime}. As in the DWT, the DWPT can also be expressed as a matrix multiplication of an orthogonal matrix, $W_{P}$. Thus $\textbf{w}_{P}$ can be shown to be $\textbf{w}_{P}=\textbf{x}W_{P}$. For both the DWT and DWPT, the resulting set of wavelet-packet coefficients will have the same length as the original signal.

Similar to the DWT, the wavelet-packet coefficients comprising $\textbf{w}_{P}$ are indexed by a scale and location. The scale indexes the bin the coefficient is in at the final level of the decomposition and thus corresponds to the second subscript from the $j=3$ level in Figure 2 while the location denotes the position of the coefficient within the set of coefficients at a given scale. Since the number of elements in $\textbf{w}_{P}$ is the same as in x, within a scale, the time-ordering of the observed signal is preserved in the ordering of the coefficients. Suppose we have a second function, $y(t)$, that we discretely observe on a grid such that \textbf{y}=\left[\begin{array}[]{ccc}y_{1}&\cdots&y_{T}\end{array}\right]^{\prime}. If we perform the DWPT on y, we will also obtain a set of wavelet-packet coefficients that preserve the time-ordering of the original signal within each scale. Provided the elements of x are sampled in time concurrently or before the elements of y, we can use the location index from their respective DWPT decompositions to constrain the surface of estimation within each scale in the wavelet-packet space. When we apply the inverse DWPT (IDWPT), the constraint in then maintained in the data-space. For more details on DWPTs, see Percival and Walden (2000, chap. 6), Misiti et al. (2007), and Nason (2008, chap. 2).

3 Bayesian Historical Functional Linear Model

We begin with the model in (2) where we constrain the estimation to the region defined by $\{v,\in\mathcal{V},t\in\mathcal{T}:v\leq t\}$. We assume the within-function errors come from a Gaussian process. Thus, $E_{i}(t)\sim\mathcal{GP}\left(0,\Sigma_{E}\right)$, where $\Sigma_{E}$ is an unstructured covariance matrix. Because the data, $y_{i}(t)$ and $x_{i}(v)$, arrive sampled on a grid of equally spaced time points $t=[t_{1},\cdots,t_{T}]^{\prime}$ and $v=[v_{1},\cdots,v_{V}]^{\prime}$, we use the discrete version of the model: $\textbf{y}_{i}=\textbf{x}_{i}\boldsymbol{\beta}+\textbf{e}_{i}$ for the vectors $\textbf{y}_{i}$, $\textbf{x}_{i}$, and $\textbf{e}_{i}$ and matrix of coefficients $\boldsymbol{\beta}$. We recommend centering and scaling both the outcome and predictor functions first. Thus, without loss of generality, we drop the intercept function from the model. Stacking the response vectors and predictor vectors into matrices gives

[TABLE]

where for $N$ total curves, E and Y are $N\times T$ while X is $N\times V$. The constrained region of integration in Model (2) restricts the form of the functional regression coefficients so that $\beta(v_{k},t_{k^{\prime}})=0$ if $v_{k}>t_{k^{\prime}}$. If $T=V$ and $t_{1}=v_{1},t_{2}=v_{2},\ldots,t_{T}=v_{V}$, then the discrete version of $\boldsymbol{\beta}$ is an upper triangular matrix of the form

[TABLE]

with zeros below the main diagonal. We propose the use of wavelet-packets to enforce this constraint in the wavelet-packet space which, given that the time ordering is preserved, will ensure the constraint is maintained in the data-space as well.

3.1 Historical Constraint via Wavelet-Packets

Working with (3), we apply the DWPT separately to each row of Y and to each row of X. Performing this transformation is equivalent to the post-multiplication of the approximately orthonormal projection matrices resulting from the DWPT (Percival and Walden, 2000). The resulting decompositions have the form $\textbf{Y}=\textbf{Y}^{W_{P}}\textbf{W}_{P,Y}$ and $\textbf{X}=\textbf{X}^{W_{P}}\textbf{W}_{P,X}$ where $\textbf{W}_{P,Y}$ and $\textbf{W}_{P,X}$ are orthogonal matrices containing the wavelet packet basis functions. Then for the two dimensional decomposition on $\boldsymbol{\beta}=\textbf{W}_{P,X}^{\prime}\boldsymbol{\beta}^{W_{P}}\textbf{W}_{P,Y}$ Model (3) in the wavelet-packet space is $\textbf{Y}^{W_{P}}\textbf{W}_{P,Y}=\textbf{X}^{W_{P}}\textbf{W}_{P,X}\textbf{W}_{P,X}^{\prime}\boldsymbol{\beta}^{W_{P}}\textbf{W}_{P,Y}+\textbf{E}^{W_{P}}\textbf{W}_{P,Y}$ for $\textbf{E}=\textbf{E}^{W_{P}}\textbf{W}_{P,Y}$. Post-multiplying by $W_{P,Y}$ and recognizing the orthogonality of the wavelet-packet basis matrices, this model reduces to $\textbf{Y}^{W_{P}}=\textbf{X}^{W_{P}}\boldsymbol{\beta}^{W_{P}}+\textbf{E}^{W_{P}}$ and subject-specific model $\textbf{y}_{i}^{W_{P}}=\textbf{x}_{i}^{W_{P}}\boldsymbol{\beta}^{W_{P}}+\textbf{e}_{i}^{W_{P}}$. This model fits into the generalized basis expansion framework for function-on-function regression considered in Meyer et al. (2015). The key difference, however, lies in our novel use of wavelet-packets as a basis function to induce the historical constraint whereas Meyer et al. (2015) is concerned only with basis functions to estimate a full, unconstrained surface or set of surfaces.

We enforce the constraint in the wavelet-packet space via our prior specification on the elements of $\boldsymbol{\beta}$. Let the DWP transformations be indexed by scales $j=1,\ldots,J^{y}$ and $s=1,\ldots,S^{x}$ and locations $k=1,\ldots,K_{j}^{y}$ and $\ell=1,\ldots,L_{s}^{x}$ in the Y and X wavelet-packet spaces respectively. Consistent with previous work on wavelet-based models in function regression, we place spike-and-slab priors on model coefficients. To restrict the surface in wavelet-packet space, we set coefficients where $\ell>k$ to zero. Thus our prior on the elements of $\boldsymbol{\beta}^{W_{P}}=\left[\beta^{W_{P}}_{s\ell,jk}\right]$ is $\beta_{s\ell,jk}^{W_{P}}\sim\ 1(\ell\leq k)\gamma_{s\ell,jk}N(0,\tau_{s\ell,j\cdot})+[1-\gamma_{s\ell,jk}]d_{0}$, where $\gamma_{s\ell,jk}\sim Bern(\pi_{s\ell,j\cdot})$ and $d_{0}$ is a point-mass distribution at zero. The regularization parameters, $\tau_{s\ell,j\cdot}$ and $\pi_{s\ell,j\cdot}$, smooth over locations $k$ which we denote using the “dot” notation in the subscript. We assume inverse-gamma and beta hyper priors, respectively, for the regularization parameters with hyper-parameters fixed and based on the data.

Morris and Carroll (2006) show that after a wavelet transformation, assuming independence in the wavelet-space does not imply independence in the data-space and therefore wavelets accommodate a wide range of covariances in the data-space. As wavelet-packets share the same whitening properties of wavelets, we assume independence in the wavelet-packet space (Percival and Walden, 2000). Thus we assume $\textbf{e}_{i}^{W_{p}}\sim N(0,\boldsymbol{\Sigma}^{W_{p}})$ where $\boldsymbol{\Sigma}^{W_{p}}=\text{diag}\left\{\sigma^{2}_{jk}\right\}$, which varies by the scale and location of the Y wavelet-packet coefficients. We place an inverse gamma prior on $\sigma^{2}_{jk}$. The independence assumption allows us to sample the coefficients corresponding to different $j$ and $k$ combinations separately.

Using the prior specifications, we now describe our sampling algorithm. For the $jk$th wavelet-packet space coefficient from the $\textbf{Y}^{W_{P}}$ decomposition and the $s\ell$th column of $\textbf{X}^{W_{P}}$, the conditional posterior distribution is a mixture of a point-mass at zero and a normal distribution of the form

[TABLE]

where $\mu_{s\ell,jk}=\hat{\beta}_{s\ell,jk}^{W_{P}}(1+\Lambda_{s\ell,jk}/\tau_{s\ell,j})^{-1}$ and $\epsilon_{s\ell,jk}=\Lambda_{s\ell,jk}(1+\Lambda_{s\ell,jk}/\tau_{s\ell,j})^{-1}$ for the OLS and variance estimates $\hat{\beta}_{s\ell,jk}^{W_{P}}$ and $\Lambda_{s\ell,jk}$ at the current step. The conditional for $\gamma_{s\ell,jk}$ is

[TABLE]

were $\alpha_{s\ell,jk}=O_{s\ell,jk}/\left(O_{s\ell,jk}+1\right)$ for $O_{s\ell,jk}=\pi_{s\ell,j}/(1-\pi_{s\ell,j})\text{BF}_{s\ell,jk}$, $\text{BF}_{s\ell,jk}=\left(1+\frac{\tau_{s\ell,j}}{\Lambda_{s\ell,jk}}\right)^{-1/2}$ $\exp{\left\{\frac{1}{2}\zeta^{2}_{s\ell,jk}\left(1+\frac{\Lambda_{s\ell,jk}}{\tau_{s\ell,j}}\right)\right\}}$, and $\zeta_{s\ell,jk}$ equal to the ratio of the current values of $\beta^{W_{P}}_{s\ell,jk}$ to the current estimate of the standard deviation of $\beta^{W_{P}}_{s\ell,jk}$. Through the indicator function, $1(\ell\leq k)$, we enforce the historical constraint by forcing coefficients for which $\ell>k$ to come from the point-mass density $d_{0}$ and thus be set to zero.

The conditionals for the diagonal elements of the wavelet-packet space variance components, $\sigma^{2}_{jk}$, have the form

[TABLE]

where $\pi\left(\sigma^{2}_{jk}\right)$ is the prior density on $\sigma^{2}_{jk}$ which we take to be inverse gamma with parameters $a_{\sigma^{2}}$ and $b_{\sigma^{2}}$ both set to the empirical Bayes estimates. To ensure the variance components are not too close to zero, we employ a Metropolis-Hastings step to sample them. The proposal densities are independent Gaussians, truncated at zero and centered at the previous value in each chain. The full conditionals for the regulation parameters are

[TABLE]

with $a_{\tau}$, $b_{\tau}$, $a_{\pi}$, and $b_{\pi}$ set to the empirical Bayes estimates. The prior specifications are consistent with the previous work on wavelet-based functional models including Meyer et al. (2015), Malloy et al. (2010), and Morris and Carroll (2006). For more details on the empirical Bayes estimates, see Morris and Carroll (2006). Our sampler then iterates between draws from (5) to (9) until convergence. Upon completion of the algorithm, we apply the inverse DWPT to the posterior samples of $\boldsymbol{\beta}^{W_{P}}$ to obtain estimates in the data space of $\boldsymbol{\beta}$, the upper triangular matrix of historically constrained coefficients.

3.2 Thresholding and Wavelet Details

Our sampler becomes computationally intensive as $V$ increases with computation time increasing almost linearly. In previous work on wavelet-based FFR, Meyer et al. (2015) address this issue by reducing the dimension of the wavelet transformed X. In their formulation, the authors use wavelet-Principal Components (wPC) retaining columns containing a large amount of the variability in X. The wPC decomposition involves first performing a DWT on X and then performing a singular value decomposition (SVD). While such an approach reduces computation time and achieves additional denoising, it does not work in the historical framework, as performing an SVD on $\textbf{X}^{W_{P}}$ would break the time ordering maintained by the transformation. To remedy the computational concerns and simultaneously achieve additional denoising, we propose a simplified thresholding procedure where we threshold all coefficients in the larger scale levels to zero. Thus, we consider keeping the scales that comprise either 25% or 50% of the wavelet-packet coefficients. For example, using the $J=3$ level partial DWPT, the first approach would result in retaining the coefficients $\textbf{a}_{3,0}$ and $\textbf{d}_{3,1}$ in Figure 2. The second approach, keeping 50%, would correspond to retaining the coefficients $\textbf{a}_{3,0}$, $\textbf{d}_{3,1}$, $\textbf{a}_{3,2}$, and $\textbf{d}_{3,3}$, again from Figure 2. Our simulation setting compares the performance of both approaches in terms of computational efficiency and estimation. It is important to note that we do not threshold the outcome and retain all coefficients of $\textbf{Y}^{W_{P}}$.

The choice of mother wavelet will depend largely on the data context. In our application and simulation, we use Daubechies wavelets with 3 vanishing moments for the separate DWPTs on both Y and X. Wavelet-packets also require a choice of boundary padding, we select zero padding—padding the boundary with zeros. Finally, we must choice the number of levels to take our partial DWPT. Because of our reliance on the maintained time ordering within scales but between decompositions, we use the same number of levels of decomposition for both Y and X. In simulation, we use $J=3$ levels for sampling densities of $T=V=64$ and $128$. For denser sampling rates, we recommend increasing $J$ and potentially adjusting the thresholding.

3.3 Posterior Functional Inference

Previous implementations of the HFLM do not discuss inferential procedures instead, in some cases, developing measures of model fit (Malfait and Ramsay, 2003; Harezlak et al., 2007). In the wavelet-based functional literature, Meyer et al. (2015) propose the use of a Bayesian False Discovery Rate (BFDR) procedure also used by Morris et al. (2008) and Malloy et al. (2010), as well as joint credible bands as discussed in Ruppert, Wand, and Carroll (2003). We consider the use of both methods in the context of Bayesian wavelet-packet HFLM.

The BFDR procedure utilizes the MCMC samples to estimate the posterior probability of a given coefficient being greater than a $\delta$-fold intensity change. Once these values are determined they are ranked and a cut-off selected to control the overall FDR at a pre-specified global $\alpha$-bound. Suppose we have $M$ MCMC samples and $\beta^{(m)}(v,t)$ is $m^{\text{th}}$ draw from the posterior estimated surface. Then for $\left\{v\in\mathcal{V}\text{ and }t\in\mathcal{T}\text{ s.t. }v\leq t\right\}$, we find the probability $P_{B}(v,t)=Pr\left\{\left|\beta(v,t)\right|>\delta|y\right\}\approx\frac{1}{M}\sum_{m=1}^{M}1\left\{\left|\beta^{(m)}(v,t)\right|>\delta\right\}.$ We then flag the set of coefficients on the historical surface satisfying $\psi=\left\{(v,t):v\leq t\text{ and }P_{B}(v,t)\geq\phi_{\alpha}\right\}$. Here, $\phi_{\alpha}$ is determined by first ranking the values of $P_{B}$ in descending order across all coefficients to obtain the set $\left\{P_{(r)},r=1,\ldots,R\right\}$ where $R$ is the total number of coefficients satisfying the historical constraint. We then define the cut-off value as \lambda=\max\big{[}r^{*}:\frac{1}{r^{*}}\sum_{r=1}^{r^{*}}\left\{1-P_{(r)}\right\}\leq\alpha\big{]}. Thus, we select coefficients with $P_{B}$ great than or equal to $\phi_{\alpha}=P_{(\lambda)}$ as significant.

For interval estimation, we consider both point-wise credible intervals as well as joint credible intervals. Point-wise credible intervals are constructed by taking the $\alpha/2$ and $1-\alpha/2$ quantiles of the posterior samples taken at each coefficient for some choice of $\alpha$. As in Meyer et al. (2015), we construct joint credible intervals using $I_{\alpha}(v,t)=\hat{\beta}(v,t)\pm q_{(1-\alpha)}\left[\widehat{\text{St.Dev}}\left\{\hat{\beta}(v,t)\right\}\right]$, where $\hat{\beta}(v,t)$ and $\widehat{\text{St.Dev}}\left\{\hat{\beta}(v,t)\right\}$ are the mean and standard deviation of the posterior samples, respectively, and $q_{(1-\alpha)}$ is the $(1-\alpha)$ quantile taken over all posterior samples of the quantity

[TABLE]

Such an interval satisfies $Pr\left\{L(v,t)\leq\beta(v,t)\leq U(v,t)\ \forall\ v\in\mathcal{V},t\in\mathcal{T}\text{ s.t. }v\leq t\right\}\geq 1-\alpha$, where $L(v,t)$ and $U(v,t)$ are the corresponding upper and lower interval bounds. This procedure yields a joint $100(1-\alpha)$ interval for the historical association surface.

4 Simulation Study

We consider four different historical surface scenarios representing plausible relationships between $x(v)$ and $y(t)$: a lagged effect of $x(v)$ on $y(t)$, a cumulative effect, a time-specific effect, and a delayed time-specific effect. Figure 3 displays the surfaces. We present their mathematical expressions in the Supplementary Material. For each scenario, we assume the data is sampled on a grid between $t_{1}=v_{1}=1$ and $t_{T}=v_{V}=64$ and vary the sampling rate such that $T=V=64$ or $T=V=128$. We also vary the sample size, considering $N=50$ and $N=200$. For the sparser, $T=64$ setting, we vary the percent of the columns of $\textbf{X}^{W_{P}}$ we retain, first retaining 25% of coefficients and then retaining 50%. In the $T=128$ setting, we only retain 25% of coefficients since retaining more considerably increases the computational burden. Thus in total, we present the results from 24 different simulation settings.

To evaluate each setting, we calculate the root mean integrated squared error (RMISE) and determine coverage probabilities for both point-wise and joint credible intervals. We also graphically explore the properties of the BFDR procedure. For each “true” historical surface $\boldsymbol{\beta}$, we generate $N$ $\textbf{x}_{i}$ curves from a mean zero Gaussian Process with a first-order auto-regressive (AR1) covariance structure. We base the variance and correlation parameters of the AR1 covariance off of the PM${}_{2}.5$ data from the Journeyman data, setting $\sigma^{2}_{AR,X}=3.5$ and $\rho_{X}=0.75$. Next we generate within-subject error functions, $\textbf{e}_{i}$, from a separate mean zero Gaussian Process with an AR1 covariance structure. Once again we base the parameters of the covariance matrix off of the Journeyman data, $\sigma^{2}_{AR,E}=0.1$ and $\rho_{E}=0.5$. We then simulate the outcome functions using $\textbf{y}_{i}=\textbf{x}_{i}\boldsymbol{\beta}+\textbf{e}_{i}$. For each setting, we repeat this data generation process 200 times and obtain 2000 posterior samples, discarding the first 1000, for each simulated dataset. For each simulated dataset we perform a $J=3$ level DWPT on both the outcome and predictor using Daubechies wavelets with 3 vanishing moments. All computation is done using MATLAB version R2017a on a desktop with a 3.2 GHz Intel Core i5 processor and 16 GB of memory.

Table 1 presents RMISE averaged across all 200 simulated datasets for each setting under consideration while Figure 4 displays a single estimated surface with near average RMISE for the case when $N=50$, $T=64$ and only 25% of the $\textbf{X}^{W_{P}}$ coefficients are retained. From Table 1, we see that increasing either the sample size or the percent of retained coefficients tends to decrease RMISE, although all settings produce similar results. Denser sampling also decreases the RMISE. Regardless of sample size or sampling density, retaining 25% of the $X$-space wavelet-packet coefficients preserves roughly 75% of the energy while retaining 50% of coefficients preserves 88% of the energy in $X(v)$. To calculate average preserved energy, we find it first for each simulated subject and then average across all subjects. From Figure 4, we see that a single simulated dataset with near average RMISE accurately estimates the true relationship, even when only retaining an average of 75% of the energy in the predictor function.

We assess inference using Table 2 and Figure 5. Table 2 contains mean point-wise coverage probabilities, calculated by first averaging over the historical surface and then over all 200 simulated datasets. We calculate coverage for both the point-wise credible interval and the joint credible interval. From Table 2, we see that in general the joint credible intervals provide higher coverage for each setting. Coverage also tends to increase, regardless of interval type, as the percent of retained coefficients increase, though in some settings the differences are negligible. For some settings, coverage of the joint intervals is affected by sampling density and sample size. However, for the joint intervals, coverage is above the nominal level for most settings. This is not the case for the point-wise intervals where coverage sits below nominal level, in some bases considerably below. All intervals are at the 95% level.

Figure 5 shows the local BFDR values using a cut-off of $\delta=0.5$ which is equal to half the largest signal, in absolute value, from each surface. Each plotted value corresponds to the posterior probability that $|\beta(v,t)|>0.5$ for each surface from Figure 4, that is for one simulated dataset with near average RMISE. We see that using this $\delta$-intensity change, the BFDR highlights the features of each setting with elevated (or depressed) association. In the Supplementary Material, we consider one other value of $\delta$, which, while smaller, produces similar results. Choice of a specific $\delta$ is data-dependent, but several different values should be considered.

Figures 4 and 5 present the results from the simulation for $N=50$, $T=64$, and retaining 25% of the columns of $\textbf{X}^{W_{P}}$. Graphical results from simulated datasets with near average RMISE for the other settings produce similar results. Consequently, we present limited additional BFDR results in the Supplementary Material. Specifically, we only show BFDR graphs for all remaining 25% retained coefficients and $T=64$ settings for $\delta=0.25$ and 0.50.

5 Analysis of Journeyman Boilermaker Data

Harezlak et al. (2007) and Cavallari et al. (2008) analyze data from 14 Journeyman boilermakers relating their five-minute SDNN to five-minute PM2.5 levels over the course of the workday. Exposure to PM2.5 came from two sources: residual oil fly-ash, a by-product of the manufacturing process, and cigarette smoke during mandatory breaks. Harezlak et al. (2007) found elevated levels association during the morning hours that corresponded to mandatory breaks which were then followed by depressions in the surface estimated over the trapezoidal region. To try to better estimate these spikes in association, we focus our analysis on the first three hours of the workday beginning at 8:30am and going until 11:30am. Further, we estimate joint intervals and calculate the BFDR to determine the significance of these peaks and troughs.

Prior to analysis, we log transform and center and scale both SDNN and PM2.5. Thus changes in the estimated surface correspond to one unit changes in $\log($PM${}_{2.5})$ and result in changes in $\log$(SDNN). To investigate the morning hours only, we take $T=V=34$ measurements. As in the simulation, we take $J=3$ levels of decomposition for both wavelet-packet transformations and use Daubechies wavelets with 3 vanishing moments. Given the results of the simulation, we retain only the first two levels of wavelet-packet coefficients, resulting in only 25% of the columns of $\textbf{X}^{W_{P}}$ being kept. We let this model run longer than the simulated settings, taking 2000 posterior samples after a burnin of 2000. We monitor convergence with running mean and trace plots to ensure all chains converged (see Supplementary Material).

Figure 6 contains the posterior estimated historical surface for the Journeyman data whereas Figure 7 presents the lower and upper bounds of the joint credible intervals. From Figure 6, moving across the $v$-axis from left to right, we see a pattern of time-specific depressions followed by elevations, which culminate in a large peak around $v=8$ and $t=20$. This suggests that an exposure 40 minutes into the day resulted in elevated SDNN that was sustained until, and peaked around, 100 minutes. After this peak, there is another depression that occurs around $v=16$ and $t=20$. These patterns are consistent with the analysis in Harezlak et al. (2007).

The joint bands experience edge effects along the $t$-axis, which somewhat distort the scale of the intervals in Figure 7. However, from this figure, we are able to see the elevation of the lower bound around the major peak around $v=8$ and $t=20$ as well as the depression around $v=16$ and $t=20$. The ridge that runs along $v=3$ in Figure 6 also appears to have an elevated lower bound. In the Supplementary Material, we present the point-wise intervals which are smaller in width, but exhibit similar features to the joint intervals. The intervals suggest regions of both increased and decreased association between SDNN and PM2.5 throughout the morning. To more accurately pinpoint these regions, we turn to the BFDR procedure.

Figure 8 presents the local BFDR values, $P_{B}(v,t)$, for each coefficient (top row) as well as the flagged coefficients from the set $\psi_{\alpha}$ (bottom row) using $\alpha=0.05$. The columns correspond to two different $\delta$ intensity cut-offs representing corresponding changes in $\log$(SDNN) of 0.5 and 0.75. The lower cut-off value highlights significant coefficients in the regions where we see the strongest effect sizes from Figure 6. These regions do not extend from the diagonal nor to the edge of the surface, suggesting a series of delayed time-specific effects that result from exposures within the first 100 minutes of the workday. Some of these effects are strong, as evidenced by the sustained significance at the higher $\delta$ level. In the Supplementary Material, we present additional results which shows the main peak around $v=8$ and $t=20$ is significant at an intensity change of $\delta=1.00$.

6 Discussion

One of the difficulties of estimating a historical functional effect is in maintaining the historical constraint. Previous authors have used tent shaped basis functions to achieve the constraint. When using wavelets coefficients from a DWT as a basis function, constraining the surface is difficult due to the lack of a clear relationship between the wavelet-space coefficients and the time domain. However, wavelet-packet coefficients from a DWPT have an exploitable relationship between their location indices and the original time scale that allows us to constrain the surface in the wavelet-packet space and maintain that constraint when projecting back to the data-space. In this work, we show that when performing the DWPT on both $y(t)$ and $x(v)$, we can use wavelet-packet basis functions, for a given choice of mother wavelet, to model historical effects.

The current literature on HFLMs is limited, largely focusing on the estimation of surface effects and the determination of model fit criterion. Most methods implement spline-based basis functions and, while one method does allow for different basis expansions, the authors only present results for spline-based models. Further, inference procedures are not discussed in the existing literature. To the best of our knowledge, the model we present here represents the first work in wavelet-based modeling of historical effects as well as the first Bayesian HFLM. Additionally, our method employs a novel use of wavelet-packets, which have not previously been used in functional regression models, to constrain the surface of estimation. Finally, we perform inference using two established multiplicity adjusted Bayesian inferential procedures for describing the uncertainty in the surface estimate through the use of joint credible intervals and identifying significant regions of coefficients using the BFDR.

We demonstrate in simulation that the wavelet-packet HFLM can accurately estimate several realistic historical surface settings. In particular, we show that regardless of the percent of $\textbf{X}^{W_{P}}$ coefficients we retain, the wavelet-packet HLFM has good RMISE levels that are similar across sample size and sampling density. Further, the joint credible intervals used in Meyer et al. (2015) for the wavelet-based FFR outperform point-wise credible intervals in terms of average coverage. And while their coverage is larger than the nominal, we prefer this to the alternative. Thus we suggest the use of the joint intervals for quantifying uncertainty in the wavelet-packet HFLM. To identify regions of significant coefficients, we implement the BFDR to show that for two reasonable choices of $\delta$ relative to the data and max signal, our model can identify regions of association in the surface. Finally, we apply the proposed model to analyze data on the association between HRV and PM2.5 exposure in a panel of journeymen, focusing on the first three hours of exposure during the workday. Using the wavelet-packet HFLM, we are able to not only estimate regions of association but clearly identify them as representing significant changes in SDNN using the BFDR.

SUPPLEMENTARY MATERIAL

Supplementary Material to Bayesian HFLM with Wavelet-packets:

Contains additional results from the simulation and application. (pdf)

MATLAB Code:

Available at https://github.com/markjmeyer/WPHFLM.