Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2915799

Formats

Article sections

- Abstract
- 1. INTRODUCTION
- 2. ALGORITHM
- 3. SIMULATION DETAILS AND RESULTS
- 4. EMPIRICAL EXAMPLES
- 5. DISCUSSION
- REFERENCES

Authors

Related links

Can J Stat. Author manuscript; available in PMC 2010 August 4.

Published in final edited form as:

Can J Stat. 2010 June; 38(2): 256–270.

doi: 10.1002/cjs.10062PMCID: PMC2915799

NIHMSID: NIHMS184162

Department of Statistics Texas A&M University College Station, Texas USA, 77843-3143

Josue G. Martinez: ude.umat.tats@zenitramgj

Faming Liang: ude.umat.tats@gnailf

Lan Zhou: ude.umat.tats@uohzl

Raymond J. Carroll: ude.umat.tats@llorrac

See other articles in PMC that cite the published article.

The authors consider the analysis of hierarchical longitudinal functional data based upon a functional principal components approach. In contrast to standard frequentist approaches to selecting the number of principal components, the authors do model averaging using a Bayesian formulation. A relatively straightforward reversible jump Markov Chain Monte Carlo formulation has poor mixing properties and in simulated data often becomes trapped at the wrong number of principal components. In order to overcome this, the authors show how to apply Stochastic Approximation Monte Carlo (SAMC) to this problem, a method that has the potential to explore the entire space and does not become trapped in local extrema. The combination of reversible jump methods and SAMC in hierarchical longitudinal functional data is simplified by a polar coordinate representation of the principal components. The approach is easy to implement and does well in simulated data in determining the distribution of the number of principal components, and in terms of its frequentist estimation properties. Empirical applications are also presented.

Because longitudinal data can be viewed as sparsely observed functional data (Rice, 2004), functional principal components analysis becomes an important tool for analyzing longitudinal data. Since the early contributions of Besse & Ramsay (1986), Ramsay & Dalzell (1991), and Rice & Silverman (1991), a substantial literature has developed for functional principal components analysis, including but not limited to the work of Silverman (1996), Cardot (2000), James, Hastie & Sugar (2000), Rice & Wu (2001), Boente & Fraiman (2002), Ramsay & Silverman (2005, Chapters 8–10), Yao, Müller & Wang (2005a,b), Hall & Hosseini-Nasab (2006), Hall, Müller & Wang (2006), Jank & Shmueli (2006), Ocaña, Aguilera & Escabias (2007), Reiss & Ogden (2007), and Huang, Shen & Buja (2008).

Let *Y*_{i} be the ^{th} response for the *i ^{th}* subject at time

$${Y}_{i\ell}=\mu \left({t}_{i\ell}\right)+{g}_{i}\left({t}_{i\ell}\right)+{\u03f5}_{i\ell ,}$$

(1)

where μ(·) is the fixed effects population-level function, *g _{i}*(·) represents random functions associated with each individual, and ϵ

It is standard in the functional data analysis literature to model the correlations among responses at different times through a mixed effects model, e.g., via the Karhunen-Lóeve expansion used in a number of the papers referenced above. As a modification of work by James et al. (2000) and Rice & Wu (2001), to fit models such as (1), Zhou, Huang & Carroll (2008) introduced penalized spline estimation of functional principal components for longitudinal data using a mixed effects model framework. Suppose that there are *i* = 1, …, *n* subjects, each of which contribute = 1, …, *m _{i}* observations at times or locations

$$\begin{array}{cc}\hfill {Y}_{i}& ={B}_{i}{\theta}_{\mu}+{B}_{i}{\Theta}_{f}{\alpha}_{i}+{\u03f5}_{i};\hfill \\ \hfill E({\u03f5}_{i}\mid {B}_{i})& =0;\hfill \\ \hfill \mathrm{cov}({\u03f5}_{i}\mid {B}_{i})& ={\sigma}_{\u03f5}^{2}{I}_{{m}_{i}};\hfill \\ \hfill {\alpha}_{i}& =\mathrm{Normal}\{0,{D}_{\alpha}=\mathrm{diag}({d}_{1},\dots ,{d}_{M})\},\hfill \end{array}$$

(2)

