Sparse Estimation of Historical Functional Linear Models with a Nested Group Bridge Approach

TL;DR

This paper introduces a novel nested group bridge approach for estimating and identifying the relevant recent past in historical functional linear models, improving lag detection and coefficient estimation.

Contribution

It develops a new estimation procedure using finite element methods and a nested group bridge penalty to simultaneously estimate the coefficient function and determine the historical lag.

Findings

Effective lag detection demonstrated in simulations.

Improved estimation accuracy over traditional methods.

Application to EMG data shows practical utility.

Abstract

The conventional historical functional linear model relates the current value of the functional response at time t to all past values of the functional covariate up to time t. Motivated by situations where it is more reasonable to assume that only recent, instead of all, past values of the functional covariate have an impact on the functional response, we investigate in this work the historical functional linear model with an unknown forward time lag into the history. Besides the common goal of estimating the bivariate regression coefficient function, we also aim to identify the historical time lag from the data, which is important in many applications. Tailored for this purpose, we propose an estimation procedure adopting the finite element method to conform naturally to the trapezoidal domain of the bivariate coefficient function. A nested group bridge penalty is developed to provide simultaneous estimation of the bivariate coefficient function and the historical lag. The method is demonstrated in a real data example investigating the effect of muscle activation recorded via the noninvasive electromyography (EMG) method on lip acceleration during speech production. The finite sample performance of our proposed method is examined via simulation studies in comparison with the conventional method.

Tables

Table 1. Table 1: Summary of the 100 simulation replications for each scenario, including the root mean squared error (RMSE), percent bias (%bias) and standard deviation (SD) of δ ^ ^ 𝛿 \widehat{\delta} , and the square root of the mean integrated squared errors (RISE) of β ^ ​ ( s , t ) ^ 𝛽 𝑠 𝑡 \widehat{\beta}(s,t) along with the corresponding standard deviation in parentheses. The two methods implemented are the proposed nested group bridge approach (NGB) and the penalized linear model approximation approach (PLMA).

	RMSE of $\hat{\delta }$			%bias of $\hat{\delta }$			SD of $\hat{\delta }$			RISE of $\hat{\beta }(s,t)$
	NGB	PLMA		NGB	PLMA		NGB	PLMA		NGB	PLMA
Scenario 1	0.0173	0.0206		1.2	1.7		0.0162	0.0188		1.29 (0.25)	0.83 (0.35)
Scenario 2	0.0480	0.1587		-5.4	-31.5		0.0396	0.0192		69.86 (7.04)	69.33 (7.54)
Scenario 3	0.0548	0.2426		-9.2	-48.3		0.0297	0.0235		8.74 (1.33)	9.75 (1.34)

Equations

$$y_i(t)=\alpha (t)+{\int }_0^T\beta (s,t)x_i(s)ds+{\epsilon }_i(t),t\in [0,T],$$

$$y_i(t)=\alpha (t)+{\int }_{t-\delta }^t\beta (s,t)x_i(s)ds+{\epsilon }_i(t),t\in [0,T],$$

$$y_i(t)=\alpha (t)+{\int }_0^t{x}_i(s)\beta (s,t)ds+{\epsilon }_i(t),\text{mbox}fort\in [0,T].$$

$$y_i^{\ast }(t)={\int }_0^t{x}_i^{\ast }(s)\beta (s,t)ds+{\epsilon }_i^{\ast }(t),t\in [0,T].$$

$$\hat{\alpha }(t)=\overline{y}(t)-{\int }_0^t\overline{x}(s)\hat{\beta }(s,t)ds.$$

$$y_i(t)={\int }_0^t{x}_i(s)\beta (s,t)ds+{\epsilon }_i(t),t\in [0,T].$$

$$\beta (s,t)\approx \sum _{k=1}^K{b}_k{\varphi }_k(s,t).$$

$$y_i(t)$$

$$\frac{1}{N}{\int }_0^T\sum _{i=1}^N{\left\{y_i(t)-\sum _{k=1}^K{b}_k{\psi }_{ik}(t)\right\}}^{2}dt+\lambda \sum _{j=1}^{M+1}{c}_j||{\mathbf{b}}_{A_j}|{|}_1^{\gamma }+{\mathbf{b}}^T\mathbf{R}\mathbf{b},$$

$$\mathbf{R}={\omega }_H{{\mathbf{D}}_H}^T{\mathbf{D}}_H+{\omega }_V{{\mathbf{D}}_V}^T{\mathbf{D}}_V+{\omega }_P{{\mathbf{D}}_P}^T{\mathbf{D}}_P,$$

$$\frac{T}{Q}\sum _{q=1}^Q\sum _{i=1}^N{\left\{y_i(t_q)-\sum _{k=1}^K{b}_k{\psi }_{ik}(t_q)\right\}}^2.$$

$$\frac{T}{Q}{\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right)}^T\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right).$$

