Longitudinal functional models with structured penalties

1 Introduction

Technological advancements and increased availability of storage of large datasets have allowed for the collection of functional data as part of time course or longitudinal studies. In the cross-sectional setting, there have been many proposed methods for estimating a regression function in a so-called functional linear model (fLM). This function is a continuous analogue of a vector of (discrete) regression coefficients; it connects the scalar response, y to a functional covariate, w ≡ w ( s ) . Although estimation for fLMs has been well studied, the extension to longitudinally collected functions has not received much attention. Only recently have longitudinal penalized functional regression (LPFR) and longitudinal functional principal component regression (LFPCR) approaches have been proposed to extend the cross-sectional fLM to a longitudinal setting by incorporating subject-specific random intercepts (Goldsmith et al., 2012; Gertheiss et al., 2013). However, a basic assumption of LPFR and LFPCR is that the regression function remains constant over time. Consequently, these methods are not suited for situations in which the association between a functional predictor and scalar response may evolve over time. We propose a technique that extends the analysis of fLMs by relating a scalar outcome to a functional predictor—both observed longitudinally—and estimates a time-dependent regression function.

The method fits into a generalized ridge regression framework by incorporating an scientifically informed context-dependent quadratic penalty term into the estimation process. Our extension of the fLM framework to the longitudinal setting has two major advantages: (a) the regression function is allowed to vary over time and (b) extrinsic, scientifically relevant information about the structure of the regression function can be incorporated directly into the estimation process. We formulate the estimation procedure within a mixed-model framework, making the method computationally efficient and easy to implement.

Ramsay and Dalzell (1991) introduced the term functional data analysis (FDA) in the statistical literature. The cross-sectional fLM with scalar response can be stated as follows (see, e.g., Yao and Müller, 2010):

E ( y | W ) = μ y + ∫ Ω W ( s ) γ ( s ) ds ,