where the principal components matrix Θ_{f} is *p* × *M* with *p* > *M*, the main effects parameter θ_{μ} is *p* × 1, and the vector of principal components scores α_{i} is *M* × 1. In addition to the orthogonality constraint on the basis functions, for identifiability we need Θ_{f} to be orthogonal, ${\Theta}_{f}^{\mathrm{T}}\Theta ={I}_{M}$. Thus, comparing (1) and (2), we are making the approximations that μ(*t*) ≈ *b*^{T}(*t*)θ_{μ} and *g _{i}*(

It is interesting to note that in common with the Karhunen-Lóeve approaches, as in ours, if the number of principal components is *M* = 1, then the correlation structure of responses at different times is equicorrelated, but this is not the case when *M* > 1.

As in any kind of spline-based analysis, it is useful to think about penalizing the likelihood to ensure smoothness of fit and to avoid overfitting. In model (2), we have two functions that need to be penalized: the fixed effects functions represented as *b*^{T}(*t*)θ_{μ} and the random effects functions represented as *b*^{T}(*t*)Θ_{f}. For the fixed effects function, for example, the usual device is to have a penalty parameter, here called λ_{μ}, and to penalize the likelihood by ${\lambda}_{\mu}{\theta}_{\mu}^{\mathrm{T}}K{\theta}_{\mu}$, where *K* is a so-called penalty matrix. The typical form of the penalty matrix is the integrated square of the second derivative matrix, i.e., *K* = ∫ *b*″(*t*)*b*″(*t*)^{T}*dt*. Methods for constructing the orthogonal basis functions *B _{i}* and for computing the penalty matrix

Enforcing smoothing of the random effects functions is somewhat more involved, because *b*^{T}(*t*)Θ_{f} has *M* columns, each of which needs to be penalized. Here is the device used by Zhou, et al. (2008). If _{j} is the *j ^{th}* column of Θ

Of course, penalty parameters need to be estimated, and for this purpose there are a variety of possibilities. Zhou, et al. (2008) used 10-fold crossvalidation to estimate (λ_{μ}, λ_{f}) for any given *M*. To choose the number of principal components, they use an informal ad hoc procedure, see their Section 5.2. James et al. (2000) used the informal %-variation explained and a likelihood ratio test, although they found that the performance of the latter was “ambiguous” and they recommended the use of the former.

Especially when the number of subjects is small, there will be ambiguity as to the choice of the number of principal components. We propose to acknowledge this ambiguity by taking a Bayesian approach to the problem.

An outline of this paper is as follows. In Section 2 we present the basic algorithm. Section 3 describes simulation studies, while Section 4 investigates empirical examples. Concluding remarks are presented in Section 5. Technical details are given in an appendix.

The model is as described at (2).

In what follows, let ${\sigma}_{\mu}^{2}={\lambda}_{\mu}^{-1}$ and ${\sigma}_{f}^{2}={\lambda}_{f}^{-1}$. It is convenient to work with the parameterization ξ_{1} = log(σ_{ϵ}), ξ_{2} = log(σ_{μ}) and **ξ**_{3} = log{diag(*D*_{α})}. We write the number of principal components as *M*.

As described in Section 2.2, because of the restriction that ${\Theta}_{f}^{\mathrm{T}}{\Theta}_{f}=I$, we will parameterize Θ_{f} in terms of auxiliary parameters Ω_{f}. If α = (α_{1}, …, α_{n}), the main effects vector is θ_{μ} and ** Y** denotes the response data, then except for constants of proportionality, we write the posterior distribution of (Ω

In formulating this problem, we face a number of challenges.

- By definition, Θ
_{f}is constrained to be orthogonal. We achieve this orthogonality by combining the polar coordinate approach of Zhang and Davidian (2001) with a Gram-Schmidt transformation, see Section 2.2. - Even though Θ
_{f}is orthogonal by construction, it still needs to be smoothed. - Because the models are of different dimensions, a basic reversible jump type algorithm is attractive (Green, 1995). However, we observed difficulty in accepting the death step of such an algorithm, and in order to explore the model space chose instead to implement the Stochastic Approximation Monte Carlo or SAMC algorithm of Liang, Liu & Carroll (2007).

In this section, we show how to write Θ_{f} as a function of other parameters Ω_{f}, in such a way that Θ_{f} is orthogonal, and all elements of its first row are positive, as required by the identifiability result in Zhou et al. (2008).

Because Θ_{f} is an orthogonal *p* × *M* matrix with *p* > *M*, we turn to a polar coordinate representation, handily provided by Zhang and Davidian (2001). Let Θ_{f} = (_{1}, …, _{M}) where *p* > *M* + 1, each _{j} is *p* × 1, and define the *p* × 1 vector *s*_{j} as follows. Its first element is sin(γ_{1j}), and its *k*^{th} element for *k* = 2, …, *p* − 1 is $\mathrm{sin}\left({\gamma}_{kj}\right){\prod}_{\ell =1}^{k-1}\mathrm{cos}\left({\gamma}_{\ell j}\right)$. Finally, its *p*^{th} element is ${\prod}_{\ell =1}^{p-1}\mathrm{cos}\left({\gamma}_{\ell j}\right)$. These are thus polar coordinates in *p* − 1 parameters. Take −π/2 ≤ γ_{kj} ≤ π/2 for all (*k*, *j*) except *k* = *j* = 1, and let 0 ≤ γ_{11} ≤ π/2 to force identifiability, see Zhou et al. (2008) for details. The matrix *s _{j}* is of unit length: ${s}_{j}^{\mathrm{T}}{s}_{j}={\parallel {s}_{j}\parallel}^{2}=1$. Define ${\gamma}_{j}^{\mathrm{T}}=({\gamma}_{1j},\dots ,{\gamma}_{(p-1)j})$ and Ω

To form an orthogonal *p*×*M* matrix, we use a Gram-Schmidt transformation. Define proj_{e}(*s*) = (*s*^{T}*e*)*e*. Then _{1} = *s*_{1}, for *k* = 2, …,*M*, define ${\vartheta}_{k\ast}=(I-{\sum}_{j=1}^{k-1}{\vartheta}_{j}{\vartheta}_{j}^{\mathrm{T}}){s}_{k}$, and _{k} = _{k*}/_{k*}. Finally, take Θ_{f} = (_{1}, …, _{M}).

Write *e*^{T} = (1, 0, …0). Note that the first element of _{k*}, *e*^{T}_{k*} is ${e}^{\mathrm{T}}(I-{\sum}_{j=1}^{k-1}{\vartheta}_{j}{\vartheta}_{j}^{\mathrm{T}}){s}_{k}$. To force the first element of _{k} to be positive, we set _{k} = sign(*e*^{T}_{k*})_{k*}/_{k*}.

Let *M* denote the number of principal components, and assume that it has a uniform prior distribution. We marginalized the posterior by integrating out θ_{μ} and α, resulting in the function *g*(Ω_{f}, **ξ**_{3}, ξ_{2}, ξ_{1}, ${\sigma}_{f}^{2}$ |** Y**) given in equation (6) in the Appendix. Let (${\Omega}_{f}^{\left(v\right)}$, ${\xi}_{3}^{\left(v\right)}$, ${\xi}_{2}^{\left(v\right)}$, ${\xi}_{1}^{\left(v\right)}$, ${\sigma}_{f}^{2\left(v\right)}$) denote the sample obtained at iteration

In this step, we decide which component of the model to update at iteration *v* + 1. We let “MH” here denote a Metropolis-Hastings step, see Section 2.3.2. Birth and death steps are described in Sections 2.3.3 and 2.3.4, respectively.

- Draw
*U*= Uniform(0, 1). - If
*M*= 1 and*U*< 0.5 do “birth” step and set ω_{12}= 0.5; if*M*= 1 and*U*≥ 0.5 do “MH” step. - If
*M*= 2 and*U*< 1/3 do “birth” step and set ω_{23}= 1/3; If*M*= 2 and 1/3 ≤*U*< 2/3 do “death” step and set ω_{21}= 1/3; if*M*= 2 and*U*≥ 2/3 do “MH” step. - If
*M*= 3 and*U*< 0.5 do “death” step and set ω_{32}= 0.5; if*M*= 3 and*U*≥ 0.5 do “MH” step.

In this step, we keep the dimensions of ${\Omega}_{f}^{\left(v\right)}$ and ${\xi}_{3}^{\left(v\right)}$ unchanged while the parameters are updated. Draw *U* = Uniform(0, 1). If *U* ≤ 1/5 make Move-1; if 1/5 < *U* ≤ 2/5, make Move-2; if 2/5 < *U* ≤ 3/5, make Move-3; if 3/5 < *U* ≤ 4/5, make Move-4, and if 4/5 < *U* ≤ 1 make Move-5. In what follows, we set τ_{1} = τ_{2} = τ_{3} = τ_{4} = 0.1.

- (a)(Move-1) Updating ${\xi}_{1}^{\left(v\right)}$ and ${\xi}_{2}^{\left(v\right)}$. Let
*U*_{1}and*U*_{2}be independent Uniform(−0.5, 0.5) random variables. Draw ${\xi}_{1}^{\ast}={\xi}_{1}^{\left(v\right)}+0.2{U}_{1}$ and ${\xi}_{2}^{\ast}={\xi}_{2}^{\left(v\right)}+0.5{U}_{2}$. Calculate the MH ratio$$r=\frac{f({\Omega}_{f}^{\left(v\right)},{\xi}_{3}^{\left(v\right)},{\xi}_{2}^{\ast},{\xi}_{1}^{\ast},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})}{f({\Omega}_{f}^{\left(v\right)},{\xi}_{3}^{\left(v\right)},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})}.$$Accept (${\xi}_{1}^{\ast}$, ${\xi}_{2}^{\ast}$) with probability min(1,*r*). - (b)(Move-2) Updating ${\Omega}_{f}^{\left(v\right)}$. Draw a vector $e=\mathrm{Normal}(0,{\tau}_{2}^{2})I$ and add
*e*to a column of ${\Omega}_{f}^{\left(v\right)}$ selected at random. Denote the new Ω_{f}matrix by ${\Omega}_{f}^{\ast}$. Calculate the MH ratiowhere $\mathcal{A}$ denotes the set ${(-\pi \u22152,\pi \u22152)}^{\mathit{dim}\left({\Omega}_{f}^{\ast}\right)}$. Accept ${\Omega}_{f}^{\ast}$ with probability min(1,$$r=\frac{f({\Omega}_{f}^{\ast},{\xi}_{3}^{\left(v\right)},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})\mathrm{I}({\Omega}_{f}^{\ast}\in \mathcal{A})}{f({\Omega}_{f}^{\left(v\right)},{\xi}_{3}^{\left(v\right)},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})},$$*r*). - (c)(Move-3) Updating ${\xi}_{3}^{\left(v\right)}$. Let
*U*_{3}= Uniform(−0.5, 0.5). Draw the*j*^{th}component of ${\xi}_{3}^{\ast}$, ${\xi}_{3,j}^{\ast}={\xi}_{3,j}^{\left(v\right)}+0.1{U}_{3}$. Calculate the MH ratio$$r=\frac{f({\Omega}_{f}^{\left(v\right)},{\xi}_{3}^{\ast},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})}{f({\Omega}_{f}^{\left(v\right)},{\xi}_{3}^{\left(v\right)},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})}.$$Accept ${\xi}_{3}^{\ast}$ with probability min(1,*r*). - (d)(Move-4) Updating ${\sigma}_{f}^{2\left(v\right)}$. This can be done by a Gibbs step, see the appendix.
- (e)(Move-5) Put the elements of ${\xi}_{3}^{\left(v\right)}$ is descending order and exchange columns of ${\Omega}_{f}^{\left(v\right)}$ and ${\Theta}_{f}^{\left(v\right)}$ accordingly. Because identifiability requires that the elements of
*D*_{α}= diag{exp(ξ_{3})} be in descending order, this step avoids the label switching problem. Accept ${\Omega}_{f}^{\ast}$ and ${\xi}_{3}^{\ast}$ with probability min(1,*r*).