$$\frac{1}{N}{\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right)}^T\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right)+\lambda \sum _{j=1}^{M+1}{c}_j||{\mathbf{b}}_{A_j}|{|}_1^{\gamma }+{\mathbf{b}}^T\mathbf{R}\mathbf{b},$$

$$\frac{1}{N}{\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right)}^T\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right)+\sum _{j=1}^{M+1}{\theta }_j^{1-1/\gamma }{c}_j^{1/\gamma }||{\mathbf{b}}_{A_j}|{|}_1+\tau \sum _{j=1}^{M+1}{\theta }_j+{\mathbf{b}}^T\mathbf{R}\mathbf{b},$$

$$s.t.{\theta }_j\ge 0,j=1,\dots ,M+1,$$

$$\frac{1}{N}{\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right)}^T\left(\mathbf{y}-\mathbf{\Psi }\mathbf{b}\right)+\sum _{k=1}^K{g}_k|b_k|+\tau \sum _{j=1}^{M+1}{\theta }_j+{\mathbf{b}}^T\mathbf{R}\mathbf{b},$$

$$s.t.{\theta }_j\ge 0,j=1,\dots ,M+1,$$

$$\frac{1}{N}\left\{{\left({\mathbf{y}}^{\ast }-{\mathbf{\Psi }}^{\ast }{\mathbf{b}}^{\ast }\right)}^T\left({\mathbf{y}}^{\ast }-{\mathbf{\Psi }}^{\ast }{\mathbf{b}}^{\ast }\right)+\sum _{k=1}^K|b_k^{\ast }|\right\}+\tau \sum _{j=1}^{M+1}{\theta }_j,$$

$$s.t.{\theta }_j\ge 0,j=1,\dots ,M+1,$$

$${\theta }_j^{(s)}=c_j{\left(\frac{1-\gamma }{\tau \gamma }\right)}^{\gamma }||{\mathbf{b}}_{A_j}^{(s-1)}|{|}_1^{\gamma },j=1,\dots ,M+1,$$

$$g_k^{(s)}=\sum _{j=1}^{\ell (k)}{\left({\theta }_j^{(s)}\right)}^{1-1/\gamma }{c}_j^{1/\gamma },k=1,\dots ,K,$$

$${\mathbf{G}}^{(s)}=N^{-1}\text{mbox}diag\left(1/g_1^{(s)},\dots ,1/g_K^{(s)}\right),$$

$${\mathbf{\Psi }}^{\ast (s)}=[{\mathbf{\Psi }}^T,\sqrt{{\omega }_{H}N}{\mathbf{D}}_H^T,\sqrt{{\omega }_{V}N}{\mathbf{D}}_V^T,\sqrt{{\omega }_{P}N}{\mathbf{D}}_P^T{]}^T{\mathbf{G}}^{(s)}.$$

$${\mathbf{b}}^{(s)}={\mathbf{G}}^{(s)}\text{mbox}{argmin}_{{\mathbf{b}}^{\ast }}\left\{{\left(\mathbf{y}-{\mathbf{\Psi }}^{\ast (s)}{\mathbf{b}}^{\ast }\right)}^T\left(\mathbf{y}-{\mathbf{\Psi }}^{\ast (s)}{\mathbf{b}}^{\ast }\right)+\sum _{k=1}^K|b_k^{\ast }|\right\}.$$

$${\hat{\beta }}_n(s,t)=\sum _{k=1}^K{\hat{b}}_{n,k}{\varphi }_k(s,t),{\hat{\delta }}_n=\frac{T}{M}\underset{j}{\min }\{1\le j\le M:{\hat{\beta }}_n(s,t)=0\text{mbox}on{D}_j\}.$$

$$df(\lambda ,\mathbf{\kappa })=\text{mbox}trace({\mathbf{\Psi }}_s({\mathbf{\Psi }}_s^T{\mathbf{\Psi }}_s+N{\mathbf{R}}_s{)}^{-1}{\mathbf{\Psi }}_s^T),$$

$$BIC(\lambda ,\mathbf{\kappa })=N\log (||\mathbf{y}-\mathbf{\Psi }{\hat{\mathbf{b}}}_n(\lambda ,\mathbf{\kappa })|{|}_2^2/N)+\log (N)df(\lambda ,\mathbf{\kappa }),$$

Full text

Key words: Finite element; Functional data analysis; Function-on-function regression; Historical linear model; Nested group bridge penalty.

1 Introduction

We are interested in functional linear regression models, where a functional predictor $x_{i}(t)$ is used to explain a functional response $y_{i}(t)$ over a time interval $[0,T]$, $i=1,\ldots,N$. There are a number of possible ways to build up the model. For instance, Ramsay and Dalzell (1991) introduced the functional linear model

[TABLE]

