Estimating Functional Linear Mixed-Effects Regression Models

0.49 MB·

by  Baisen Liu

Estimating Functional Linear Mixed-Effects Regression
Models
Baisen Liua, Jiguo Caob,∗
aSchool of Statistics, Dongbei University of Finance and Economics, Dalian 116025, China
bDepartment of Statistics and Actuarial Science, Simon Fraser University, Burnaby, BC
6 V5A1S6, Canada
1
0
2

n
a Abstract
J

6 The functional linear model is a popular tool to investigate the relationship

] betweenascalar/functionalresponsevariableandascalar/functionalcovariate.
E
M We generalize this model to a functional linear mixed-effects model when re-
. peated measurements are available on multiple subjects. Each subject has an
t
a
individual intercept and slope function, while shares common population in-
t
s
[ tercept and slope function. This model is flexible in the sense of allowing the

1 slope random effects to change with the time. We propose a penalized spline
v
9 smoothing method to estimate the population and random slope functions. A
3
0 REML-based EM algorithm is developed to estimate the variance parameters
1
for the random effects and the data noise. Simulation studies show that our es-
0
.
1 timation method provides an accurate estimate for the functional linear mixed-
0
effectsmodelwiththefinitesamples. Thefunctionallinearmixed-effectsmodel
6
1
is demonstrated by investigating the effect of the 24-hour nitrogen dioxide on
:
v
i the daily maximum ozone concentrations and also studying the effect of the
X
daily temperature on the annual precipitation.
r
a
Keywords: EM algorithm, Functional Linear Regression, Penalized Splines,
Random Effects Model
∗Correspondingemail: jiguo [email protected]
Preprint submitted to Computational Statistics & Data Analysis January 7, 2016
1. Introduction
When a random variable is measured or observed at multiple time points
or spatial locations, the data can be viewed as a function of time or spatial
locations. This type of data is generally called as functional data (Ramsay
and Silverman, 2005). In the current big data era, functional data analysis
(FDA) has become very popular in statistical methodology and applied data
analysis. Functionallinearmodels(FLMs)isoneofthemostpopularmodelsin
FDA.Itmodelstherelationshipbetweenfunctionalvariablesand/orpredictsthe
scalar response from the functional input. FLMs have been studied extensively
since Ramsay and Dalzell (1991) introduced them. With the developments
of modern technology, FLMs have been popularly applied to model functional
data in many fields such as economics, medicine, environment, climate [see for
instance, Ramsay and Silverman (2002), Ramsay and Silverman (2005), and
Ferraty and Vieu (2006), for several case studies].
There is extensive literature studying estimations and properties of FLMs.
For example, Chiou et al. (2003) applied a quasi-likelihood approach to study
a FLM with a functional response and a finite-dimensional vector of scalar
predictors. Yao et al. (2005) studied FLMs for sparse longitudinal data and
suggestedanonparametricestimationmethodbasedonthefunctionalprincipal
components analysis (FPCA). Their proposed functional regression approach is
flexible to allow for different measurement time points of functional predictors
andthefunctionalresponse. CaiandHall(2006)discussedthepredictionprob-
lem in FLMs based on the FPCA technique. Crambes et al. (2009) proposed
a smoothing spline estimator for the functional slope parameter, and extended
it to covariates with measurement-errors. Yuan and Cai (2010) suggested a
smoothness regularization method for estimating FLMs based on the reproduc-
ing kernel Hilbert space (RKHS) approach. They provided a unified treatment
for both the prediction and estimation problems by developing a tool on si-
multaneous diagonalization of two positive-definite kernels. Wu et al. (2010)
proposedavarying-coefficientFLMwhichallowsfortheslopefunctiondepend-
2
ing on some additional scalar covariates. A systematic review on FLMs can be
found in Morris (2015).
OnepopularFLMistolinkascalarresponsevariableY ,j =1,...,m,with
j
a functional predictor X (t) through the following model
j
(cid:90)
Y =α+ β(t)X (t)dt+(cid:15) , (1)
j j j
S
where α is the intercept, β(t) is a smooth slope function, (cid:15) ’s are independent
j
and identically distributed (i.i.d.) random variables with mean 0 and variance
σ2, and S is often assumed to be a compact subset of an Euclidean space such
(cid:15)
as [0,1]. The slope function, β(t), represents the accumulative effect of the
functional covariate X (t) on the scalar response Y .
j j
For purposes of illustration, we take the air pollution data as an example.
This data is from the R package NMMAPSdata (Peng and Welty, 2004). The
data have hourly measurements of ozone and nitrogen dioxide NO concentra-
2
tions for some U.S. cities. Our aim is to study the relationship of the daily
maximum ozone concentration and the functional predictors nitrogen dioxide
NO (t) measured during 24 hours (from 0 am to 11 pm) of that day. For the
2
i-thcity, thescalarresponseY isthemaximumozoneconcentrationduring24
ij
hours (from 0 am to 11 pm) in the j-th day, and the functional covariate X (t)
ij
isthehourlyNO (t)concentrationmeasuredduring24hours. Inapreliminary
2
analysis,weperformedafunctionallinearregressionmodel(Cardotetal.,2007)
on each individual city and found that there was a dramatic variation of the
estimated βˆ(t). This indicates that each city has different effects of the hourly
NO (t) concentration on the daily maximum ozone. Therefore, it may not be
2
appropriate to pool all the data of U.S. cities together and provide only one av-
erageeffectofhourlyNO (t)onthedailymaximumozone. Ontheotherhand,
2
we may not use all of the data information available if we fit the functional
linear model for each individual city separately.
To address this dilemma, we generalize the FLM (1) to incorporate random
effects into the slope function, and call it the functional linear mixed-effects
model (FLMM). Assume that we repeatedly observed a distinct functional pre-
3
dictorandscalaroutcomeforeachsubjectoverseveralvisits. Thentheobserved
data has the structure {Y ,X (t)} for 1≤i≤n and 1≤j ≤m , where Y is
ij ij i ij
the j-th repeated measurement of the scalar response for the i-th subject, and
X (t) is the corresponding functional predictor. The functional linear mixed-
ij
effects model can be expressed as following:
(cid:90)
Y =α +a + [β(t)+b (t)]X (t)dt+(cid:15) , i=1,2,...,n,j =1,2,...,m , (2)
ij 0 i i ij ij i
S
where α is the population intercept, a is the intercept random effect, β(t)
0 i
represents the population effect of X (t) on Y , b (t) stands for the random
ij ij i
effect of X (t) on Y for the i-th subject, and (cid:15) is the i.i.d. random variable
ij ij ij
with mean 0 and variance σ2. In this article, we assume that a ∼ N(0,σ2),
(cid:15) i a
(cid:15) ∼N(0,σ2), and b (t) follows a Gaussian stochastic process with mean 0 and
ij (cid:15) i
covariance function γ(s,t), that is, b (t)∼GP(0,γ(s,t)). We also assume that
i
a , (cid:15) , b (t), and X (t) are mutually independent. The above functional linear
i ij i ij
mixed-effects model is very attractive, because it can estimate the population
effectandrandomeffectofthefunctionalpredictorX(t)(e.g. thehourlyNO (t)
2
in the air pollution study) on the scalar response Y(e.g. daily maximum ozone
concentrations) as well as the population intercept and intercept random effect
simultaneously. The application of the proposed functional linear mixed-effects
modelontheairpollutionproblemisnotunique,andmanysimilarapplications
can be found in environmental or biological problems.
The proposed functional linear mixed-effect model (2) is different from the
following functional mixed model (Goldsmith et al., 2011, 2012):
(cid:90)
Y =Z b + β(t)X (t)dt+(cid:15) , (3)
ij i i ij ij
S
where b ∼ N(0,σ2I) accounts for correlations in the repeated outcomes for
i b
the i-th subject. The highlight of the distinction between (2) and (3) is: the
subject-specificrandomeffectb of(3)remainsthesameacrossvisits,whilethe
i
random effect b (t) of (2) allows for varying with time. The including of the
i
randomeffectb (t)in(2)cancharacterizethedifferenttrendeffectoffunctional
i
predictor on scalar outcomes for different subjects.
4
Many nonparametric smoothers used for the FLMs can be applied to fit
the model (2). In this article, we use the idea of penalized splines smoothers of
RamsayandSilverman(2005)toestimateβ(t)andb (t)in(2). Then,themodel
i
(2)istransformedbyalinearmixed-effectsmodel(LMM).ThenaREML-based
EM algorithm is proposed to fit the LMMs, and its efficiency is illustrated by
examples.
The remainder of this article is organized as follows. Section 2 introduces a
smoothing spline method to estimate the above functional linear mixed-effects
model. Section 3 implemented some simulations to evaluate the finite sample
performanceofthesmoothingsplinemethod. Thenthefunctionallinearmixed-
effectsmodelisdemonstratedbytworealapplicationsinSection4. Conclusions
are given in Section 5.
2. Method
Without giving any parametric assumption on the slope functions, β(t) and
b (t), we estimate them as linear combinations of splines basis functions
i
J K
(cid:88) (cid:88)
β(t)= c φ (t)=φ(cid:48)(t)c, b (t)= b ψ (t)=ψ(cid:48)(t)b .
j j i ik k i
j=1 k=1
where φ(t) = (φ (t),...,φ (t))(cid:48) and ψ(t) = (ψ (t),...,ψ (t))(cid:48) are two vectors
1 J 1 K
of basis functions with dimensions J and K, respectively, and c = (c ,...,c )(cid:48)
1 J
and b = (b ,...,b )(cid:48) are the corresponding vectors of basis coefficients to
i i1 iK
estimate. LetD=Cov(b )=E(b b(cid:48))denotethevariance-covariancematrixof
i i i
random-effects, then b ∼ N(0,D), and the covariance function γ(s,t) for the
i
random effect b (t) can be expressed as γ(s,t)=ψ(cid:48)(s)Dψ(t). We first consider
i
the scenario of the functional predictors X (t) observed without measurement
ij
errors. WhenX (t)isobservedwithmeasurementerrors,manynonparametric
ij
smoothing approaches can be applied to reconstruct the underlying functional
predictors X (t), such as the functional principal component analysis method
ij
(Yao et al., 2005), which is introduced in Subsection 2.5.
5
2.1. Estimating Fixed and Random Effects
Let θ = (α ,c(cid:48))(cid:48) and ξ = (ξ(cid:48),...,ξ(cid:48) )(cid:48) with ξ = (a ,b(cid:48))(cid:48),i = 1,...,n. The
0 1 n i i i
fixedeffects{α ,β(t)},andrandomeffects{a ,b (t)}areestimatedbyminimiz-
0 i i
ing
(cid:88)n (cid:88)mi 1 (cid:18) (cid:90) (cid:19)2 1(cid:88)n
H(θ,ξ) = Y −α −a − [β(t)+b (t)]X (t)dt + b(cid:48)D−1b
2σ2 ij 0 i i ij 2 i i
i=1j=1 (cid:15) S i=1
+(cid:20)λβ (cid:90) (cid:26)d2β(t)(cid:27)2dt+ λb (cid:88)n (cid:90) (cid:26)d2bi(t)(cid:27)2dt(cid:21)+ 1 (cid:88)n a2 . (4)
2 dt2 2 dt2 2σ2 i
S i=1 S a i=1
Define three vectors Y = (Y ,...,Y )(cid:48), W = [w ,··· ,w ](cid:48) and Z =
i i1 imi i i1 imi i
[z ,··· ,z ](cid:48)wherew =(cid:0)1,(cid:82) φ(cid:48)(t)X (t)dt(cid:1)(cid:48)andz =(cid:0)1,(cid:82) ψ(cid:48)(t)X (t)dt(cid:1)(cid:48).
i1 imi ij S ij ij S ij
Define a J ×J matrix G = (cid:82) (d2φ(t)/dt2)(d2φ(t)/dt2)(cid:48)dt and a K ×K ma-
S
trix G = (cid:82) (d2ψ(t)/dt2)(d2ψ(t)/dt2)(cid:48)dt. Then H(θ,ξ) can be expressed in a
b S
matrix form
n n
(cid:88) 1 1(cid:88)
H(θ,ξ) = (cid:107)Y −W θ−Z ξ (cid:107)2+ b(cid:48)D−1b
2σ2 i i i i 2 i i
i=1 (cid:15) i=1
n n
+(λβc(cid:48)Gc+ λb (cid:88)b(cid:48)G b )+ 1 (cid:88)a2 .
2 2 i b i 2σ2 i
i=1 a i=1
Then the estimates for θ and ξ are obtained by minimizing H(θ,ξ):
(cid:16) (cid:17)−1
θˆ = W(cid:48)V(cid:101)−1W+λG(cid:101) W(cid:48)V(cid:101)−1Y,
(5)
ξˆ = (In(cid:78)D(cid:101)ξ)Z(cid:48)V(cid:101)−1(Y−Wθˆ).
where Y = (Y1(cid:48),...,Yn(cid:48))(cid:48), W = (W1(cid:48),...,Wn(cid:48))(cid:48), V(cid:101) = diag(V˜1,...,V˜n) with
V˜i = ZiD(cid:101)ξZ(cid:48)i +σ(cid:15)2Imi,i = 1,...,n, G(cid:101) = diag(0,G), D(cid:101)ξ = (D−ξ1 +λbGξ)−1,
D = diag(σ2,D), G = diag(0,D ), Z = diag(Z(cid:48),...,Z(cid:48) ), and (cid:78) denotes the
ξ a ξ b 1 n
kronecker product.
Once obtaining the estimates θˆ = (αˆ0,(cid:98)c(cid:48))(cid:48) and ξˆi = (aˆi,b(cid:98)(cid:48)i)(cid:48), the estimates
of β(t) and b (t),i=1,...,n, can be given by
i
βˆ(t)=φ(cid:48)(t)(cid:98)c, ˆbi(t)=ψ(cid:48)(t)b(cid:98)i, i=1,...,n. (6)
6
2.2. The REML-based EM algorithm
To estimate the fixed-effects, θ, the random-effects, ξ, and the variance
parameters, σ2, σ2 and D, we recommend an EM algorithm procedure called
a (cid:15)
the REML-based EM-algorithm. It was proposed by Wu and Zhang (2006)
for estimating nonparametric mixed-effects regression models with longitudinal
data. The REML-based EM-algorithm has three steps, which are outlined as
follow.
Initializing. Initializing the starting values for σ2, σ2 and D, denoted by
a (cid:15)
σ2(0),σ2(0)andD(0) ,respectively. Forexample,wecanchooseσ2(0) =σ2(0) =1
a (cid:15) a (cid:15)
and D(0) as an identity matrix.
Step 1. Set r =r+1. Compute
D(cid:101)(ξr−1) = [{Dξ(r−1)}−1+λbGξ]−1,
V˜i(r−1) = ZiD(cid:101)ξ(r−1)Z(cid:48)i+σ(cid:15)2(r−1)Imi,i=1,2,...,n.
Denote V˜(r−1) =diag(V˜(r−1),...,V˜(r−1)). Then estimate θˆ(r) and ξˆ(r) by
1 n i
θˆ(r) = [W(cid:48){V˜(r−1)}−1W+λ G˜]−1W(cid:48){V˜(r−1)}−1Y,
β
ξˆ(ir) = D(cid:101)(ξr−1)Z(cid:48)i{V˜i(r−1)}−1(Yi−Wiθˆ(r)),i=1,2,...,n.
Step 2. Compute the residualsˆ(cid:15)(r) =Y −W θˆ(r)−Z ξˆ(r) and the matrix
i i i i i
H(r−1) ={V˜(r−1)}−1−{V˜(r−1)}−1W [W(cid:48){V˜(r−1)}−1W+λ G˜]−1W(cid:48){V˜(r−1)}−1.
i i i i β i i
Then the updates of σ2(r) and D(r) =diag(σ2(r),D(r)) are given by
(cid:15) ξ a
n (cid:110) (cid:111)
σ2(r) = N−1 (cid:80) {ˆ(cid:15)(r)}(cid:48)ˆ(cid:15)(r)+σ2(r−1)[m −σ2(r−1)trace(H(r−1))] ,
(cid:15) i i (cid:15) i (cid:15) i
i=1
D(r) = n−1 (cid:80)n (cid:110)ξˆ(r){ξˆ(r)}(cid:48)+[D(r−1)−D(r−1)Z(cid:48)H(r−1)Z D(r−1)](cid:111).
ξ i i ξ ξ i i i ξ
i=1
Step 3. Repeat Steps 2 and 3, until some convergence conditions are satis-
fied.
2.3. Smoothing Parameter Selection
Thesmoothnessofβ(t)andb (t),i=1,...,n,arecontrolledbythesmooth-
i
ing parameter λ and λ , respectively. Define λ = (λ ,λ )(cid:48), then the optimal
β b β b
7
value for λ is chosen by minimizing the generalized cross-validation (GCV) cri-
terion defined as follows
SSE(λ)
GCV(λ)= ,
(N −df(λ))2
n
where SSE(λ) = (cid:80)(cid:107)Yi − Wiθ(cid:98) − Zi(cid:98)ξi(cid:107)2, N = (cid:80)ni=1mi, and df(λ) is the
i=1
effective degrees of freedom, which is calculated as df(λ)=trace(Q), where Q
is given by
    −1 
1 W(cid:48)W W(cid:48)Z λβG˜ 0 W(cid:48)
Q=(W,Z)  +   .
σ(cid:15)2 Z(cid:48)W Z(cid:48)Z 0 In(cid:78)D(cid:101)ξ Z(cid:48)
Define the matrix S = W [W(cid:48)V˜−1W+λ G˜]−1W(cid:48)V˜−1. The predictor Yˆ =
i i β i
W θˆ+Z ξˆ can then be expressed as Yˆ =Q Y with
i i i i i
Qi =Si+ZiD(cid:101)ξZ(cid:48)iV˜−1(Imi −Si).
Note that the smooth matrix Q=(Q(cid:48),...,Q(cid:48) )(cid:48).
1 n
2.4. Constructing the Confidence Intervals
To construct the confidence intervals of α and the point-wise confidence
0
bands of β(t), we need to calculate the covariance matrix of θˆ:
(cid:32) n (cid:33)−1(cid:32) n (cid:33)
Cov(θˆ) = (cid:88)Wi(cid:48)V˜i−1Wi+λβG(cid:101) (cid:88)Wi(cid:48)V˜i−1Cov(Yi)V˜i−1Wi
i=1 i=1
(cid:32) n (cid:33)−1
(cid:88)Wi(cid:48)V˜i−1Wi+λβG(cid:101) . (7)
i=1
In (7), Cov(Y ) can be replaced by V˜ to account for our roughness penalty on
i i
b (t),i=1,...,n. For simplicity, instead of using Cov(θˆ), we use
i
(cid:32) n (cid:33)−1
Cov(θˆ)= (cid:88)Wi(cid:48)V˜i−1Wi+λβG(cid:101) . (8)
i=1
Let C(cid:100)ov(θˆ) be the estimator of the covariance matrix (8) and partition it as
 
σˆ2 Σˆ
C(cid:100)ov(θˆ)= Σˆ1(cid:48)1 Σˆ12 . Then the 95% confidence intervals of α0 is approx-
12 22
imately as
(αˆ −1.96σˆ , αˆ +1.96σˆ ),
0 11 0 11
8
and the 95% pointwise bands of β(t) can be approximately given by
(cid:18) (cid:113) (cid:113) (cid:19)
βˆ(t)−1.96 V(cid:100)ar[βˆ(t)], βˆ(t)+1.96 V(cid:100)ar[βˆ(t)] , for all t∈S,
where V(cid:100)ar[βˆ(t)]=Φ(cid:48)(t)Σˆ22Φ(t). Moreover, the estimate of γ(s,t) can be given
as
γˆ(s,t)=Ψ(cid:48)(s)D(cid:98)bΨ(t),
(cid:16) (cid:17)−1
where we use D(cid:98)b = (0,IK) D(cid:98)−ξ1+λbGξ (0,IK)(cid:48) instead of D(cid:98) in order to
account for our roughness penalty on b (t),i=1,...,n in our method.
i
2.5. Reconstructing the predictors X (t)
ij
When covariates X (t) in the functional linear mixed-effects model (2) are
ij
not be exactly observable but measured with errors, the estimators and infer-
ence may be biased if one ignores these measurement errors. Hence, we need to
adjusttheresultingbias. Inthisarticle,wesuggesttoreconstructthefunctional
predictorsX (t)byusingalargenumberoffunctionalprincipalcomponentsob-
ij
tained from a smooth estimator of the covariance matrix estimator (Goldsmith
et al., 2012) firstly. Then, we treat the estimated Xˆ (t) as the true predictors
ij
and applying the REML-based EM algorithm.
Define the covariance function of X(t) as
C(s,t)=Cov(X(t),X(s)).
Mercer’s theorem (Ash and Gardner, 1975) states that C(s,t) has the eigen-
decomposition
∞
(cid:88)
C(s,t)= λ ϕ (s)ϕ (t),
k k k
k=1
where λ ≥ λ ≥ ... ≥ 0 satisfying (cid:80)∞ λ < ∞, and ϕ (t)’s form a complete
1 2 k=1 k k
orthonormal basis in S×S. Then, X(t) allows the Karhunen-Loeve decomposi-
tion (Rice and Silverman, 1991)
∞
(cid:88)
X(t)=µ(t)+ ξ ϕ (t)
k k
k=1
9
whereϕ (·)istheorthonormaleigenfunction,whichisalsocalledthefunctional
k
principalcomponent(FPC).Thecoefficientsξ iscalledtheFPCscoreofX(t),
k
which satisfies E(ξ )=0, E(ξ2)=λ , and E(ξ ξ )=0 for k (cid:54)=l.
k k k k l
Suppose we have the following additive measurement error model,
W (t)=X (t)+e (t),
ij ij ij
whereW (t)istheobservedvalue,X (t)istheunderlyingtruevaluefortheith
ij ij
subject at the time point t, and e (t) represents the measurement error at the
ij
time point t. We assume that e (t) is a mean zero process, and {X (t),e (t)}
ij ij ij
are mutually independent. We estimate C(s,t) by using a method-of-moments
approach, and then smooth the off-diagonal elements of this observed covari-
ance matrix to remove the ‘nugget effect’ that is caused by measurement error
(Staniswalis and Lee, 1998; Yao et al., 2005; Goldsmith et al., 2012).
We use the principal analysis by conditional estimation (PACE) algorithm
proposedbyYaoetal.(2005)toestimatethemeancurveµ (t),theFPCsϕ (t)
i ik
and the FPC scores ξ from the observations W (t ). Let µˆ (t), ϕˆ (t), and
ijk ij ijk i ik
ξˆ bethecorrespondingestimatorsofµ (t),ϕ (t)andξ ,respectively. Then
ijk i ik ijk
an estimate of X (t) is obtained as
ij
M
Xˆ (t)=µˆ (t)+(cid:88)ξˆ ϕˆ (t)
ij i ijk ik
k=1
wherethenumberofFPCs,M,canbechosenbyAIC,BIC,thecross-validation
method, or the empirical experience based on the percentage of explained vari-
ance (such as 90% or 95%).
3. Simulation studies
In this section, we perform some numerical experiments to assess the effi-
ciency of our proposed estimating procedure for the functional linear mixed-
effects model (2). The performance of our estimation method is evaluated by
the following relative mean integrated square error (RMISE) for the estimated
10