In this step, we try to append a new column to ${\Omega}_{f}^{\left(v\right)}$ and append an element to ${\xi}_{3}^{\left(v\right)}$.

- (a)(Angle sampling) Draw a random vector γ from Uniform$(-\frac{\pi}{2},\frac{\pi}{2})$, and draw ${\xi}^{\ast}=\mathrm{Normal}(0,{\tau}_{4}^{2})$. set ${\Omega}_{f}^{\ast}=({\Omega}_{f}^{\left(v\right)},\gamma )$ and ${\xi}_{3}^{\ast}=\mathrm{diag}[{\xi}_{3}^{\left(v\right)},{\xi}^{\ast}]$.
- (b)(Proposal acceptance) Calculate the MH ratio$$r=\frac{f({\Omega}_{f}^{\ast},{\xi}_{3}^{\ast},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y}+1)}{f({\Omega}_{f}^{\left(v\right)},{\xi}_{3}^{\left(v\right)},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})}\frac{{\tau}_{4}}{{\pi}^{-(p-1)}\varphi ({\xi}^{\ast}\u2215{\tau}_{4})}\frac{{\omega}_{k+1,k}}{{\omega}_{k,k+1}}.$$(3)
- (c)Accept ${\Omega}_{f}^{\ast}$ and ${\xi}_{3}^{\ast}$ with probability min(1,
*r*).

In this step, we try to delete a random column of ${\Omega}_{f}^{\left(v\right)}$ and the corresponding diagonal element of ${\xi}_{3}^{\left(v\right)}$.