where $\alpha(t)$ is the intercept function, $\varepsilon_{i}(t)$ is the residual and $\beta(s,t)$ is an unconstrained regression bivariate coefficient function of primary interest. The bivariate coefficient function $\beta(s,t)$ represents the effect of the functional predictor at time $s$ on the functional response at time $t$. According to model (1), future facts about the covariate $x(s)\mbox{ with }s>t$ is used to explain $y(t)$, the response at current time. This is reasonable if the underlying process is periodic, but illogical otherwise. If the prediction of the response $y(t)$ only depends on concurrently observed predictor $x(t)$, a restriction is then introduced on $\beta(s,t)$, leading to the functional concurrent model $y_{i}(t)=\alpha(t)+x_{i}(t)\beta(t)+\varepsilon_{i}(t)$, for $t\in[0,T]$ where a univariate coefficient function $\beta(t)$ suffices to fully describe their relationship (Ramsay and Silverman, 2005). In other situations, the response at a particular time point may depend on the recent history of the predictor, then the functional concurrent model becomes too constrained. For example, the recovery of a patient may well depend on treatment received over the past a few days.

To properly model the relationship between the historical predictor and the current response, Malfait and Ramsay (2003) introduced the historical functional linear model

[TABLE]

where $\delta\in[0,T]$ represents a forward time lag back into the history such that $x_{i}(s)$ is likely to have an impact on $y_{i}(t)$ only over the time frame $[t-\delta,t]$ rather than over the complete time frame $[0,T]$. Notice that the degenerated case with $\delta=0$ reduces to a concurrent model. The time lag $\delta$ is typically unknown and of interest. The objective of our study is to estimate the time lag $\delta$ along with the coefficient function $\beta(s,t)$ from the data .

Concerning the classic function-on-function regression model (1), one can apply basis expansion to $x_{i}(t)$ and $y_{i}(t)$, and obtain a weighted least squares type estimator, possibly with a bivariate roughness penalty (Besse and Cardot, 1996; Ramsay and Silverman, 2005). Alternatively, one can compute the functional principal component scores for $x_{i}(t)$ and $y_{i}(t)$, and base the estimation on the functional principal component scores (Chen and Wang, 2011; Ivanescu, et al., 2015; Yao, Müller, and Wang, 2005).

The functional concurrent model, as a type of varying coefficient model, has also been studied relatively well in the literature (see, Fan and Zhang, 2000; Hastie and Tibshirani, 1993; Wu, Chiang and Hoover, 1998; Wu and Liang, 2004; Zhou, Huang, and Carroll, 2008). For a comprehensive introduction to functional regression models, we refer to monographs by Ferraty and Vieu (2006), Hsing and Eubank (2015), Kokoszka and Reimherr (2017), Ramsay and Silverman (2005), as well as the review papers by Morris (2015) and Wang, Chiou and Müller (2016) and references therein.

There also exist intermediate models between the functional linear model (1) and the functional concurrent model which try to model the effect of past values of the predictor on current response. A class of time-varying functional regression models were discussed by Müller and Zhang (2005), Şentürk and Müller (2008; 2010), etc. Brockhaus et al. (2017) and Greven and Scheipl (2017) discussed a general framework of functional regression, where many aforementioned models are included as special cases, and proposed gradient-boosting-based estimation procedures. Kim, Şentürk and Li (2011) proposed an estimation procedure that was geared towards sparse longitudinal data with irregular observation times and small amount of measurements per subject. Assuming the dependence of $y(t)$ on $x(s)$ for $s\in[t-\delta,t]$ does not change over time, Asencio, Hooker and Gao (2014) introduced the functional convolution model which is a functional extension of distributed lag models in time series, and proposed a penalized ordinary least squares estimator for the regression coefficient given the historical time lag $\delta$.

The historical functional linear model (2) is identifiable but not estimable, with effectively an infinite number of covariates, therefore regularization or roughness penalty on $\beta(s,t)$ is necessary in the estimation process (Ramsay and Silverman, 2005). Malfait and Ramsay (2003) regularized the fit by approximating $\beta(s,t)$ using an expansion with a finite number of basis functions. Such approach bears the well known limitation that a small number of basis functions would decrease the goodness of fit while a large number of basis functions would lead to unstable estimation. Harezlak et al. (2007) improved the fit of model (2) by imposing a discrete roughness penalty forcing neighboring coefficients to be similar, which is an extension of the P-spline by Eilers and Marx (1996). Both work focus on the estimation of coefficient function $\beta(s,t)$ for a given historical lag $\delta$, while the lag $\delta$ is decided in a separate manner. Malfait and Ramsay (2003) considered deciding $\delta$ as a model selection problem. Briefly put, they divided the whole time course $[0,T]$ into $M$ sub-intervals of equal length, fitted the model with $\delta=k/M$ for $k=0,\ldots,M$, and chose the best model from the $M+1$ candidate models according to a certain criterion.

In this work, we focus on the problem of estimating the historical lag and propose a tailored procedure for the simultaneous estimation of the coefficient function $\beta(s,t)$ and the lag $\delta$ in model (2). Considering the following “full” model,

[TABLE]

Figure 1 illustrates a scenario when the bivariate coefficient function $\beta(s,t)$ becomes zero over the triangular region defined by vertices $(0,\delta),(0,T),(T-\delta,T)$. This scenario corresponds to that $x(s)$ does not affect $y(t)$ for $s<t-\delta$. In other words, the full model (3) is equivalent to the historical functional linear model (2) when the support domain of the coefficient function $\beta(s,t)$ is a trapezoidal area $S$ defined by vertices $(0,0),(T,T),(T-\delta,T),(0,\delta)$, given that the lag $\delta>0$.

We propose a nested group bridge shrinkage method to estimate the historical lag $\delta$ and the coefficient function $\beta(s,t)$ with sparsity, which tackles the problem from a completely different perspective. Our method has two features. First of all, we utilize the triangular basis functions from the finite element method theory, which respects the non-rectangular support domain of the coefficient function $\beta(s,t)$. Such basis was also used by Malfait and Ramsay (2003) and Harezlak et al. (2007), as well as in spatial data analysis, where data distribute over irregularly shaped spatial domains with features like complex boundaries, strong concavities and interior holes (e.g. Ramsay, 2002; Sangalli, Ramsay and Ramsay, 2013). Then under model (3), we organize the basis coefficients in such a way that the nested group bridge penalty is able to shrink specifically the coefficient function $\beta(s,t)$ over the upper triangular region with vertices $(0,T)$, $(0,\delta)$ and $(T-\delta,T)$ towards zero. The group bridge shrinkage was originally proposed by Huang et al. (2009) for variable selection. Such penalty has been utilized by Wang and Kai (2015) for locally sparse estimation in nonparametric regression for a scalar-on-function historical functional linear model with B-spline basis function expansion. The major advantage of our proposed approach is that we can estimate $\delta$ automatically without predetermining the candidate values of $\delta$. Our simulation studies show that our estimator of $\delta$ has a better finite sample performance than the conventional methods.

The rest of the paper is organized as follows. We introduce in Section 2 the proposed estimation procedure as well as the computational details. Simulation studies are presented in Section 4, and an application to speech production data in Section 3, and finally a summary in Section 5.

2 Method

To ease the notation, the intercept function $\alpha(t)$ can be dropped from model (3) without loss of generality. Let $y_{i}^{*}(t)=y_{i}(t)-\overline{y}(t)$ and $x_{i}^{*}(t)=x_{i}(t)-\overline{x}(t)$ denote the pointwise centered response curves and predictor curves, respectively. A centered model without intercept is

[TABLE]

Upon obtaining an estimate $\widehat{\beta}(s,t)$, the intercept function is then estimated as

[TABLE]

We drop the asterisk from model (4) from now on and focus the discussion on the estimation of the coefficient function $\beta(s,t)$ in the model:

[TABLE]

2.1 Approximation with the Finite Element Method

We propose to approximate the coefficient function $\beta(s,t)$ with the triangular basis, which originates from the finite element method and is widely used in the numerical solution for the boundary-value problems of partial differential equations. Noticing that the support of the coefficient function $\beta(s,t)$ is non-rectangular, approximation with a commonly used bivariate spline generated via tensor product would result in a jagged shape along the boundary $t=s$. Therefore, finite elements that can approximate arbitrary regions serve as a natural alternative.

Let $\phi_{1}(s,t),\ldots,\phi_{K}(s,t)$ denote the $K$ known triangular basis functions. The coefficient function $\beta(s,t)$ is approximated by the expansion,

[TABLE]

Plugging the above approximation into model (4), we obtain

[TABLE]

where $\psi_{ik}(t)=\int_{0}^{t}x_{i}(s)\phi_{k}(s,t)ds$ is known, $\varepsilon_{i}^{{}^{\prime}}(t)$ includes both the residual $\varepsilon_{i}(t)$ and the approximation error.

Divide the interval $[0,T]$ on each axis into $M$ subintervals with equidistant nodes $0=t_{0}<t_{1}<\ldots<t_{M}=T$. Further split each square into two triangles by the diagonal parallel to the line $t=s$. This divides the triangular region over which the coefficient function $\beta(s,t)$ is possibly non-zero into $M^{2}$ congruent triangles (i.e. the triangular elements, or finite elements) with $K=(M+1)(M+2)/2$ nodes. Each node corresponds to a basis function $\phi_{k}(s,t)$ and a coefficient $b_{k}$, both of which are indexed from bottom to top and row-wise from left to right. For example, the left panel in Figure 2 shows the triangulation and indexation of the nodes when $M=5$. There are 25 triangular elements and 21 nodes corresponding to $b_{k}$ and $\phi_{k}(s,t)$ for $k=1,\ldots,21$.

Each basis function $\phi_{k}(s,t)$ has a compact support, defined over the hexagon centered at the $k$-th node. Basis functions of degree 1 are tent shaped, piecewise linear and continuous, with value 1 at node $k$ and value 0 at the boundary of the hexagon. For instance, the right panel in Figure 2 shows a tent-shaped basis function corresponding to the 9th node at $(s,t)=(0.4,0.6)$. Though the $\phi_{k}(s,t)$’s are not orthogonal basis, their compact support property brings in certain computational advantage. We refer to Larson and Bengzon (2013) for a comprehensive introduction to the finite element basis system.

2.2 The Nested Group Bridge Approach

The coefficient function $\beta(s,t)$ in model (5) is formally defined over the triangular region in Figure 1 corresponding to the lag $\delta=T$, while we are trying to identify whether the upper left triangular region is zero or not. Therefore, a proper penalty should be able to shrink the upper left triangular region towards zero specifically, while respecting the nested group structure among the basis coefficients $b_{k}$’s.

We start by defining a sequence of nested triangular regions. Let $\Delta=T/M$ and $D_{j}$ denote the triangle with vertices $(0,T),(0,(j-1)\Delta),((M+1-j)\Delta,T)$, for $j=1,\ldots,M$, and let $D_{M+1}=\{(0,T)\}$ contain a single point. Notice that $D_{1}\supset\ldots\supset D_{M+1}$. In correspondence to these regions, we define a sequence of decreasing groups $A_{1}\supset\ldots\supset A_{M+1}$, where $A_{j}$ consists of the indices of the nodes contained in triangle $D_{j}$ according to the node indexation described in Section 2.1. Taking $M=3$ as an example, there are 10 nodes and the index sets are $A_{1}=\{7,4,8,2,5,9,1,3,6,10\},A_{2}=\{7,4,8,2,5,9\},A_{3}=\{7,4,8\}$, and $A_{4}=\{7\}$. Furthermore, we denote by $\bm{b}_{A_{j}}=\{b_{k}:k\in A_{j}\}$ the vector of basis coefficients whose indices belong to set $A_{j}$, and follow the conventional notations that $||\cdot||_{1}$ and $||\cdot||_{2}$ represent the $L_{1}$ and $L_{2}$ norms, respectively.

We adapt the group bridge approach proposed by Huang et al. (2009) and propose to estimate the vector of basis coefficients $\bm{b}=(b_{1},\ldots,b_{K})$ for $\beta(s,t)$ in model (5) by minimizing the penalized least squares criterion

[TABLE]

with a fixed $\gamma\in(0,1)$, a nonnegative tuning parameter $\lambda$, known weights $c_{j}$ to offset the effect of different dimensions of $A_{j}$, and a known smoothness penalty matrix $\bm{R}$ to be introduced shortly. Following Huang et al. (2009), a simple choice for the weights is $c_{j}\propto|A_{j}|^{1-\gamma}$.

The first term in criterion (6) is the ordinary least squares taking into account the whole time course $[0,T]$. The second term is the so called nested group bridge penalty, which originated from Huang et al. (2009) for simultaneous variable selection at both the group and within-group levels. Consider any two sets $A_{j}$ and $A_{k}$, with $j<k$. Due to the nested nature, the vector of coefficients $\bm{b}_{A_{k}}$ is always a subvector of $\bm{b}_{A_{j}}$, hence the coefficients corresponding to the nodes in region $D_{k}$ appear in more groups than the ones corresponding only to $D_{j}$. This suggests that the nested group bridge penalty shrinks more heavily the coefficients corresponding to regions in a closer proximity to $D_{M+1}$, as desired. The last term imposes a discrete roughness penalty on the basis coefficients $\bm{b}$, extending the idea of Eilers and Marx (1996) and Harezlak et al. (2007). Concerning the smoothness of the coefficient function $\beta(s,t)$, there are three interpretable directions in which we desire to control the changes in the adjacent coefficients, i.e., horizontal, vertical and diagonal directions, as shown in the left panel of Figure 2. Along each direction, the rows of a penalty matrix correspond to the differences in the basis coefficients of adjacent nodes. The discrete smoothness penalty matrix is defined as

[TABLE]

which combines the three directions with nonnegative weights $\omega_{H},\omega_{V}$, and $\omega_{P}$. An example for $M=3$ is as follows, where the columns of $\bm{D}_{H}$, $\bm{D}_{V}$ and $\bm{D}_{P}$ from left to right correspond to the basis coefficients $b_{1},\ldots,b_{10}$,

\bm{D}_{H}=\left[{\begin{array}[]{*{20}{cccccccccc}}0&-1&1&0&0&0&0&0&0&0\\ 0&0&0&-1&1&0&0&0&0&0\\ 0&0&0&0&-1&1&0&0&0&0\\ 0&0&0&0&0&0&-1&1&0&0\\ 0&0&0&0&0&0&0&-1&1&0\\ 0&0&0&0&0&0&0&0&-1&1\\ \end{array}}\right],

\bm{D}_{V}=\left[{\begin{array}[]{*{20}{cccccccccc}}-1&1&0&0&0&0&0&0&0&0\\ 0&-1&0&1&0&0&0&0&0&0\\ 0&0&-1&0&1&0&0&0&0&0\\ 0&0&0&-1&0&0&1&0&0&0\\ 0&0&0&0&-1&0&0&1&0&0\\ 0&0&0&0&0&-1&0&0&1&0\\ \end{array}}\right],

\bm{D}_{P}=\left[{\begin{array}[]{*{20}{cccccccccc}}-1&0&1&0&0&0&0&0&0&0\\ 0&-1&0&0&1&0&0&0&0&0\\ 0&0&-1&0&0&1&0&0&0&0\\ 0&0&0&-1&0&0&0&1&0&0\\ 0&0&0&0&-1&0&0&0&1&0\\ 0&0&0&0&0&-1&0&0&0&1\\ \end{array}}\right].

2.3 Computation

The key difficulty in optimizing the objective function (6) is due to its non-convexity. We can work on an equivalent constrained optimization problem, which is convex and easier to solve (Huang et al. 2009). The objective function (6) involves an integration of the time dependent ordinary least squares criterion over time. Dividing the time course $[0,T]$ into $Q$ equally-spaced subintervals at time points $t_{q}$, $q=0,1,\ldots,Q$, the integrated least squares criterion can be approximated by the finite sum as

[TABLE]

With a sufficiently large $Q$, the above approximation with Riemann sums over a regular partition of the integral domain provides a satisfactory level of precision with acceptable computational load. Define the $Q$-by-$1$ vector $\bm{y}_{i}=(y_{i}(t_{1}),\ldots,y_{i}(t_{Q})^{T}$, and the $Q$-by-$K$ matrix $\bm{\Psi}_{i}$ whose $(q,k)$-th element is $\psi_{ik}(t_{q})$. By stacking the vectors $\bm{y}_{i}$ into a long vector $\bm{y}$ of length $NQ$ and the matrices $\bm{\Psi}_{i}$ into a tall matrix $\bm{\Psi}$ of size $NQ$ by $K$, the approximated integrated least squares criterion is rewritten in the matrix expression as

[TABLE]

The objective function (6) is then rewritten as the familiar penalized least squares

[TABLE]

where $\lambda$ differs from earlier notation by a constant factor $Q/T$.

Define a vector $\bm{\theta}=(\theta_{1},\ldots,\theta_{M+1})$. Then $\widehat{\bm{b}}_{n}$ minimizes criterion (7) if and only if $(\widehat{\bm{b}}_{n},\widehat{\bm{\theta}})$ minimizes

[TABLE]

[TABLE]

for $0<\gamma<1$ and $\tau=[\lambda\gamma^{\gamma}(1-\gamma)^{1-\gamma}]^{1/(1-\gamma)}$ (Huang et al. 2009). Notice that $||\bm{b}_{A_{j}}||_{1}=\sum_{k\in A_{j}}|b_{k}|$ by definition of the $L_{1}$ norm, and that if $b_{k}$ appears in $A_{j}$ then it also appears in all $A_{s}$ with $s<j$. Thus, exchange the order of summation in the second term of (8), rearrange and collect all the individual terms involving $b_{k}$, the objective function becomes

[TABLE]

[TABLE]

where $g_{k}=\sum_{j=1}^{\ell(k)}\theta_{j}^{1-1/\gamma}c_{j}^{1/\gamma}$ with $\ell(k)=\max\{j:b_{k}\in A_{j}\}$. Intuitively, $\ell(k)$ is the number of times that coefficient $b_{k}$ appears in the nested group bridge penalty.

Define $\bm{G}$ as a $K\times K$ diagonal matrix with the $k$th diagonal element $(Ng_{k})^{-1}$, $\bm{b}^{*}=\bm{G}^{-1}\bm{b}$, $\bm{\Psi}^{*}=[\bm{\Psi}^{T},\sqrt{\omega_{H}N}\bm{D}_{H}^{T},\sqrt{\omega_{V}N}\bm{D}_{V}^{T},\sqrt{\omega_{P}N}\bm{D}_{P}^{T}]^{T}\bm{G}$, and $\bm{y}^{*}=(\bm{y}^{T},\bm{0}^{T})^{T}$ of proper length. Finally, the objective function (8) is expressed as

[TABLE]

[TABLE]

where $b_{k}^{*}$ is the $k$th element of vector $\bm{b}^{*}$. The following iterative algorithm is used to compute $\widehat{\bm{b}}_{n}$.

• •

Step 1: Obtain an initial estimate $\bm{b}^{(0)}$ with ordinary least squares.

• •

Step 2: At each iteration $s$, update $\bm{\theta}$ based on ${\bm{b}}^{(s-1)}$ as

[TABLE]

and correspondingly

[TABLE]

[TABLE]

[TABLE]

• •

Step 3: At each iteration $s$, update $\bm{b}$ based on $\bm{\theta}^{(s)}$ by recognizing this is a LASSO problem

[TABLE]

• •

Repeat Step 2 and Step 3 until convergence.

Once obtaining $\widehat{\bm{b}}_{n}$, the estimators for $\beta(s,t)$ and $\delta$ are defined as

[TABLE]

After obtaining $\widehat{\delta}_{n}$, we then refine the estimation for $\beta(s,t)$ by minimizing criterion (7) with $\lambda=0$, i.e. excluding the nested grouped bridge penalty.

When implementing the algorithm, there are two types of tuning parameters, i.e. the shrinkage parameter $\lambda$ (or equivalently $\tau$) and the smoothness parameters as one type, and the number of grids on the time interval $[0,T]$ as the other type. Since the precision of $\widehat{\delta}_{n}$ obviously depends on $M$, we desire a relatively large $M$ in order to achieve a reasonably good estimate of the historical lag $\delta$, and at the same time to capture enough local features of the coefficient function $\beta(s,t)$. This is also consistent with the common strategy when applying penalized least squares approach, where a relatively large number of nodes is preferred and potential overfitting caused by such a choice would be offset via the sparsity and roughness penalty. The shrinkage parameter and smoothness parameters can be chosen via, for example, the Bayesian information criterion (BIC). The effective degrees of freedom for given $\lambda$ and $\bm{\kappa}=(\omega_{H},\omega_{V},\omega_{P})$ can be approximated by

[TABLE]

where $\bm{\Psi}_{s}$ consists of the columns of $\bm{\Psi}$ corresponding to the nonzero coefficients in $\widehat{\bm{b}}_{n}$, and $\bm{R}_{s}$ is obtained in the same way by properly selecting the columns of $\bm{D}_{H}$, $\bm{D}_{V}$ and $\bm{D}_{P}$. The BIC is then approximated by

[TABLE]

and the tuning parameters are selected as the minimizer of $BIC(\lambda,\bm{\kappa})$.

3 Analysis of Real Data

In this analysis, we apply the proposed nested group bridge approach to data from a speech production experiment. It is known that there are over 100 muscles that must be controlled centrally during speech production, such as the muscles of the thoracic and abdominal walls, the neck and face, and the oral cavity, etc. The timing of activation of different muscle groups is a central issue in anatomical and physiological research of speech. Due to the fact that muscle contractions generate electrochemical changes, the noninvasive electromyography (EMG) method is used to collect data on muscle activation. With electrodes attached to the skin over the muscle, potentials as a result of muscle stimulation can be picked up. In this experiment, a subject said the syllable “bob” 32 times. By Newton’s second law, the accelerations of the center of the lower lip reflect the force applied to tissue by muscle contraction, which was recorded as the response curves. Ramsay and Silverman (2002) provided a comprehensive introduction of the background knowledge about this physiology study.

A random subset of 10 smoothed observations for EMG recordings (top) and corresponding lower lip accelerations (bottom) is presented in Figure 3, where the time range is $[0,0.64]$. Considering that the current lip acceleration is very likely a result of recent muscle movement, this motivates us to fit the historical functional linear model (2) with an unknown lag. The EMG recordings are the functional covariates and the lip accelerations are the functional response. The goal of this analysis is to explore the association between EMG recording and lip acceleration, and identify a historical lag if there is any. To apply the nested group bridge approach, we divides the time range [0,0.64] into 20 subintervals with equidistant nodes, i.e. $M=20$. This leads to $K=231$ nodes and triangular basis functions. By minimizing criterion (6) as described in Section 2, we obtain the estimates for both the regression coefficient and the historical lag.

The estimated historical lag is $\widehat{\delta}=0.352$ seconds, indicating there is a historical effect of muscle activation on the lip movement. A 95% confidence interval for the historical lag $\delta$ is $(0.348,0.406)$, constructed via the bootstrap method by re-sampling the residuals and then re-estimating the model. Figure 4 shows the corresponding estimate of the bivariate coefficient function ${\beta}(s,t)$. Dividing $\widehat{\beta}(s,t)$ pointwise by its standard error obtained via the bootstrap method, the plot of the test statistics looks very similar to Figure 4, with slight difference in the magnitude of the values. We observe that most of the coefficients are significant, using a two-sided $z$-test at 5% significance level. The estimate of $\beta(s,t)$ has extraordinarily large values along the diagonal direction of two regions, $(s,t)\in(0,0.06)\times(0,0.06)$ and $(s,t)\in(0.33,0.44)\times(0.33,0.44)$. This is as expected, corresponding to the two productions of the /b/ syllable when the lip is closed and the muscle activation is most influential. The width of the diagonal band with large values varies approximately from 50 to 60 milliseconds, which corresponds to the delay for a neural signal to be transduced into muscle contraction. And peaks at larger lags suggest possible covariation of EMG and lip movement that requires solid physiological knowledge for interpretation.

4 Simulation Studies

Simulation studies are conducted to evaluate the proposed estimation procedure with a nested group bridge penalty (NGB) in comparison with the penalized approach by Harezlak et al. (2007) based on linear model approximation (PLMA). Set $T=1$. Let $S$ denote the trapezoidal support with vertices $(0,0)$, $(1,1)$, $(1-\delta,1)$, $(0,\delta)$ as shown in Figure 1, and $I(\cdot)$ the indicator function. We have considered the following three scenarios, and the true coefficient function $\beta(s,t)$ is shown in Figure 5.

Scenario 1. $\beta(s,t)=10I\{(s,t)\in S_{\varepsilon}\}+10(\frac{\delta}{\varepsilon}+\frac{s}{\varepsilon}-\frac{t}{\varepsilon})I\{(s,t)\in S\backslash S_{\varepsilon}\}$, where $S_{\varepsilon}$ is the trapezoid with vertices $(0,0)$, $(1,1)$, $(1-\delta+\varepsilon,1)$, $(0,\delta-\varepsilon)$. The coefficient function $\beta(s,t)$ is constant over the region $S_{\varepsilon}$, linear over $S\backslash S_{\varepsilon}$ and vanishes at the line $t=s+\delta$.

Scenario 2. $\beta(s,t)=10(1+\frac{s}{\delta}-\frac{t}{\delta})\times I\{(s,t)\in S\}$. The coefficient function $\beta(s,t)$ is linear over its support $S$, with maximum along the line $t=s$ and vanishes at the line $t=s+\delta$.

Scenario 3. Some “holes” are created randomly on the coefficient function $\beta(s,t)$ of Scenario 2, such that $\beta(s,t)$ can vanish inside $S_{\varepsilon}$.

In Scenario 1, we take a small value of $\varepsilon=0.05$, which leads to a sharp drop in the coefficient function $\beta(s,t)$ towards zero when going outside from the region $S_{\varepsilon}$ toward the line $t=s+\delta$. We consider this scenario as the easiest situation to determine the historical lag $\delta$ due to the sharp change, and expect both approaches to perform well and similarly. Scenario 2 is slightly more difficult, since the boundary $t=s+\delta$ becomes more blurred than that in Scenario 1 when data are contaminated with errors. Scenario 3 adds further difficulty with irregularly shaped coefficient function $\beta(s,t)$ and admits the more general assumption that the dependence of response on the historical predictor can vary over time. For all scenarios, the observed covariate curves in the real data example presented in Section 3 are rescaled to time interval $[0,1]$ and are used along with the true coefficient function $\beta(s,t)$ to generate the mean response curves, possibly to mimic real situations. A random subset of 10 covariate curves are shown in the top panel of Figure 3. We refer to Section 3 for details about the example. The errors are additive and generated pointwise from a $N(0,0.5^{2})$ distribution.

For the proposed NGB approach, we set $\gamma=0.5$ and use the ordinary least squares estimate by Malfait and Ramsay (2003) as an initial value for the algorithm. We choose $c_{j}=|A_{j}|^{1-\gamma}/||\bm{b}_{A_{j}}^{(0)}||_{2}^{\gamma}$ similarly to the idea of adaptive LASSO (Zou, 2006). For each scenario, 100 independent replications are simulated. The results are summarized in Table 1. Concerning the estimated root mean squared error of $\widehat{\delta}$, the proposed NGB estimator outperforms the PLMA estimator in all three scenarios, with a much smaller bias and comparable standard deviation. While the PLMA estimator performs well in Scenario 1 with a small bias and a reasonably small standard deviation, the performance deteriorates badly in more difficult scenarios and the estimator exhibits a tendency to underestimate $\delta$. To assess the accuracy of the estimated bivariate coefficient function, we calculate the square root of the mean integrated squared errors (RISEs) over a meshgrid, also presented in Table 1. The NGB and PLMA approaches do not shown discernible difference in terms of the estimated RISE for $\widehat{\beta}(s,t)$.

5 Summary

We have considered in this article the estimation of the historical functional linear model (2) with an unknown forward lag. We propose a nested group bridge approach, tailored for the simultaneous estimation of the historical lag and the regression coefficient function. The nested group bridge penalty is able to shrink not only a group of coefficient corresponding to the designated area towards zero, but also individual coefficients corresponding to the remaining area towards zero. We adopt the triangular basis from the finite element method in order to conform naturally to the non-rectangular domain of the regression coefficient function. The triangular basis system is computationally efficient in the sense that the compact support of the basis functions leads to a sparse design matrix.

Under the historical functional linear model (2), we assume that the historical lag $\delta$ is independent of $t$. It may be of interest to describe more precisely the historical effect of the functional covariate on the functional response via a time-dependent lag $\delta(t)$. Another point of interest is the alternative choice of the smoothness penalty. The discrete penalty originated from P-spline (Eilers and Marx, 1996) involves a large matrix and its inversion. It is of interest to explore other possibilities. We leave the above mentioned points for future work.