As there is no unique γ that solves this equation, some form of regularization or constraint is required. A common constraint is to impose smoothness on γ . One approach to this is to expand both the regression function γ and predictor functions W in terms of B-splines and then obtain the regularized estimate of γ (Ramsay and Silverman, 1997). Another approach is to express the regression function γ in terms of the empirical orthonormal basis obtained by the eigenfunctions of the covariance of W (i.e., a Karhunen—Loève (K—L) expansion (e.g., Müller, 2005). A third approach, known as penalized functional regression (PFR) (Goldsmith et al., 2011), combines the above two methods. In PFR, a spline basis is used to represent γ and a subset of empirical eigenfunctions is used to represent each W . Another approach is to use a wavelet basis, instead of splines or eigenfunctions, to represent the predictor functions (Morris and Carroll, 2006).

Here we adopt an approach by Randolph et al. (2012) that does not begin by explicitly projecting onto a pre-specified subspace of functions. Instead, prior information about this subspace and its inherent functional structure is incorporated into the estimation process by way of a penalty operator, L . This approach of partially empirical eigenvectors for regression (PEER) exploits the fact that a penalized least squares regression estimate mathematically arises as a series expansion in terms of a set of basis functions determined jointly by the covariance (empirical functional structure) and the penalty (imposed structure). This naturally extends ridge regression (non-structured penalty) and smoothing penalties such as a second-derivative penalty (presuming a smooth regression function). Here we extend the scope of the PEER approach to the longitudinal setting in a manner that allows the estimated regression function γ to vary with time. Our interest lies primarily in the estimation of a regression function β rather than the prediction of an outcome y ; see, for example, Cai and Hall, (2006) for a discussion regarding the substantially different characteristics of these two goals.

The problem we address involves repeated observations from each of N subjects. For each subject, i , at each observation time, t , we collect data on a scalar response variable, y , and a (idealized) predictor function, W . We are interested in longitudinal regression models of the form:

y it = x it ⊤ β + ∫ 0 1 W it ( s ) γ ( t , s ) ds + z it ⊤ b i + ε it .

(1.1)

Our work is motivated by a study in which magnetic resonance (MR) spectra have been collected longitudinally from late-stage human immunodeficiency virus (HIV) patients (Harezlak et al., 2011). We consider global deficit score (GDS) as a scalar response variable, y , and MR spectra as predictor functions, W . The association of GDS with MR spectra and how this association evolves with time are of interest. One MR spectrum is shown in the left panel of Figure 1: The amplitude, W ( s ) , is plotted against the frequency s , normalized to the [ 0 , 1 ] interval ( x -axis). The pattern and amplitudes of the peaks contain information about the concentration of metabolites present in tissue. Each metabolite has a unique spectrum and so one MR spectrum is a mixture of spectra from each individual metabolite (plus background and random noise); see the right panel in Figure 1 that displays spectra from nine metabolites. Consequently, one expects an observed spectrum from tissue to lie near a functional subspace, Q , spanned by the spectra of pure metabolites. The regression function, γ , models the association between y and W and hence, in principle, should also lie near Q . Hence, the subspace Q should be more informative than B-splines or cosine functions that are in some sense ‘external’ to the problem. For this reason, we adopt a methodology that encourages the estimate of γ to be in or near the preferred subspace, Q . The approach is implemented using a decomposition-based penalty, L , which penalizes the estimate of γ lightly if it belongs to Q and strongly if it does not (Randolph et al., 2012).

OPEN IN VIEWER Figure 1

Left panel: An observed MR spectrum from tissue. Right panel: The nine pure metabolite spectra. In each plot, the y -axis represents amplitude and x -axis the frequency of nucleus, s , normalized to [ 0 , 1 ] interval.

The cross-sectional fLM with scalar response has been a focus of various investigations (Faraway, 1997; Ramsay and Silverman, 1997; Fan and Zhang, 2000; Cardot et al., 1999, 2003; Cai and Hall, 2006; Cardot et al., 2007; Reiss and Ogden, 2009), many of which estimate a regression function in two steps. For example, Cardot et al. (2003) first perform principal component regression (PCR), which projects the observed predictor curves onto an empirical basis to obtain an estimate and then use B-splines to smooth the result. Reiss and Ogden, (2009) study several of these methods along with modifications that include versions of PCR using B-splines and second-derivative penalties (cf. Ramsay and Silverman, 1997; Silverman 1996). Extensions of fLM have been made towards generalized linear model with functional predictors (James, 2002; Müller and Stadtmüller, 2005) and quadratic functional regression (Yao and Müller, 2010). We are interested in extending the fLM to a longitudinal setting. To our knowledge, the only published methods addressing the longitudinal functional predictor framework are LPFR (Goldsmith et al., 2012) and LFPCR (Gertheiss et al., 2013). However, both LPFR and LFPCR assume that the regression function in equation (1.1), γ , remains constant over time. In contrast, proposed estimate is based on modelling the regression function γ as a time-dependent linear combination of several time-invariant component functions, { γ d } d = 0 D , each of which is estimated via a penalty operator, L , that is informed by the context of the scientific question.

An important concern for any regularized estimation method, either cross-sectional or longitudinal, is identifiability of the estimate, that is, the lack of uniqueness or, possibly, its instability. In the fLM, this arises from the lack of invertibility of the empirical covariance operator. A finite number of predictor curves means the dimension of the range of this operator is finite and so, as an operator on a infinite-dimensional domain, it has a non-trivial null space. The philosophy behind a penalty-operator approach is that estimation is constrained to the subspace spanned by functions that are jointly determined by W and L . A sufficient condition for uniqueness of this estimate is to assume Null ( W ) ∩ Null ( L ) = { 0 } (Bjorck, 1996; Engl et. al., 2000). We assume this throughout.

Longitudinal functional regression model with time-varying regression function is described in Section 2. In Section 3.1, estimate of regression function is obtained as a generalized ridge estimate (Tikhonov, 1963; Hoerl and Kennard, 1970). Concept of decomposition-based penalty is reviewed in Section 3.2. We have described mixed model approach of estimation as best linear unbiased predictors (BLUP) in Section 4.1 and expressions for the precision of the estimates are derived in Section 4.2. In the supplemental material (available online in journal's website), we present how our longitudinal penalized estimate, along with its bias and precision, can be obtained, under some weak assumptions, in terms of generalized singular vectors.

Simulation results are presented in Section 5. An application to real magnetic resonance spectroscopy (MRS) data using, and a summary of our findings, is presented in Section 6. The methods discussed in this paper have been implemented in the R package refund (Crainiceanu et al., 2012) via the peer() and lpeer() functions.

Let W denote a continuous predictor function defined on the closed interval Ω , without loss of generality Ω = [ 0 , 1 ] . Let W it denote a predictor function from the ith subject ( i = 1 , … , N ) at the ith timepoint ( t = t 1 , … , t n i ). Technically, an observed predictor arises as a discretized sampling from an idealized function and we will assume that each observed predictor is sampled at the same p locations, s 1 , … , s p ∈ [ 0 , 1 ] , with sampling that is appropriately regular and dense enough to capture informative functional structure, as seen, for instance, in the MRS data in Section 6. Let w it : = [ w it ( s 1 ) , … , w it ( s p ) ] ⊤ be the p × 1 vector of values sampled from the realized function W it . Then, the observed data are of the form { y it ; x it ; w it } , where y it is a scalar outcome, x it is a K × 1 column vector of measurements on K scalar predictors and w it is the sampled predictor from the ith subject at time t . Denoting the true regression function at time t by γ ≡ γ ( t , · ) , the longitudinal functional outcome model of interest is

y it = x it ⊤ β + ∫ 0 1 W it ( s ) γ ( t , s ) ds + z it ⊤ b i + ε it ,

(2.1)

The functional structure, indexed by s , and time structure, indexed by t , have somewhat unequal roles in our model, as we assume that the longitudinal observations are more limited in the amount of information relative to the densely sampled s index. For example, γ may vary linearly with time, γ ( t , · ) = γ 0 ( · ) + t γ 1 ( · ) , or quadratically, γ ( t , · ) = γ 0 ( · ) + t γ 1 ( · ) + t 2 γ 2 ( · ) . This is similar in spirit to a linear mixed effects model with linear or quadratic time slope (see, e.g., Fitzmaurice et al., 2004). In general, we assume that γ can be decomposed into several time-invariant component functions γ 0 , … , γ D as

γ ( t , · ) = γ 0 ( · ) + f 1 ( t ) γ 1 ( · ) + … + f D ( t ) γ D ( · ) ,

y it = x it ⊤ β + ∫ 0 1 W it ( s ) { γ 0 ( s ) + f 1 ( t ) γ 1 ( s ) + … + f D ( t ) γ D ( s ) } ds + z it ⊤ b i + ε it .

y = X β + W γ + Z b + ε .

(2.2)

W = W 1 ⋮ W N W i = w it 1 ⊤ f 1 ( t 1 ) w it 1 ⊤ … f D ( t 1 ) w it 1 ⊤ ⋮ ⋮ ⋮ w it n i ⊤ f 1 ( t n i ) w it n i ⊤ … f D ( t n i ) w it n i ⊤

The proposed estimation approach builds on intuition from single-level functional regression that encourages an estimate of γ ≡ γ ( t , · ) to be in or near a ‘preferred’ space via choice of penalty operator (Randolph et al., 2012). To describe the effect of a general penalty operator, L , it is useful to consider the familiar example of a Laplacian penalty, L . The typical heuristic for this arises by viewing γ as a function whose local ‘smoothness’ is informative. In this case, the term | | L γ | | 2 penalizes sharp changes in γ . For our perspective, it is helpful to recall that the dominant eigenvectors of L (those corresponding to the largest eigenvalues) are sharply oscillatory while the least dominant eigenvectors are very smooth. Hence a linear-algebraic view of this is that rather than penalizing sharp changes, smoothness in the estimate is inherited from the eigenproperties of L . More specifically, structure in the estimate arises from the joint eigenproperties of W and L ; this is given by the generalized singular value decomposition (GSVD), as described in the Supplemental section. In general, the least-dominant eigenfunctions of a penalty L will have the largest effect on the estimate. This property can be used to construct a ‘preferred subspace’ by defining a penalty L whose least dominant (or perhaps zero-associated) eigenfunctions are preferred. The estimates of γ 0 , … , γ D are obtained as follows: (a) identify the ‘preferred’ functional subspace Q where γ is expected to belong; (b) define a decomposition-based penalty (see Section 3.2) to encourage the estimate to be close to Q and (c) obtain a penalized estimate of γ 0 , … , γ D . We say ‘encourage’ rather than ‘constrain’ because the estimate is not forced to be in a subspace Q , but rather it is jointly determined by Q and W (as illustrated in Section 5.4). Proposed approach allows different preferred subspace for each of γ d s. In the longitudinal (or t ) dimension, γ is more explicitly and severely constrained by the choice of f 1 , … , f D .

3.1 Generalized ridge estimate

The model described in Section 2 can be written as

y = X β + W γ + ε * ,

(3.1)

where ε * = Z b + ε ∼ N ( 0 , V ) and V = Z Σ b Z ⊤ + σ ε 2 I . For each d = 0 , … , D , let L d be the penalty operator for γ d and let λ d 2 be the associated tuning parameter. The corresponding penalized estimates of β and γ are minimizers of:

| | y − X β − W γ | | V − 1 2 + λ 0 2 | | γ 0 | | L 0 ⊤ L 0 2 + … + λ D 2 | | γ D | | L D ⊤ L D 2 ,

(3.2)

where | | a | | B 2 = a ⊤ B a for some symmetric, positive definite matrix B . A generalized ridge estimate of β and γ minimizing equation (3.2) is obtained as (see Ruppert et al., 2003)

β ̂ γ ̂ = ( C ⊤ V − 1 C + D ) − 1 C ⊤ V − 1 y ,

(3.3)

where, C = [ X W ] , D = blockdiag { 0 , L ⊤ L } and L = blockdiag { λ 0 L 0 , … , λ D L D } . In the supplemental material, we have derived an expression for the generalized ridge estimate γ ̂ , explicitly in terms of the GSVD.

3.2 Decomposition-based penalty

Let γ ̃ d ≡ γ ̃ L d , λ d be the estimate obtained from the penalty operator L d and tuning parameter λ d 2 , for each d = 0 , … , D . For example, L d may denote I p (a ridge penalty) or a second-order derivative penalty (giving an estimate having continuous second derivative). Alternatively, with prior knowledge about potentially relevant structure in a regression function, a targeted decomposition-based penalty can be defined in terms of a subspace defined by such structure (Randolph et al., 2012). To be precise, if it is appropriate to impose scientifically informed constraints on the ‘signal’ being estimated by γ , this prior may be implemented by encouraging the estimate to be in or near a subspace, Q L 2 ( Ω ) . We further assume that the subspace Q is spanned by the functions { q j } j = 1 J .

Returning to our notation that reflects functional predictors observed at p sampling points, we represent Q by the range of a p × J matrix Q whose columns are { q j } j = 1 J , where each q j represents the discretized q j . Consider the orthogonal projection P Q = Q Q + onto the range ( Q ) , where Q + is Moore—Penrose inverse of Q . Then a decomposition penalty is defined as

L Q = ϕ b P Q + ϕ a ( I − P Q )

(3.4)

for scalars ϕ a and ϕ b . To see how L Q works, let γ ̃ d be any estimate of γ d . When γ ̃ d ∈ span ( Q ) , we have L Q γ ̃ d = ϕ b γ ̃ d , but when γ ̃ d ∉ span ( Q ) , we have L Q γ ̃ d = ϕ a γ ̃ d . The condition ϕ a > ϕ b imposes more penalty for γ ̃ d ∉ span ( Q ) , compared to when γ ̃ d ∈ span ( Q ) . The weights ϕ a and ϕ b determine the relative strength of emphasizing Q in the estimation process. Note that L Q is invertible, provided ϕ a and ϕ b are nonzero, and that ϕ a = ϕ b results in ordinary ridge estimate.

We aim to construct an appropriate mixed model that minimizes the expression in equation (3.2). In general, the penalty, L , is not required to be invertible but for simplicity, this will be assumed here. The mixed model approach provides an automatic selection of tuning parameters λ 1 , … , λ D restricted maximum likelihood (REML)’based estimation of the tuning parameters has been shown to perform as well as other criteria and under certain conditions, it is less variable than generalized cross validation (GCV)’based estimation (Reiss and Ogden, 2009).

4.1 Estimation of parameters

Using Henderson's justification (Henderson, 1950), one can show that, for each d = 0 , … , D , model y = X β + W γ + ε * , where ε * ∼ N ( 0 , V ) and γ d ∼ N ( 0 , 1 λ d 2 ( L d ⊤ L d ) − 1 ) minimizes the expression (3.2) to obtain the BLUP. Thus, the generalized ridge estimate of β and γ corresponds to the BLUP from the following model:

y = X β + W * γ * + ε ,

where, W * = [ W Z ] , γ * = [ γ ⊤ b ⊤ ] ⊤ ∼ N [ 0 , Σ γ * ] and ε ∼ N ( 0 , σ ε 2 I ) with Σ γ * = blockdiag { ( L ⊤ L ) − 1 , Σ b } and Σ b = blockdiag { Σ b 1 , … , Σ b N } . This representation allows us to estimate fixed and functional predictors, simply by fitting a linear mixed model (e.g., using nlme package in R or PROC MIXED in SAS).

4.2 Precision of estimates

Our ridge estimate is the BLUP from an equivalent mixed model, hence the variance of the estimate depends on whether the parameters are random or fixed. Randomness of γ is a device used to obtain the ridge estimate while ε and b in our case are truly random. With this in mind, we assume that the variance of the estimates is conditional on γ , but not on b . The BLUP of β , γ and b can be expressed as (see, e.g., Robinson, 1991; Ruppert et al., 2003):

β ̃ = X ⊤ V 1 − 1 X − 1 X ⊤ V 1 − 1 y γ ̃ = ( L ⊤ L ) − 1 W ⊤ V 1 − 1 ( y − X β ̃ )

b ̃ = Σ b Z ⊤ V 1 − 1 ( y − X β ̃ ) ,

where V 1 = V + W ( L ⊤ L ) − 1 W ⊤ . β ̃ is an unbiased estimator of β , but γ ̃ is not unbiased. It is trivial to see that Cov ( y | γ ) = V . Thus, the conditional variances of β ̃ and γ ̃ are:

Cov ( β ̃ | γ ) = X ⊤ V 1 − 1 X − 1 X ⊤ V 1 − 1 V V 1 − 1 X X ⊤ V 1 − 1 X − 1

Cov ( γ ̃ | γ ) = A γ V A γ ⊤ A γ = ( L ⊤ L ) − 1 W ⊤ V 1 − 1 { V 1 − X ( X ⊤ V 1 − 1 X ) ⊤ } V 1 . − 1

(4.1)

To obtain the unconditional variance, one must replace V by V 1 in the above expressions, but this will overestimate the variance of the estimates. The 95% confidence bands can be constructed as Estimate ± 3.201 × ( Standard Error ) , using the expression for the standard error derived here. The constant 3.201 was obtained by Salganik et al. (2004) through a large simulation study to ensure 95% simultaneous coverage.

Expressions for the predicted value of y and its variance are:

y ̃ = X β ̃ + W γ ̃ + Z b ̃ Cov ( y ̃ | γ ) = A y V A y ⊤ ,

where A y = [ { V 1 − W L ⊤ L W − Z Σ b Z ⊤ } − 1 X X ⊤ V − 1 X − 1 X ⊤ V − 1 + W L ⊤ L W ⊤ + Z Σ b Z ⊤ ] V 1 − 1 . Let, T = [ 1 f 1 ( t ) … f d ( t ) ] ⊗ I K . Then the discretized version of regression function at time t is γ t = [ γ ( t , s 1 ) , … , γ ( t , s K ) ] = T γ . Therefore, the estimate of γ t is γ ̃ t = T γ ̃ and the estimate of its variance is T Cov ( γ ̃ | γ ) T ⊤ .

5 Simulation

The proposed method of estimating regression function using decomposition-based penalty will be referred to as longitudinal partially empirical eigenvectors for regression (LongPEER). We pursue several simulations to evaluate its properties. The first simulation study (Section 5.1) compares the performance of LongPEER with the LPFR approach. In the remaining simulation studies, only LongPEER is considered due to a lack of comparable alternatives. The purpose of the second simulation study is to evaluate the influence of sample size and the contribution of prior information about the functional structure (as determined by the tuning parameters ϕ a and ϕ b in equation (3.4)) on the LongPEER estimate. In the third simulation study, we evaluate the coverage probabilities of the confidence bands constructed using the formula presented in Section 4.2. Finally, in Section 5.4, we evaluate the performance of a LongPEER estimate when information on some features is missing. These studies are performed in the context of MRS data and hence the simulated predictor functions are constructed based on the structure of MRS data. All results summarized in this section are based on 100 simulated datasets.

For each subject and visit, predictor functions were simulated using equation (5.3). Predictor functions were flat with bumps of varied widths at a number of pre-specified locations. White noise was added to the predictor functions to account for the instrumental measurement noise. These ‘bumpy’ regression functions were generated with bumps at some (but, not all) of the bump locations of the predictor function. For the simulation in Section 5.1, the regression function is assumed to be independent of time, whereas it varies with time in the simulation of Section 5.2. For both predictor and regression functions, 100 equispaced sampling points in [0,1] are used.

For the decomposition-based penalty (3.4), the matrix L d is defined as follows: (a) select the discretized functions q 1 , … , q J spanning the ‘preferred’ subspace, where each q j is defined to have a single bump corresponding to a region in the simulated predictor functions (see Figure 2) and (b) compute L d = Q Q + , where Q = col [ q 1 , … , q J ] . The columns of Q need not be orthogonal (cf. Figure 9). For ridge and second-order penalties, the matrices L d are defined as I and D 2 , respectively, where D 2 = [ d i , j ] is a square matrix with entries d i , i = d i , i + 2 = 1 , d i , i + 1 = − 2 and d i , j = 0 otherwise.

Our primary interest lies in estimation of regression function; therefore, estimation error is summarized in terms of the mean squared error (MSE) defined as | | γ − γ ̃ | | 2 , where γ ̃ denotes the estimate of γ . We also calculated the sum of squared prediction errors (SSPE) as | | y − y ̃ | | 2 / N , where y ̃ denotes the estimate of the true (noiseless) y to evaluate the prediction performance. The estimates based on the proposed methods, including the LongPEER estimate, were obtained as BLUPs from the mixed model formulation described in Section 4.1. OPEN IN VIEWER

Figure 2

Average estimates of γ for the simulation as described in Section 5.1 with ϕ a = 10 and ϕ b = 1 . Top panel: Prior information used in LongPEER estimation. Bottom panel: True γ and the average estimates from 100 simulations. LongPEER, ridge and second-order difference estimates have been obtained using decomposition-based penalty, ridge penalty and second-order difference penalties, respectively, as described in third paragraph of Section 5. LPFR estimates were obtained using lpfr() with default options in refund package.

5.1 Comparison with LPFR

As mentioned, LPFR estimates a regression function that does not vary with time. Therefore, in order to compare LongPEER with LPFR, in the first set of simulations, we generated outcomes using a time-invariant regression function (i.e., γ ( t , · ) = γ 0 ( · ) , for all t ). The following model was used to generate the outcome data for 100 individuals, each at 4 timepoints ( t = 0 , 1 , 2 , 3 ):

y it = β 0 + ∫ 0 1 W it ( s ) γ 0 ( s ) ds + b i + ε it , i = 1 , … , 100 .

(5.1)

We uniformly sampled ” W ( s ) ” at s = 0.01 , 0.02 , … , 0.99 , 1.00 . γ 0 was defined with bumps centred at H γ 0 = { 0.15 , 0.50 , 0.80 } :

γ 0 ( s ) = ∑ h ∈ H γ 0 a 0 , h exp − c 0 , h * ( s − h 100 ) 2 , for s ∈ [ 0 , 1 ] ,

where a 0 , h and c 0 , h correspond to amplitude and degree of curvature, respectively, at the bumps, as specified in Table 1. We constructed W to consist of wide, moderate and narrow bumps centred at: H wide = { 0.50 } , H moderate = { 0.15 , 0.30 , 0.70 , 0.80 , 0.90 } and H narrow = { 0.05 , 0.95 } . Specifically, we generated the subject-specific functional predictors w it ( s ) (i.e., the observed values of W) of the form

w it ( s ) = ∑ h ∈ H wide ( ξ it w , h + a w , h ) exp − c w , h * s − h 100 2 + ∑ h ∈ H moderate ( ξ it m , h + a m , h ) exp − c m , h * s − h 100 2 + ∑ h ∈ H narrow ( ξ it n , h + a n , h ) exp − c n , h * s − h 100 2 ,

(5.2)

where { a w , h , c w , h } , { a m , h , c m , h } and { a n , h , c n , h } correspond to the amplitude and degree of curvature at the wide, moderate and narrow bumps, respectively. These are defined in Table 1. Note that the amount of curvature in W differs considerably from the amount imposed in γ 0 and γ 1 . ξ it w , h , ξ it m , h and ξ it n , h were drawn independently from uniform(0, 0.1). Also, β 0 = 0.06 , ε it ∼ N [ 0 , ( 0.02 ) 2 ] and b i ∼ N [ 0 , ( 0.05 ) 2 ] .

OPEN IN VIEWER

Table 1

The amplitude and curvature used to generate the predictor function W and regression function γ 0 in simulation studies, as descibed in Sections 5.1–5.4.

	Wide		Moderate		Narrow		H γ 0		H γ 1
	amp.	curv.	amp.	curv.	amp.	curv.	amp.	curv.	amp.	curv.
h	( a w , h )	( c w , h	( a m , h )	( c m , h )	( a n , h )	( c n , h )	( a 0 , h )	( c 0 , h )	( a 1 , h )	( c 1 , h )
5					0.1	2500
15			0.6	1000			0.2	500
30			0.6	1000					0.2	500
50	0.5	250					− 0.15	250
70			0.5	1000					− 0.2	1000
80			0.5	1000			0.15	500
90			0.4	1000
95					0.1	2500

Notes: amp. Amplitude; curv. Curvature Parameters defining features in equation (5.3) are specified in columns Wide, Moderate and Narrow. H_γ0 and H_γ1 represent the set of sampling points, corresponding to bumps in γ0 and γ1, respectively. Note that the degrees of curvature in γ0 and γ1 are different from the ones assumed in W.

Both LPFR (using lpfr() available in the refund package in R (Crainiceanu et al., 2012)) and the LongPEER estimates were obtained. For LPFR estimate, the dimension of both principal components and truncated power series spline basis were set to 30. The prior information used here is represented by bumps of various scalings and locations, as shown in the top panel of Figure 2. This structure is loosely consistent with that present in the predictors. We used ϕ a / ϕ b = 10 , a choice motivated by our findings in Sections 5.2 and 5.4.

Table 2 displays the MSE and prediction error obtained for LongPEER and LPFR estimates. The SSPE was similar for both methods (1.1572 and 1.1527); however, the LongPEER estimate has smaller MSE. Both bias and variance are higher for the LPFR estimate and consequently it has the greater MSE. Figure 2 displays the estimates of the regression function. It should be emphasized that any comparison of these methods is not entirely fair since LongPEER is designed to exploit presumed structural information while LPFR is not. We, note that the ability to exploit such information may be limited and so in this simulation we used imprecise information about the shapes of features; see top panel in Figure 2.

OPEN IN VIEWER

Table 2

Estimation and prediction errors for LPFR and LongPEER estimates based on 100 simulated datasets. The sample size is N = 100 and the number of longitudinal observations is n i = 4 .

	LongPEER	LPFR
MSE ( γ 0 )	0.0257	0.2229
Trace of variance ( γ 0 )	0.0057	0.0653
| | Bias ( γ 0 ) | | 2	0.0200	0.1576
SSPE of Y	1.1572	1.1527

5.2 Simulation with a time-varying regression function

Here and in subsections 5.3 and 5.4, the regression function varies parametrically with time. Due to the lack of existing alternatives for time-varying functional regression for comparison, we instead evaluate the performance of our approach in several ways. First, our primary goal was to assess the effects of sample size, fraction of variance explained by the model and the relative contribution of extrinsic information (as determined by ϕ a and ϕ b in equation (3.4)) on estimate. Without loss of generality, we set ϕ b = 1 and vary ϕ a on an exponential scale. Larger values of ϕ a indicate greater emphasis of prior information on the estimation process. The model considered here was similar to that described in Section 5.1, with the exception that γ ( t , · ) = γ 0 ( · ) + t γ 1 ( · ) . The function γ 0 is defined in equation (5.1) and γ 1 was defined with bumps centred at H γ 1 = { 0.30 , 0.70 }

γ 1 ( s ) = ∑ h ∈ H γ 1 a 1 , h exp − c 1 , h * ( s − h 100 ) 2

where a 1 , h and c 1 , h correspond to amplitude and degree of curvature, respectively, at the bumps as specified in Table 1. As before, it was assumed that β 0 = 0.06 . Realizations of functional predictors were generated, as described in Section 5.1. For each simulation, an appropriate σ ε 2 was chosen to ensure that the squared multiple correlation coefficient R 2 = s y 2 / [ s y 2 + σ ε 2 ] is 0.6 and 0.9 . Here, s y 2 = 1 4 ∑ t = 0 3 1 N − 1 ∑ i = 1 N ( y it − y ̅ . t ) 2 denotes the average sample variance in the set { y it − ε it : i = 1 , … , N ; t = 0 , … , 3 } with y ̅ . t = 1 N ∑ i = 1 N y it .

We have repeated the simulation for four scenarios: (a) N = 100 , R 2 = 0.6 ; (b) N = 100 , R 2 = 0.9 ; (c) N = 200 , R 2 = 0.6 ; and (d) N = 200 , R 2 = 0.9 . The prior information used to obtain LongPEER estimates of γ 0 and γ 1 is plotted in the top panel of Figure 5. Results for AIC, MSE and SSPE are displayed graphically in Figure 3. The standard deviation of MSE is plotted in Figure 4. As the sample size and R 2 increased, both MSE ( γ 0 ) and MSE ( γ 1 ) decreased, providing empirical evidence that the LongPEER estimates were consistent. In general, MSE( γ 0 ) (except for N = 100 , R 2 = 0.60 ) and MSE( γ 1 ) were minimized at ϕ a = 10 and plateaued after that, that is, an increase in ϕ a up to 10 resulted in improvement in estimation of both γ 0 and γ 1 and thereafter estimation performance remained almost unchanged. Finally, note that the value of ϕ a that maximized AIC also minimized MSE ( γ 0 ) and MSE ( γ 1 ) . These results suggest that (a) AIC can be used to guide the choice of ϕ a , and (b) with a sufficiently large ratio of ϕ a / ϕ b , there is minimal impact on the estimation performance. OPEN IN VIEWER

Figure 3

Average AIC, SSPE and MSE for simulations, as described in Section 5.2 over 100 simulations. At ϕ a = 10 , average AIC were maximized and MSE( γ 0 ) and MSE( γ 1 ) were minimized. In general, average AIC increased with the increase in sample size and R 2 whereas SSPE, MSE( γ 0 ) and MSE( γ 1 ) decreased.

OPEN IN VIEWER Figure 4

Standard deviation of MSE( γ 0 ) and MSE( γ 1 ) over 100 simulations, as described in Section 5.2. Both standard deviations decreased with increasing sample size and R 2 and minimized at ϕ a = 10 and plateaued after that.

OPEN IN VIEWER Figure 5

Average estimates of the components of regression functions for simulations, as described in Section 5.2 with ϕ a = 10 and ϕ b = 1 . Top panel: prior information used in LongPEER estimation. Middle and bottom panels: the average estimates of γ 0 and γ 1 ; these improve as N and/or R 2 increase.

The average LongPEER estimates of γ 0 and γ 1 using a decomposition penalty are displayed in Figure 5, with ϕ a = 10 and ϕ b = 1 . For smaller sample sizes and R 2 , the LongPEER estimate may: (a) oversmooth (i.e., negatively bias) the estimated regression function at locations of a true feature and (b) be positively biased in locations corresponding to features in Q but the true γ is zero. However, by increasing the sample size to 200 and/or R 2 to 0.9, we observe that the average LongPEER estimates γ 0 and γ 1 approach the true functions.

5.3 Coverage probability OPEN IN VIEWER

Figure 6

Coverage probabilities of LongPEER estimates in 100 simulations with ϕ a = 10 and ϕ b = 1 , as discussed in Section 5.3. Top panel: Prior information used in LongPEER estimation. Middle and bottom panels: 95% confidence band (shaded region) and coverage proportions (the dotted line) based on N = 100 , and N = 400 subjects, respectively. The left column displays the cross-sectional function γ 0 and the right column displays the longitudinal function γ 1 . The horizontal line in each plot marks the nominal coverage of 95%.

In this section, data were generated according to simulation scenario described in Section 5.2 with R 2 = 0.9 . The prior information used in LongPEER estimation along with the confidence bands and the coverage probabilities of LongPEER estimates are displayed in Figure 6. The 95% confidence bands are constructed as discussed in Section 4.2.

A notable improvement in coverage probabilities for both γ 0 and γ 1 was observed with the increase of N . For N = 400 , the coverage for both γ 0 and γ 1 was close to nominal level of 95%. We can see the dip in coverage probability in the first feature in γ 0 at s = 0.15 due to negative estimation bias (i.e., over-smoothing); however, the coverage probability improves when N increases. We also explored the influence of ϕ a on the confidence band and coverage probability (not shown here). The higher values of ϕ a led to the confidence band shrinkage and this in turn resulted in under-coverage of both γ 0 and γ 1 .

5.4 Estimation in the presence of incomplete information OPEN IN VIEWER

Figure 7

Top panel: true regression functions (solid lines) γ 0 (left) and γ 1 (right) and prior information about features (dashed lines). Note that there is no information about the peak at s = 0.5 in γ 0 . Middle panel: Average ridge penalty estimate from 100 simulations. Bottom panel: Average LongPEER estimate from 100 simulations using prior information is displayed in the top panel and ϕ a = 10 0.75 .

Since the LongPEER estimate uses extrinsic scientific information in the estimation process, it is of interest to evaluate its estimation performance when only partial information is available. For this, we considered a simulation scenario similar to Section 5.3, but LongPEER estimate was obtained disregarding the information about the peak at s = 0.5 . LongPEER estimates with ϕ b = ϕ a = 1 (i.e., a ridge penalty) and ϕ b = 1 , ϕ a = 10 0.75 are displayed in Figure 7. LongPEER estimates of γ 0 ( s ) have an appropriate structure at s = 0.5 , on average. Indeed as with an ordinary ridge penalty, this structure is inherited from the empirical eigenfunctions of W . This highlights the advantage of an estimate obtained from the jointly determined eigenfunctions of W and L (see supplemental material available online in journal's website); the estimate depends on the relative contributions of predictor function W and penalty L (i.e., relative contribution of extrinsic information), controlled by ϕ a / ϕ b . The relative increase in the contribution of extrinsic information resulted in shrinkage towards zero at s = 0.5 in LongPEER estimate.

OPEN IN VIEWER

Figure 8

Prediction performance of model, as described in equation (6.1). Left panel: observed GDS score ( y ) and predicted value ( y ̃ ). Right panel: observed y ̃ and residuals ( y − y ̃ ) .

We applied LongPEER to investigate potential associations of metabolite spectra, obtained from basal ganglia, and the global deficit score (GDS) in a longitudinal study of late-stage HIV patients. Of particular interest is how such an association evolves over time. The study cohort was comprised of chronically HIV-infected patients enrolled in the HIV Neuroimaging Consortium (HIVNC), a longitudinal study of HIV-associated brain injury, at multiple sites in the United States. At the time of enrolment, patients were on a stable anti-retroviral treatment. More detailed study description is available in Harezlak et al. (2011). We treat GDS as our scalar continuous response variable and MR spectrum (sampled at K = 399 distinct frequencies) as functional predictor. GDS is often used as a continuous measure of neurocognitive impairment (e.g., Carey et al., 2004) and a large GDS score indicates a high degree of impairment. The MRS spectra are comprised of pure metabolite spectra, instrument noise and a background profile. We collected a total of n • = 306 observations from N = 114 subjects. The longitudinal observations for each subject were within three years from baseline. The number of observations per subject ranged from 1 to 5 with a median equal to 3. Spectral information of following nine pure metabolites was used as prior information for the LongPEER estimation: Creatine (Cr), Glutamate (Glu), Glucose (Glc), Glycerophosphocholine (GPC), myo-Inositol (Ins), N-Acetylaspartate (NAA), N-Acetylaspartylglutamate (NAAG), scyllo-Inositol (Scyllo) and Taurine (Tau). These spectra are displayed in Figure 1.

Information available on demographic factors includes: age at baseline, gender and race. We relied on AIC to choose (a) scalar covariates in the model, (b) ϕ a (while setting ϕ b = 1 ) for determining the relative contribution of prior information and (c) the time structure of γ ≡ γ ( t , · ) . Based on the AIC (see Table 3), Models 1, 3, 5 and 6 are almost identical and appear to be better than the remaining models. In this illustration, we choose to consider relatively simple model. We fit Model 1 (with ϕ a = 100 , ϕ b = 1 ) as follows:

y it = β 0 + β 1 t + ∫ Ω W it ( s ) γ ( t , s ) ds + b i + ε it ,

(6.1)

The estimates of λ (tuning parameter) associated with γ 0 and γ 1 were 1.152 and 2.242 , respectively, and the estimates of σ ε 2 and σ b 2 were 0.0786 and 0.3332, respectively. The GDS score, fitted values and residual plot are displayed in Figure 8 for the purpose of model checking. The residuals do not show any obvious pattern.

OPEN IN VIEWER Figure 9

LongPEER estimates of γ 0 and γ 1 with ϕ a = 100 , ϕ b = 1 as described in Section 6. Shaded region in both the plots represent 95% pointwise confidence bands. Selected (scaled) pure metabolite spectra are also shown on both plots.

OPEN IN VIEWER

Table 3

Comparison of AIC for selection of scalar covariates, ϕ a ( ϕ b = 1 ) and time structure in γ ≡ γ ( t , · ) , as descibed in Section 6.

	Scalar covariates	Time structure in γ ( t , s )	ϕ a	AIC
Model 1	t	γ 0 ( s ) + t γ 1 ( s )	10	− 395.2335
Model 2	Age, t	γ 0 ( s ) + t γ 1 ( s )	10	− 405.2796
Model 3	Gender, t	γ 0 ( s ) + t γ 1 ( s )	10	− 395.9040
Model 4	Race, t	γ 0 ( s ) + t γ 1 ( s )	10	− 398.5607
Model 5	t , t 2	γ 0 ( s ) + t γ 1 ( s ) + t 2 γ 2 ( s )	10	− 394.5752
Model 6	t	γ 0 ( s ) + t γ 1 ( s )	100	− 395.367 0

Figure 9 displays the estimates of γ 0 and γ 1 with pointwise 95% confidence bands. To aid interpretation, selected pure metabolite spectra are displayed. These figures reveal that estimated γ 0 (the ‘baseline’ part of the regression function) is different from zero at the locations where at least one of the pure metabolites Cr, Glu, NAA, NAAG and Scyllo has a bump. Similarly, each non-zero part of estimated γ 1 (the ‘longitudinal’ part of the regression function) coincides with bump locations of one or more pure metabolite profiles of Cr, Glu, NAA, GPC and Ins.

Pointwise confidence intervals for γ 0 and γ 1 contain the 0 line over large intervals. The estimated γ 0 is significant in the region s ∈ ( 0.4 , 0.5 ) ∪ ( 0.6 , 0.8 ) ) and estimated γ 1 is significant in a region s ∈ ( 0.5 , 0.6 ) . To be precise, peaks in both γ ̂ 0 ( · ) and γ ̂ 1 ( · ) are significant at locations where at least one of the pure metabolite profiles NAA or Glu have bumps. The observation of negative ‘longitudinal’ effect of NAA is worth commenting; it suggests that GDS increases as NAA concentration decreases in basal ganglia, a finding consistent with several studies in which a reduced concentration of NAA is seen to be associated with a decrease in neuronal mass (Christiansen et al., 1993; Lim and Spielman, 1997; Soares and Law, 2009).

Finally, we considered other forms of f ( t ) , such as exp ( t ) − 1 or log ( t + 1 ) . However, relative improvement in AIC due to use of these complex structure was very minimal. For example, AIC associated with γ ( t , · ) = γ 0 ( · ) + [ exp ( t ) − 1 ] γ 1 ( · ) was observed as − 394.56 in comparison to AIC of − 395.23 observed with γ ( t , · ) = γ 0 ( · ) + t γ 1 ( · ) .

We have proposed a novel estimation method for longitudinal functional regression and derived some properties of the estimated regression function. A valuable contribution of this framework is that it allows this estimate to vary with time as it extends the scope of penalized regression to the realm of longitudinal data. The approach may be viewed as an extension of longitudinal mixed effects models, replacing scalar predictors by functional predictors. Advantages of this framework include: estimating a time-dependent regression function, ability to incorporate structural information into the estimation process and easy implementation through linear mixed model equivalence.

The first simulation study of Section 5.1 illustrates the potential advantage of LongPEER estimate in exploiting an informed structured penalty, as compared to the more generic smoothness or spline-based constraints. The simulation in Section 5.3 suggests that coverage probabilities of the confidence bands for the true regression function are close to the nominal level. However, for small sample sizes, the naive confidence bands do not seem to be sufficient and an alternative solution that takes into account the estimation bias is needed. In the case when only partial information is available, the proposed method can be still useful, if we limit the relative contribution of the ‘informed’ space and/or increase the sample size (see Section 5.4). In the absence of prior information, one may impose more vaguely defined constraints—such as identity penalties, smoothing penalties or re-weighted projections onto empirical subspaces—to estimate the regression function.

Estimation in generalized ridge regression can be expressed in many forms. Clearly, one natural way to view this is via a Bayesian equivalence formulation (see, e.g., Robinson, 1991) with the informative priors quantifying the available scientific knowledge. In our formulation, the linear mixed model equivalence provides an easy computational implementation as well as an automatic choice of the tuning parameters using REML criterion. As detailed in the Supplemental material, GSVD provides intuition and insight about the role of the penalty operator via an explicit algebraic formulation of the estimate. Indeed, estimate of γ is obtained from a simultaneous diagonalization of W and L . The GSVD also provides explicit derivations of the bias and variance for the estimate.

A possible extension of this work is to incorporate multiple functional predictors. Consider we have two predictor functions W ( 1 ) and W ( 2 ) and the associated regression functions are γ ( 1 ) ≡ γ ( 1 ) ( t , · ) and γ ( 2 ) ≡ γ ( 2 ) ( t , · ) . Now, express γ ( 1 ) = γ 0 ( 1 ) + f 1 ( 1 ) ( t ) γ 1 ( 1 ) + … + f d ( 1 ) ( 1 ) ( t ) γ d ( 1 ) ( 1 ) and γ ( 2 ) = γ 0 ( 2 ) + f 1 ( 2 ) ( t ) + γ 1 ( 2 ) + … + f d ( 2 ) ( 2 ) ( t ) γ d ( 2 ) ( 2 ) . Then the component functions can be estimated as BLUP from the mixed model:   y = X β + W ( 1 ) γ ( 1 ) + W ( 2 ) γ ( 2 ) + Z b + ε , where W ( 1 ) , and W ( 2 ) represent design matrices for the two functional predictors. Another possible extension is when response variable is binary or count. Indeed, an important problem that arises is to understand the neurocognitive impairment status of HIV patients, based on MRS collected over time. Estimation in these general settings appears to be possible with the proposed framework.