$$r=\frac{f({\Omega}_{f}^{\ast},{\xi}_{3}^{\ast},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y}-1)}{f({\Omega}_{f}^{\left(v\right)},{\xi}_{3}^{\left(v\right)},{\xi}_{2}^{\left(v\right)},{\xi}_{1}^{\left(v\right)},{\sigma}_{f}^{2\left(v\right)}\mid \mathit{Y})}\frac{{\pi}^{-(p-1)}\varphi ({\xi}^{\left(v\right)}\u2215{\tau}_{4})}{{\tau}_{4}}\frac{{\omega}_{k-1,k}}{{\omega}_{k,k-1}},$$

(4)

where ξ^{(v)} denotes the selected diagonal element of ${\xi}_{3}^{\left(v\right)}$.

The two groups of parameters, (Ω_{f}, ξ_{3}) and (${\xi}_{1},{\xi}_{2},{\sigma}_{f}^{2}$), are updated in an unbalanced way in the above algorithm. This is allowed by the Gibbs sampler.

Because of low acceptance rates in the death step of the algorithm in Section 2.3, we opt to implement the Stochastic Approximation Monte Carlo or SAMC algorithm (Liang et al., 2007, Section 5; Liang, 2009), which is a method designed to provide samples from the entire space that are then weighted to form posterior distributions. The SAMC algorithm, especially appropriate for our problem, is given in equation (13) of Liang et al. (2007). We describe the algorithm for three principal components, detailed reasoning for this is provided later.

Let *v* be the index in the MCMC calculations. Let ζ = (1/3, 1/3, 1/3) and let *v*_{0} be a tuning parameter that can be adjusted. Define κ_{v} = *v*_{0}/ max(*v*_{0}, *v*). Also define θ^{(v)} = (θ_{v,1}, θ_{v,2}, θ_{v,3})^{T}, where θ_{v,i} is the current estimate of the *i*^{th} normalizing constant of the target density, which in this case equals log(1/ζ_{i}). Let $\mathcal{T}$ be the sampling space of allowable values of θ^{(v)}. Liang et al. (2007) make this a very large space, but because of the exponentiation involved in the adjusted Metropolis ratio, we have restricted this to the cube [−10, 10]^{3}.

First start with *v* = 1 and θ^{(v)} = {log(3), log(3), log(3)}^{T}. Let *M*^{(v)} be the current number of principal components, and let *M** be a proposal for the number of principal components. If *M*^{(v)} = *M**, the Metropolis steps in Section 2.3.2 do not change. If *M*^{(v)} ≠ *M**, and if *r* is the birth or death ratio in (3) or (4), respectively, the only thing that changes is that *r* is multiplied by exp(θ_{v,M(v)} − θ_{v,M*}). Let the new model be *M*^{(v+1)}. Set θ* = θ^{(v)} + κ_{v+1}(*e*_{M(v+1)} − ζ), where the three component vector, *e*_{M}^{(v+1)}, is all zeros except with a 1.0 in location *M*^{(v+1)}. Then if ${\theta}^{\ast}\in \mathcal{T}$, we set θ^{(v+1)} = θ*. Otherwise, we set θ^{(v+1)} = θ* + *c**, where *c** is chosen so that ${\theta}^{\ast}+{c}^{\ast}\in \mathcal{T}$.

Having run SAMC, we are now in a position to do posterior inference. Thus for example, assuming equal prior probabilities for the three possible models, the posterior probability that the number of principal components, *M* = *j*, is given as

$$\mathrm{pr}(M=j\mid \mathit{Y})=\sum _{v}^{T}\mathrm{exp}\left({\theta}_{v,{M}^{\left(v\right)}}\right)I({M}^{\left(v\right)}=j)\u2215\sum _{v}^{T}\mathrm{exp}\left({\theta}_{v,{M}^{\left(v\right)}}\right),$$

where *T* is the number of iterations used to do posterior inference. Posterior quantities can be calculated similarly. For example, given that we have *j* principal components, let ${\sigma}_{\mu}^{\left(v\right)}$ be the value of σ_{μ} at iteration *v*. Then

$$\begin{array}{cc}\hfill E({\sigma}_{\mu}\mid M=j,\mathit{Y})& =\sum _{v}^{T}{\sigma}_{\mu}^{\left(v\right)}\mathrm{exp}\left({\theta}_{v,{M}^{\left(v\right)}}\right)I({M}^{\left(v\right)}=j)\u2215\sum _{v}^{T}\mathrm{exp}\left({\theta}_{v,{M}^{\left(v\right)}}\right)I({M}^{\left(v\right)}=j);\hfill \\ \hfill E({\sigma}_{\mu}\mid \mathit{Y})& =\sum _{v}^{T}{\sigma}_{\mu}^{\left(v\right)}\mathrm{exp}\left({\theta}_{v,{M}^{\left(v\right)}}\right)\u2215\sum _{v}^{T}\mathrm{exp}\left({\theta}_{v,{M}^{\left(v\right)}}\right).\hfill \end{array}$$

(5)

From (5), it can be seen that our approach employs model averaging to do posterior inference.

We provide details only for the case of up to 3 principal components because, in biological data, 3 principal components usually suffice, and most algorithms become unstable for higher number of principal components, see the comments of Rice and Silverman (1991) and Hall, Müller & Wang (2006) in a different context. If more than 3 principal components are used in our simulations and data analyses, the 4^{th} and higher diagonal elements of the variance matrix *D*_{α}, which are in decreasing order, become crowded around zero, making the corresponding elements of Θ_{f} essentially unidentified.

We simulated 200 data sets from model (2) with *n* = 20 subjects, *m _{i}*

We ran SAMC with 200,000 Gibbs steps, and repeated the process 5 times with different starting values. In each case, we started by assuming 1 principal component, with Ω_{f} generated completely at random, and with starting values for ξ_{1} = log(σ_{ϵ}) = log{Uniform(0.4, 0.6)}, ξ_{3} = Uniform(0.5, 1.5), ξ_{2} = log(σ_{μ}) = Uniform[4, 6] and ${\sigma}_{f}^{2}=2c$, where *c* yields 4.5 degrees of freedom for the smoother matrix trace{(*B*^{T}*B* +*K*/*c*)^{−1}*B*^{T}*B*}. We also repeated the analysis with no smoothing: as might be expected, there was slightly more variability in this latter case.

In all cases, ${\sigma}_{\u03f5}^{2}$ and θ_{μ} were estimated with little bias. In addition, the within-subject function variances diag{cov(*B*_{i}Θ_{f}α_{i})} were estimated with only minor bias. An example of this is given in Figure 1, for the case that *D*_{α} = diag(1.50, 0.75). The top left panel shows the true function and the 5^{th} and 95^{th} pointwise percentiles of the simulations: there is almost no bias so that the mean of the simulations is not displayed. The bottom right panel shows that the within-subject function variances are also essentially unbiased. The top right panel gives the posterior probabilities of the models: in this case, 2 principal components were selected for all 200 data sets. The bottom left panel gives 16 examples of the within-subject function variances.

Results of the simulation study in the case of two principal components with *D*α = diag(1.50; 0.75), as described in Section 3.1. Top left: all 200 posterior means for *B*_{i}θ_{μ}. The true value is within the 200 lines. Top right: the **...**

We also compared our analysis with the method of Zhou et al. (2008). Their informal method choose 2 principal components in 90.5% of the cases, and 1 principal component in 7% of the cases. In Table 1, we compare the mean absolute error and the mean squared error for estimating the within-subject function variances for the method of Zhou et al. (2008) versus our method. SAMC had somewhat smaller mean squared errors for estimating the mean function, and much smaller mean squared errors for estimating the variance function.

In this case, we simulated with ${\sigma}_{\u03f5}^{2}=0.25$, and in two settings with *D*_{α} = 1.00 and *D*_{α} = diag(1.50, 0.75). We used the Bspline basis functions as in Zhou et al. (2008) with 7 basis functions. In the case of one principal component, the model is

$${Y}_{i\ell}=\mathrm{sin}\left(2\pi {t}_{\ell}\right)+\sqrt{2}\mathrm{cos}\left(2\pi {t}_{\ell}\right){\alpha}_{i}+{\u03f5}_{i\ell},$$

while with two principal components, and α_{i} = (α_{i1}, α_{i2})^{T},

$${Y}_{i\ell}=\mathrm{sin}\left(2\pi {t}_{\ell}\right)+\sqrt{2}\mathrm{cos}\left(2\pi {t}_{\ell}\right){\alpha}_{i1}+\{11.65{(x-0.5)}^{2}-\mathrm{cos}\left(2\pi {t}_{\ell}\right)\u2215{\pi}^{2}\}{\alpha}_{i2}+{\u03f5}_{i\ell}.$$

We used 7 basis functions with 3 interior knots. All functions are adequately but not exactly approximated with this number of knots. The same priors were used as above, except that we centered them at the estimates from Zhou et al. (2008). We also experimented with fixed priors as in Section 3.2, with little change in results.

The results are again reported in Table 1, with graphical representation in the two principal components case given in Figure 2. In the former, we see that SAMC very similar to the method of Zhou et al. (2008) for estimating the mean function sin(2π*t _{j}*), and much better in terms of estimating the variance function. Figure 2 shows that there is more ambiguity in the posterior model probabilities across the simulated data sets, although in 97% of the cases two principal components have posterior probability > 0.50. We actually think that this ambiguity is a strength of our approach. If the model is not correct, then there should be ambiguity in the choice of the number of principal components, with some random data sets coming close to being fit by 2 principal components, and others needing 3.

The simulations presented in Section 3 suggest that the method of Zhou et al. (2008) and our SAMC method will have roughly equal performance in terms of estimating the mean function, but will tend to have different performance when estimating the variance function. The SAMC method of course comes endowed with Bayesian posterior credible intervals. In this section, we will describe two data sets, each with two responses, and then display the analysis of all four. All analyses used 7 basis functions: the estimates of ${\sigma}_{\u03f5}^{2}$, ${\sigma}_{\mu}^{2}=1\u2215{\lambda}_{\mu}$, ${\sigma}_{f}^{2}=1\u2215{\lambda}_{f}$ and the first component of *D*_{α} from Zhou et al. (2008) were used as starting values. The starting value of Θ_{f} was random. SAMC was run with 200,000 iterations twice, and the results pooled.

Our first example concerns data from an AIDS clinical trial. Both virologic and immunologic surrogate markers such as plasma HIV RNA copies (viral load) and CD4+ cell counts currently play important roles in evaluating antiviral therapies in AIDS clinical research. Both CD4+ cell counts and plasma HIV RNA copies have each been used as sole surrogate markers in AIDS clinical trials. Here we study the time course of plasma HIV RNA copies and CD4+ cell counts obtained from an AIDS clinical study conducted by the AIDS Clinical Trials Group (ACTG 315). In this study, forty-six evaluable HIV-1 infected patients were treated with potent antiviral therapy consisting of ritonavir, 3TC and AZT, see Lederman, Connick, Landay, Kuritzkes, Spritzler, St Clair, Kotzin, Fox, Chiozzi, Leonard, Rousseau, Wade, Roe, Martinez & Kessler (1998) and Wu & Ding (1999) for more details on this study. After initiation of potent antiviral treatment at day 0, patients were followed for up to 10 visits and at each visit both viral load and CD4+ cell counts were monitored simultaneously. The actual visit times are irregularly spaced and different for different patients. The visit time varies from day 0 (first visit) to day 196. We took logarithms for both CD4 data and RNA data and then standardized them to have mean zero and variance one. We also standardized follow up times to the unit interval. In this case, what we have been calling “time”, *t*_{i}, is indeed time.

We apply our method to phenotypic longitudinal data of colonic crypts in the rat model. Colonic crypts are tube like structures in the colon formed by cells which undergo development from stem cells to fully differentiated mature colonic cells. This development begins at the bottom of the colonic crypt where the mother stem cells divide and produce new colon cells which will mature as they travel up the crypt to the colon surface. Colonic DNA damage can change the behavior of any cell by fundamentally altering its normal function and potentially inducing uncontrollable growth. This can contribute to carcinogenesis, see Hong, Bancroft, Turner, Davidson, Murphy, Carroll, Chapkin & Lupton (2005).

In this model, the times *t*_{i} correspond to the cell locations, normalized to the unit interval. The responses are the DNA adduct levels, a measure of DNA damage, in the proximal and in the distal region of the colon. In this analysis, there were 30 subjects. The DNA adduct levels in the proximal and distal regions were log-transformed, and then standardized to have mean zero and variance one.

In this case, what we call “time” is not time, but actually location of the cell. Since stem cells, or colloquially “baby” cells are at the bottom of the crypts, with *t* near zero, and since mature cells at the top of the crypts, *t* near one, the functions may be thought to be functions of the age of the cells.

The method of Zhou et al. (2008) picked 1 principal component in all cases. The SAMC method had posterior probabilities of 1 principal component equal to 1.0 for the CD4 and RNA data. However, for the distal adducts, the posterior probability of two principal components was 0.92, while for proximal data 1 principal component has a posterior probability of 0.60. In all cases, the posterior probability of 3 principal components equalled 0.00.

As expected by the simulations, the estimated posterior mean function *b*^{T}(*t*)θ_{μ} was almost identical to that found by Zhou et al. (2008) (not displayed). In Figure 3, we display the posterior mean and medians of *b*^{T}(*t*)θ_{μ}, along with 90% pointwise credible intervals. In the adduct level, the distribution of *t*_{i} was almost uniform, and hence the lengths of the credible intervals might be expected to be nearly constant, a conjecture seen in the plots. On the other hand, in the CD4 and RNA data, the distribution of *t*_{i} is highly skewed right, with most of the observations being between 0.0 and 0.4. In this case, one would expect the credible intervals to be much wider for larger values of *t*_{i} than for smaller values, an expectation born out in the plot.

Results of the data analysis of Section 4. Displayed are the posterior mean function (dash-dotted magenta line), the posterior median (black solid line) and 90% pointwise posterior credible intervals (blue and red dashed lines). The labels for each graph **...**

Figure 4 displays results for estimating the marginal variances ${b}^{\mathrm{T}}\left(t\right){\Theta}_{f}{D}_{\alpha}{\Theta}_{f}^{\mathrm{T}}b\left(t\right)$ (magenta dot-dashed lines), the estimate from Zhou et al. (2008) (black solid lines) and posterior pointwise 90% credible intervals. Variance functions are of course much harder to estimate than mean functions, and as seen in the simulations here we do not expect SAMC and the method of Zhou et al. (2008) to be nearly identical, and they are not, although the differences are fairly small given the uncertainties.

Results of the data analysis of Section 4. Displayed are the posterior mean variance function (dash-dotted magenta line), the estimated variance function from the method of Zhou et al. (2008) and labeled “pdfa” (black solid line) and 90% **...**

We have shown how to use Stochastic Approximation Monte Carlo (SAMC) for Bayesian principal component selection and model estimation for longitudinal functional data. Our simulations suggest that our method is comparable to the method of Zhou et al. (2008) in terms of frequentist mean squared errors when estimating the mean function, and generally considerably more accurate when estimating the random variation about the mean function. While there are statistical gains in our approach, it is still MCMC based, and hence in terms of computing time generally much slower than the EM-based algorithm. Both methods scale up similarly when the number of subjects increases.

An interesting case to consider would be where the random functions *g _{i}*(·) in (1) or their spline equivalents

Another interesting question is how to deal with correlation when, for example, there are pairs of responses such as might occur for familial data where different family members are followed in time (Sneddon & Sutradhar, 2004). Zhou, et al. (2008) discuss a generalization of model (2), this may be one way of attacking the problem from the principal components framework: they use an EM-algorithm. It would not be a great leap to generalize the SAMC algorithm to handle this case.

We wish to thank the editor and a referee for their detailed reading of the original manuscript and their many suggestions that led to a significant improvement in the paper. We also thank the participants of the International Symposium in Statistics on Inferences in Generalized Linear Longitudinal Mixed Models, whose comments led us to rethink important parts of this work.

Martinez was supported by a post-doctoral training grant from the National Cancer Institute (R25T-CA90301), Carroll's research was supported by a grant from the National Cancer Institute (CA57030), Liang's research was supported by a grant from the National Science Foundation (DMS-0706755) and Zhou was supported by a grant from the National Science Foundation (DMS-0907170). Both Carroll and Liang were also supported by Award Number KUS-CI-016-04, made by King Abdullah University of Science and Technology (KAUST).

We provide some detail of integrations which lead to the marginal posterior density

$$f({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2}\mid \mathit{Y}).$$

Recall that ξ_{1} = log(σ_{ϵ}), ξ_{2} = log(σ_{μ}) and **ξ**_{3} = log{diag(*D*_{α})}. Our prior distributions are in terms of (ξ_{1}, ξ_{2}, **ξ**_{3}), but since they are not being integrated our displayed equations will use (σϵ, σ_{μ}, *D*_{α}).

Our prior distributions were ${\xi}_{1}=\mathrm{Normal}({\mu}_{1},{\sigma}_{1}^{2})$, ${\xi}_{2}=\mathrm{Normal}({\mu}_{2},{\sigma}_{2}^{2})$, and ${\xi}_{3}=\mathrm{Normal}({\mu}_{3},{\sigma}_{j}^{2}I)$. It is more tricky to penalize Θ_{f} as done by Zhou et al. (2008), because Θ_{f} must be orthogonal. We wrote the prior for $[{\sigma}_{f}^{2}\mid {\Omega}_{f}]$ to be $\mathrm{IG}(A+pM\u22152,B+{\sum}_{\ell =1}^{M}{\vartheta}_{j}^{\mathrm{T}}K{\vartheta}_{j})$, with *A* = *B* = 3, where _{j} is the *j*^{th} column of Θ_{f}. This has the convenient behavior that samples from this distribution are of the same form as samples for a smoothing parameter in an ordinary penalized spline regression without constraints, see Carroll, Ruppert, Tosteson, Crainiceanu & Karagas (2004). The prior distribution for Ω_{f} is independent uniforms on [−π/2, π/2], and is written as *p*(Ω_{f}) = π^{−M(p−1)}*I*(−π/2 < Ω_{f} < π/2).

Let the prior density for $({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2})=p({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2})$. Then except for a normalizing constant, the posterior density is

$$f({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2},{\theta}_{\mu},\alpha \mid \mathit{Y})=\prod _{i=1}^{n}\left\{p({Y}_{i}\mid {\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2},{\theta}_{\mu},\alpha )p({\alpha}_{i}\mid {\xi}_{3})\right\}p\left({\theta}_{\mu}\right)p({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2}).$$

For convenience in displaying the equations, where it causes no confusion we will write this posterior in terms of (σ_{ϵ}, σ_{μ}, *D*_{α}), which is given as

$$\begin{array}{cc}\hfill & {\left({\sigma}_{\u03f5}^{2}\right)}^{-N\u22152}\mathrm{exp}\{-{\left(2{\sigma}_{\u03f5}^{2}\right)}^{-1}\sum _{i=1}^{n}{({Y}_{i}-{B}_{i}{\theta}_{\mu}-{B}_{i}{\Theta}_{f}{\alpha}_{i})}^{T}({Y}_{i}-{B}_{i}{\theta}_{\mu}-{B}_{i}{\Theta}_{f}{\alpha}_{i})\}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\times {\mid 2\pi {D}_{\alpha}\mid}^{-n\u22152}\mathrm{exp}\{-(1\u22152)\sum _{i=1}^{n}{\alpha}_{i}^{T}{D}_{\alpha}^{-1}{\alpha}_{i}\}\times {\mid {\sigma}_{\mu}^{2}{K}^{-1}\mid}^{-1\u22152}\mathrm{exp}\{-(1\u22152){\theta}_{\mu}^{T}(K\u2215{\sigma}_{\mu}^{2}){\theta}_{\mu}\}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\times p({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2}),\hfill \end{array}$$

where $N={\sum}_{i=1}^{n}{n}_{i}$. Integrate with respect to α = (α_{1}, …, α_{n})

$$\begin{array}{cc}\hfill f({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2},{\theta}_{\mu}\mid \mathit{Y})=& {\left({\sigma}_{\u03f5}^{2}\right)}^{-N\u22152}\times {\mid 2\pi {D}_{\alpha}\mid}^{-n\u22152}\times {\mid {\sigma}_{\mu}^{2}{K}^{-1}\mid}^{-1\u22152}\times \left(\prod _{i=1}^{n}{\mid {\mathcal{C}}_{2i}\mid}^{1\u22152}\right)\hfill \\ \hfill & \times \mathrm{exp}\{-{\left(2{\sigma}_{\u03f5}^{2}\right)}^{-1}\sum _{i=1}^{n}{R}_{i}^{\mathrm{T}}{R}_{i}-(1\u22152){\theta}_{\mu}^{\mathrm{T}}(K\u2215{\sigma}_{\mu}^{2}){\theta}_{\mu}\phantom{\}}\hfill \\ \hfill & \phantom{\{}+(1\u22152)\sum _{i=1}^{n}{\mathcal{C}}_{1i}^{\mathrm{T}}{\mathcal{C}}_{2i}{\mathcal{C}}_{1i}\}\times p({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2}),\hfill \end{array}$$

where ${\mathcal{C}}_{2i}={({\Theta}_{f}^{\mathrm{T}}{B}_{i}^{\mathrm{T}}{B}_{i}{\Theta}_{f}\u2215{\sigma}_{\u03f5}^{2}+{D}_{\alpha}^{-1})}^{-1}$, ${\mathcal{C}}_{1i}^{\mathrm{T}}=({R}_{i}^{\mathrm{T}}{B}_{i}{\Theta}_{f}\u2215{\sigma}_{\u03f5}^{2})$, and *R _{i}* =

$$\begin{array}{cc}\hfill f({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2}\mid \mathit{Y})=& {\left({\sigma}_{\u03f5}^{2}\right)}^{-N\u22152}{\mid 2\pi {D}_{\alpha}\mid}^{-n\u22152}{\mid {\sigma}_{\mu}^{2}{K}^{-1}\mid}^{-1\u22152}\left(\prod _{i=1}^{n}{\mid {\mathcal{C}}_{2i}\mid}^{1\u22152}\right)\times {\mid {\mathcal{C}}_{4}\mid}^{1\u22152}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\mathrm{exp}\{-(1\u22152)\sum _{i=1}^{n}{Y}_{i}^{\mathrm{T}}({\sigma}_{\u03f5}^{-2}{\mathrm{I}}_{{n}_{i}}-{\Psi}_{i}){Y}_{i}+(1\u22152){\mathcal{C}}_{3}^{\mathrm{T}}{\mathcal{C}}_{4}{\mathcal{C}}_{3}\}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\times p({\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2}),\hfill \end{array}$$

(6)

where ${\mathcal{C}}_{4}={({\sigma}_{\mu}^{-2}K+{\sigma}_{\u03f5}^{-2}{\sum}_{i=1}^{n}{B}_{i}^{\mathrm{T}}{B}_{i}-{\sum}_{i=1}^{n}{B}_{i}^{\mathrm{T}}{\Psi}_{i}{B}_{i})}^{-1}$, ${\mathcal{C}}_{3}^{\mathrm{T}}=({\sigma}_{\u03f5}^{-2}{\sum}_{i=1}^{n}{Y}_{i}^{\mathrm{T}}{B}_{i}-{\sum}_{i=1}^{n}{Y}_{i}^{\mathrm{T}}{\Psi}_{i}{B}_{i})$ and ${\Psi}_{i}={\sigma}_{\u03f5}^{-4}{B}_{i}{\Theta}_{f}{\mathcal{C}}_{2i}{\Theta}_{f}^{\mathrm{T}}{B}_{i}^{\mathrm{T}}$.

Easy calculations show that samples from θ_{μ} and α_{i} given the data and the parameters can be obtained as

$$\begin{array}{cc}\hfill [{\theta}_{\mu}\mid {\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2},\mathit{Y}]& =\mathrm{Normal}({\mathcal{C}}_{4}{\mathcal{C}}_{3},{\mathcal{C}}_{4});\hfill \\ \hfill [{\alpha}_{i}\mid {\Omega}_{f},{\xi}_{3},{\xi}_{2},{\xi}_{1},{\sigma}_{f}^{2},{\theta}_{\mu},\mathit{Y}]& =\mathrm{Normal}({\mathcal{C}}_{2i}{\mathcal{C}}_{1i},{\mathcal{C}}_{2i}).\hfill \end{array}$$

- Baladandayuthapani V, Hong MY, Mallick BK, Lupton JR, Turner ND, Carroll RJ. Bayesian hierarchical spatially correlated functional data analysis with application to colon carcino-genesis. Biometrics. 2008;64:64–73. [PMC free article] [PubMed]
- Besse P, Ramsay JO. Principal components-analysis of sampled functions. Psychometrika. 1986;51:285–311.
- Boente G, Fraiman R. Kernel-based functional principal components. Statistics & Probability Letters. 2000;48:335–345.
- Cardot H. Nonparametric estimation of smoothed principal components analysis of sampled noisy functions. Journal of Nonparametric Statistics. 2000;12:503–538.
- Carroll RJ, Ruppert D, Tosteson TD, Crainiceanu C, Karagas MR. Nonparametric regression and instrumental variables. Journal of the American Statistical Association. 2004;99:736–750.
- Green PJ. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika. 1995;82:711–732.
- Hall P, Hosseini-Nasab M. On properties of functional principal components analysis. Journal of the Royal Statistical Society, Series B. 2006;68:109–126.
- Hall P, Müller HG, Wang JL. Properties of principal component methods for functional and longitudinal data analysis. Annals of Statistics. 2006;34:1493–1517.
- Huang JZ, Shen H, Buja A. Functional principal components analysis via penalized rank one approximation. Electronic Journal of Statistics. 2008;2:678–695.
- Hong MY, Bancroft LK, Turner ND, Davidson LA, Murphy ME, Carroll RJ, Chapkin RS, Lupton JR. Fish oil decreases oxidative DNA damage by enhancing apoptosis in rat colon. Nutrition and Cancer. 2005;52(2):166–175. [PubMed]
- James GM, Hastie TJ, Sugar CA. Principal component models for sparse functional data. Biometrika. 2000;87:587–602.
- Jank W, Shmueli G. Functional data analysis in electronic commerce research. Statistical Science. 2006;21:155–166.
- Lederman MM, Connick E, Landay A, Kuritzkes DR, Spritzler J, Clair M, Kotzin BL, Fox L, Chiozzi MH, Leonard JM, Rousseau F, Wade M, Roe JD, Martinez A, Kessler H. Immunological responses associated with 12 weeks of combination antiretroviral therapy consisting of Zidovudine, Lamivudine & Ritonavir: Results of AIDS Clinical Trials Group Protocal 315. J. Infectious Diseases. 1998;178:70–79. [PubMed]
- Li Y, Wang N, Hong MY, Turner ND, Lupton JR, Carroll RJ. Nonparametric estimation of correlation functions in longitudinal and spatial data, with application to colon carcinogenesis experiments. Annals of Statistics. 2007;35:1608–1643.
- Liang F. On the use of stochastic approximation Monte Carlo for Monte Carlo integration. Statistics and Probability Letters. 2009;79:581–587.
- Liang F, Liu C, Carroll RJ. Stochastic approximation in Monte Carlo computation. Journal of the American Statistical Association. 2007;102:305–320.
- Ocaña FA, Aguilera AM, Escabias M. Computational considerations in functional principal component analysis. Computational Statistics. 2007;22:449–465.
- Sneddon G, Sutradhar BC. On semiparametric familial longitudinal models. Statistics & Probability Letters. 2004;69(3):369–379.
- Ramsay JO, Dalzell CJ. Some tools for functional data analysis. Journal of the Royal Statistical Society, Series B. 1991;53:539–572. With discussion.
- Ramsay JO, Silverman BW. Functional Data Analysis. 2nd Edition Springer; New York: 2005.
- Reiss PT, Ogden T. Functional principal component regression and functional partial least squares. Journal of the American Statistical Association. 2007;102:984–996.
- Rice JA, Silverman BW. Estimating the mean and covariance structure nonparametrically when the data are curves. Journal of the Royal Statistical Society, Series B. 1991;53:233–243.
- Rice JA, Wu CO. Nonparametric mixed effects models for unequally sampled noisy curves. Biometrics. 2001;57:253–59. [PubMed]
- Rice JA. Functional and longitudinal data analysis: Perspectives on smoothing. Statistical Sinica. 2004;14:613–29.
- Silverman BW. Smoothed functional principal components analysis by choice of norm. Annals of Statistics. 1996;24:1–24.
- Wu H, Ding A. Population HIV-1 dynamics in vivo: application models and inference tools for virological data from AIDS clinical trials. Biometrics. 1999;55:410–18. [PubMed]
- Yao F, Müller HG, Wang JL. Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association. 2005a;100:577–90.
- Yao F, Müller HG, Wang JL. Functional linear regression analysis for longitudinal data. Annals of Statistics. 2005b;33:2873–903.
- Zhang D, Davidian M. Linear mixed models with flexible distributions of random effects for longitudinal data. Biometrics. 2001;57:795–802. [PubMed]
- Zhou L, Huang JZ, Carroll RJ. Joint modeling of paired sparse functional data using principal components. Biometrika. 2008;95:601–619. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |