International Journal of Clinical Biostatistics and Biometrics
Dynamic Factor Analysis for Multivariate Time Series: An Application to Cognitive Trajectories
Yorghos Tripodis^{1*} and Nikolaos Zirogiannis^{2}
^{1}Department of Biostatistics, Boston University, USA
^{2}School of Public and Environmental Affairs, Indiana University, USA
^{*}Corresponding author:
Yorghos Tripodis, Assistant Professor of Biostatistics, Boston University School of Public Health, USA, Email: yorghos@bu.edu
Int J Clin Biostat Biom, IJCBB1001, (Volume 1, Issue 1), Research Article; ISSN: 24695831
Received: July 23, 2015  Accepted: August 26, 2015  Published: August 28, 2015
Citation: Tripodis Y, Zirogiannis N (2015) Dynamic Factor Analysis for Multivariate Time Series: An Application to Cognitive Trajectories. Int J Clin Biostat Biom 1:001. 10.23937/24695831/1510001
Copyright: © 2015 Tripodis Y, et al. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Abstract
We propose a dynamic factor model appropriate for large epidemiological studies and develop an estimation algorithm which can handle datasets with large number of subjects and short temporal information. The algorithm uses a two cycle iterative approach for parameter estimation in such a large dataset. Each iteration consists of two distinct cycles, both following an EM algorithm approach. This iterative process will continue until convergence is achieved. We utilized a dataset from the National Alzheimer Coordinating Center (NACC) to estimate underlying measures of cognition based on a battery of observed neuropsychological tests. We assess the goodness of fit and the precision of the dynamic factor model estimators and compare it with a nondynamic version in which temporal information is not used. The dynamic factor model is superior to a nondynamic version with respect to fit statistics shown in simulation experiments. Moreover, it has increased power to detect differences in the rate of decline for a given sample size.
Keywords
Dynamic factor models, EM algorithm, Panel data, Statespace models, Cognition, Alzheimer's disease, Neuropsychological performance
Introduction
An increasing number of studies provide data on many variables across a large number of individuals in a longitudinal setting. These studies give us an unprecedented opportunity for epidemiological research on psychological measures over time and between subpopulations. For example, many observational and clinical trial studies of cognitive aging use neuropsychological test batteries to assess overall cognition and its specific domains [1]. Statistical tools have been developed to extract information from these evaluation tests in order to estimate a single or multiple latent cognitive indices. A common method for the estimation of such latent variables is confirmatory factor analysis (CFM) [2,3]. The repeated nature of the studies is often ignored in these models, even though recent studies attempt to capture the temporal information in order to increase performance of measures for cognitive change [46]. In psychology, the problem of estimating latent structures from repeated measurements has been addressed by the use of dynamic factor models (DFM). There is a long literature on estimation theory and applications of DFMs, either direct autoregressive factor models (DAFS) or white noise factor score models (WNFS) [710]. These methods have been applied in studies with a relatively long time series (at least 70 repeated observations per subject) of a single or a small number of individuals [1113]. On the other hand, in longitudinal epidemiological studies, we often have a large number of subjects with a very short time series, usually nonstationary, with typically 2 to 5 repeated observations per subject. Furthermore, in order to compare latent variables across subpopulation of interest crosssectionally and over time, we require that the estimated factor scores are comparable across individuals. This implies certain restrictions on the factor model. Molenaar et al. (1992) extended DFMs to nonstationary time series [14]. Although models for large number of individuals and short time series are theoretically feasible by applying standard multivariate time series methods, they are computationally restrictive. As the number of individuals becomes large, so does the number of parameters to be estimated, and direct optimisation becomes harder and more time consuming. Markov Chain Monte Carlo methods for parameter estimation have been implemented [15,16], but these methods require long time series, typically in the hundreds of observations over time. On the other hand, Zuur et al. (2003) showed that parameter estimation using the EM (ExpectationMaximazation) algorithm in a relatively short time window (30 repeated observations) for 12 crosssectional units can be easily implemented [17]. This methodological approach has its own limitations since applying the EM algorithm to a larger number of individuals and even shorter time series would substantially increase computation time [18]. The limitations of the data and complications of the estimation methods have resulted in a very narrow interest for DFM applications in large epidemiological studies.
In this paper, we provide statistical tools for analyzing DFMs using data which are typical in epidemiological studies with large number of participants and short nonstationary time series. Specifically, we develop an estimation algorithm, extending the classic EM algorithm, by developing an iterative twocycle estimation process, following the steps of the ECME (Expectation/Conditional Maximization Either) algorithm [19]. This estimation method is flexible enough to be applicable in studies with multiple individuals, and short unequally spaced temporal information.
We apply the dynamic factor model to a variety of neuropsychological tests using data from the National Alzheimers Coordinating Center (NACC) study and estimate a smooth cognitive measure for each individual's total cognition as well as measures for specific cognitive domains, such as memory, attention and language. We hypothesize that by incorporating longitudinal information into the factor models we increase the accuracy of the estimates of change over time and consequently increase power to detect differences between groups. We focus on a casecontrol sample using data where participants are selected to be cognitively normal. Cases in this study include participants that will convert to Mild Cognitive Impairment (MCI) after the period used in the analysis, while controls will continue having normal cognition for the next two followup visits after the end of the analytic period.
In the next section we describe the dynamic factor model and its estimation method. In §3, we assess its performance in estimating the underlying cognition and its domains and compare it with a factor model which ignores any longitudinal information (nondynamic model). We apply the dynamic and nondynamic versions of the factor model to clinical data collected by NACC, and compare their power of detecting differences in the rate of cognitive change for various sample sizes. Finally in §4, we conclude with a discussion of the methods and results including limitations and directions for future studies.
Methods
In this section we describe the dynamic version of the factor model and its estimation process. The difference of the dynamic factor model for panel data from the nondynamic version is that the former captures not only correlations between input variables but also autocorrelations and cross correlations of these variables of interest.
Model
We let U_{t} denote the nq×1 vector containing the unobserved cognitive indices of factors for subjects (with q < n) at time t, with t = 1,..., T. We assume that the dynamic properties of U_{t} can be captured by a Markov process. For illustration purposes, and without loss of generality, we first present the case where we have equally spaced observations and equal number of neuropsychological tests for each subject. We subsequently present the model for the general case with unequally spaced or missing observations. Hence, we form the following linear Gaussian state space model:
$${y}_{t}=B{U}_{t}+{e}_{t},\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.05em}}\text{\hspace{1em}}{e}_{t}\sim N(0,D),\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.1)$$ $${U}_{t+1}=T{U}_{t}+{\eta}_{t},\text{\hspace{1em}}\text{\hspace{1em}}{\eta}_{t}\sim N(0,{I}_{n}),\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.2)$$where B is the matrix of factor loadings with dimensions $np\times nq$ , with denoting the number of observed variables, ${y}_{t}$ is a $np\times 1$ vector of observed neuropsychological measures per individual, T is $np\times nq$ transition matrix and ${I}_{n}$ is a $nq\times nq$ identity matrix and ${e}_{t}$ and ${\eta}_{t}$ are error terms [20,21]. The state space formulation described in (2.1) and (2.2), models the behavior of the unobserved state vector U_{t} over time using the observed values y1,....., yn . The state vector U_{t} is assumed to be independent of the error terms e_{t} and ${\eta}_{t}$ for all t. In addition, the error terms e_{t} and ${\eta}_{t}$ are assumed to be independent, identically distributed (i.i.d.) [22,23]. In general, the model defined by equations (2.1) and (2.2) is not identifiable. Zirogiannis and Tripodis (2014) state the conditions for identifiability for a general dynamic factor model [24]. In order for the model in (2.1) and (2.2) to be identifiable we must impose a certain structure. We first assume that the unobserved cognitive indices follow a multivariate random walk, so that $T={I}_{n}$. This is a reasonable assumption when modeling cognition for an aging population where the spacing of the observation period is roughly annual. Similar nonstationary models for psychological constructs have been suggested by Molenaar and Campbell (2009) and used, among others, by Hekler et al. (2013) and Gu et al. (2014) [12,25,26] . We also impose a structure on the factor loading matrix B and the variance of the idiosyncratic errors D. We assume that the factor loadings for each observed variable are the same for each individual in the study. This assumption is necessary in order to have comparable estimated cognitive indices across individuals. We also assume that participants in the study are conditionally independent and that the variance of the idiosyncratic errors is the same for all individuals. These assumptions result in a block diagonal structure for D. The imposed structure results in a model that is fully identifiable. The model can be easily extended for cases where nondiagonal block elements of B and D are not zero. The choice of the structure of B and D is specific to our application of interest.
Unequally spaced and missing observations
It is very common in longitudinal observational studies to have unequally spaced or missing observations. Let ${\tau}_{it}$ be the distance between observations t and $t+{\tau}_{it}$ of the i^{th} subject, and ${\tau}_{t}$ the vector with the distances between two subsequent observations at time. Then we can rewrite the statespace form of the multivariate random walks as:
$${y}_{t}=B{U}_{t}+{e}_{t},\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.05em}}\text{\hspace{1em}}{e}_{it}\sim N(0,D)$$ $${U}_{t+\tau}={U}_{t}+{\eta}_{t},\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}{\eta}_{t}\sim N(0,{Q}_{t})$$where ${Q}_{t}={\displaystyle \sum _{i=1}^{q}{E}_{i}{\tau}_{t}{e}_{i}}$ with E_{i} a $q\times q$ matrix with 1 for the element (i, i) and 0 everywhere else and ${e}_{i}$ is a $1\times q$ vector with 1 for the element i and 0 everywhere else. This timevarying model can be used for unequally spaced and missing observations, as well as for forecasting for any ${\tau}_{n}$ steps ahead.
2step modified ECME Algorithm
The high dimensionality of the data vector y_{t} makes estimation of our model rather problematic. Moreover, in biomedical applications such as the one we explore in this paper, we deal with cases where T is very small while n is very large. Usual Newtontype gradient methods do not work in this situation creating the need for a novel estimation approach. We introduce a modified ECME algorithm that makes estimation of the model specified in (2.1) and (2.2), feasible through an iterative twocycle process. The 2cycle modified ECME algorithm is an extension of the ECME algorithm developed by Liu and Rubin (1998), which itself is an extension of the widely known EM algorithm [27]. The modified ECME algorithm starts by partitioning the vector of unknown parameters $\Psi $ into (${\Psi}_{1}$ ,${\Psi}_{2}$ ) where ${\Psi}_{1}$ contains the elements of D that need to be estimated, while ${\Psi}_{2}$ contains the relevant elements of B. We use the term "cycle" as an intermediary between a "step" and an "iteration" as in Meng and Dyk (1997) [28]. In the case of our modified ECME algorithm, every iteration is comprised of two cycles. Each cycle includes one Estep and one Mstep, where the first cycle estimates ${\Psi}_{1}$ and ${\Psi}_{2}$ given the estimates of $\Psi $ of the previous iteration, while the second cycle estimates ${\Psi}_{2}$ given the estimates of $\Psi $ of the previous cycle.
The functional form of the completedata loglikelihood at time period t is [29]:
$$\mathrm{log}\mathcal{l}\left(\Psi \right)=\frac{1}{2}\mathrm{log}\left\{\left{D}^{1}\right\right\}\frac{1}{2}\left\{{\left({y}_{t}B{u}_{t}\right)}^{\text{'}}{D}^{1}\left({y}_{t}B{u}_{t}\right){\left({u}_{t+1}{u}_{t}\right)}^{\text{'}}\left({u}_{t+1}{u}_{t}\right)\right\}$$
Since ${u}_{t}$ is unobserved, we can consider it missing and use the EM algorithm framework. In order to find the MLE, we need to calculate the distribution of the latent variable u_{t} conditional on the observed values of ${y}_{t}$. There is a long literature describing the EM procedure for factor analysis in crosssectional data starting with Rubin and Thayer (1982) [30]. Applying the EM framework for longitudinal data we need to condition not only on the concurrent observed value of ${y}_{t}$ but on all the previous observed history y_{1}, . . . , ${y}_{t}$. As we will see in the following two subsections, we use the first cycle to obtain estimates for ${u}_{t}$ by conditioning on the concurrent observed variables, ${y}_{t}$, and the second to update these estimates by conditioning on the history of the observed variable, ${y}_{1}$, . . . , ${y}_{t}$ using the Kalman filter [31]. This iterative process will continue until the likelihood function stops increasing and convergence is achieved.
First cycle: During the k^{th} iteration of the first cycle, the Estep of the 2cycle ECME algorithm is:
$${Z}_{\Psi}({\Psi}_{1},{\Psi}_{2};{\Psi}_{1}{}^{(k1)},{\Psi}_{2}{}^{(k1)})={E}_{\Psi}\left\{{\displaystyle \sum _{t=1}^{T}\mathrm{log}{\mathcal{l}}_{t}\left[\left({\Psi}_{1},{\Psi}_{2}\right){y}_{t},{\Psi}_{1}^{(k1)},{\Psi}_{2}^{(k1)}\right]}\right\}.\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.3)$$
Following the notation presented in [29], the sufficient statistics are calculated in the (k1) iteration by the following equations:
$${\gamma}^{(k1)}={\left({B}^{\left(k1\right)}{B}^{\left(k1\right)}+{D}^{(k1)}\right)}^{1}{B}^{(k1)}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.4)$$ $${\omega}^{(k1)}=I{\gamma}^{{(k1)}^{\text{'}}}{B}^{(k1)}$$
The first Mstep involves differentiating ${Z}_{\Psi}({\Psi}_{1},{\Psi}_{2};{\Psi}_{1}{}^{(k1)},{\Psi}_{2}{}^{(k1)})$ with respect to ${\Psi}_{1}$ and ${\Psi}_{2}$ in order to obtain ${\Psi}_{1}^{(k)}$ and ${\Psi}_{2}^{(k/2)}$ :
$${Z}_{\Psi}({\Psi}_{1}^{(k)},{\Psi}_{2}^{(k/2)};{\Psi}_{1}{}^{(k1)},{\Psi}_{2}{}^{(k1)})\ge {Z}_{\Psi}({\Psi}_{1},{\Psi}_{2};{\Psi}_{1}{}^{(k1)},{\Psi}_{2}{}^{(k1)}),\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.5)$$
The firstcycle Mstep is identical to the Mstep of the traditional EM algorithm for factor analysis models [32]:
$${B}^{(k/2)}={C}_{yy}{\gamma}^{(k1)\text{'}}{\left\{{\gamma}^{(k1)\text{'}}{C}_{yy}{\gamma}^{(k1)}+n{\omega}^{(k1)}\right\}}^{1},\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.6)$$ $${D}^{(k)}={n}^{1}diag\left\{{C}_{yy}{C}_{yy}{\gamma}^{(k1)}{B}^{\text{'}}\right\},\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.7)$$
Where ${C}_{yy}$ is the sample unconditional covariance matrix of ${Y}_{T}={({y}_{1}^{\text{'}},\mathrm{...},{y}_{T}^{\text{'}})}^{\text{'}},$ i.e. $E({Y}_{T}{Y}_{T}^{\text{'}})={C}_{yy}$ . At the end of the first cycle we have updated estimates for all the elements of the variance matrix of the idiosyncratic errors, D, and intermediate estimates for the matrix of factor loadings, B. We use these estimates in the second cycle to get updated estimates for the factor loadings.
Second cycle: In the Estep of the second cycle we estimate ${\Psi}_{2}^{(k)}$. We proceed by calculating:
$${Z}_{{\Psi}_{2}}({\Psi}_{2};{\Psi}_{1}^{(k)},{\Psi}_{2}^{(k/2)})={E}_{{\Psi}_{2}}\left\{{\displaystyle \sum _{t=1}^{T}\mathcal{l}[{\Psi}_{2}{Y}_{t1},{\Psi}_{1}^{(k1)},{\Psi}_{2}^{(k/2)}]}\right\}\text{\hspace{1em}}\text{\hspace{1em}}(2.8)$$
The second Estep involves forming the expected completedata log likelihood conditional on ${Y}_{t1},$ which is the set of past observations ${y}_{1,\mathrm{...},}{y}_{t1}.$ The subsequent Mstep involves differentiating ${Z}_{{\Psi}_{2}}({\Psi}_{2};{\Psi}_{1}^{(k)},{\Psi}_{2}^{(k/2)})$ with respect to ${\Psi}_{2}$. We choose ${\Psi}_{2}^{(k)}$ such that:
$${Z}_{{\Psi}_{2}}({\Psi}_{2};{\Psi}_{1}^{(k)},{\Psi}_{2}^{(k/2)})\ge {Z}_{{\Psi}_{2}}({\Psi}_{2};{\Psi}_{1}^{(k)},{\Psi}_{2}^{(k/2)}).\text{\hspace{0.05em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.9)$$
Upon maximization of ${Z}_{{\Psi}_{2}}$ , the estimate ${\Psi}_{2}^{(k)}$ is used in the Estep of the first cycle of the next iteration. We calculate and maximize ${Z}_{{\Psi}_{2}}({\Psi}_{2};{\Psi}_{1}^{(k)},{\Psi}_{2}^{(k/2)})$ by using the prediction error decomposition of the conditional likelihood [33]:
$$\mathrm{log}{\mathcal{l}}_{t}\left({\Psi}_{2}\right)=\mathrm{log}\frac{1}{2\pi}\frac{1}{2}[\mathrm{log}{F}_{t}+\text{}{\upsilon}_{t}^{\text{'}}{F}_{t}^{1}\upsilon ]\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(2.10)$$
Where ${\upsilon}_{t}$ is the prediction error conditional on past history and $F$ is its variance. Quantities, ${\upsilon}_{t}$ and $F$ can be estimated with the use of the Kalman filter, which is a set of recursions which allow information about the system to be updated every time an additional observation ${Y}_{t}$ is introduced [21]. Once ${\upsilon}_{t}$ and $F$ are calculated, (2.10) is maximized with respect to ${\Psi}_{2}$, as illustrated in (2.9).
Results
In the next section, we assess and apply the model and the estimation process described in §2. We first assess the performance of the 2cycle ECME estimator using a simulation study. We then apply the model in data from the NACC study. We also compare the dynamic factor model with a nondynamic version in which temporal information is not used. The nondynamic version is observationally equivalent to a Confirmatory Factor Analysis (CFM) model and it is defined solely by the observation equation (2.1). Estimation was done using a code written by the authors in ${O}_{x}$ Programming Language [34].
Simulation
The model from which we simulate is a variant used by Doz et al. (2011) which is based on a simulation scheme used by Stock and Watson (2002) [35,36]. We define:
$$B=\left(\begin{array}{cccc}f& 0& \cdots & 0\\ 0& f& \cdots & 0\\ \vdots & \cdots & \ddots & \vdots \\ 0& \cdots & \cdots & f\end{array}\right),\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}D=\left(\begin{array}{cccc}d& 0& \cdots & 0\\ 0& d& \cdots & 0\\ \vdots & \cdots & \ddots & \vdots \\ 0& \cdots & \cdots & d\end{array}\right),\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}(3.1)$$
With
1. f a $p\times 1$ vector of factors loadings with f_{[K]} ∼ U(0,1) subject to $\sum _{k=1}^{p}{f}_{[K]}=\text{}{1}_{,}$
2. d a $p\times p$ diagonal matrix of variances for the idiosyncratic elements, with ${d}_{\left[k\right][k]}={f}_{[k]}\frac{{\beta}_{k}}{1{\beta}_{k}}$ with ${\beta}_{k}\text{}~u\text{(0}\text{.1,0:9)}$
where $k=1,\mathrm{...},p.$ We generate 1000 replicates from the model defined by (2.1), (2.2) and (3.1) with ${\text{U}}_{0}\text{}~{\text{N(0,I}}_{n}\text{),}$ for different combinations of sizes for observed tests, p, number of subjects n, and time points, T. Specifically, we use p = 5,10,15, n = 10,50,100,200,300 and T = 3,5,7,10,15.
The choice of these values corresponds to our specific application. We specify factor loadings which are the same across individuals who do not share any familial or other relationship. The coefficient ${\beta}_{k}$ is the ratio between the variance of the idiosyncratic component, e_{t} , and the total variance of the corresponding observed variable, Y_{t}. In the simulation, this ratio is drawn from a uniform distribution between the interval of (0.1, 0.9). This interval was chosen in order to avoid parameters at the boundary of the parameter space.
Estimation was done using the 2cycle modified ECME and yielded estimates of the factor $\widehat{U}={({\widehat{U}}_{1},{\widehat{U}}_{2},\mathrm{....}{\widehat{U}}_{T})}^{\text{'}}$. Performance was measured by the trace statistic: $$\frac{tr({\text{U}}^{\text{'}}\hat{U}{(\hat{U}\text{'}\hat{U})}^{1}\hat{U}\text{'U})}{tr({\text{U}}^{\text{'}}\text{U})}.$$
The trace statistic is a multivariate version of the R^{2} of the regression of the true factors on the estimated factors [35]. A number close to 1 implies a good approximation of the estimated latent variable to the true factor. We used the trace statistic, TR_{DFM}, as a performance measure for the Dynamic factor model. We also obtained estimates of the latent factor using a nondynamic factor model defined only by equation (2.1). In order to have comparable results, we also used the 2cycle modified ECME for the nondynamic version of the factor model. We then calculated the equivalent trace statistic, TR_{CFM} for the nondynamic model. We use the ratio of the two trace statistics, $\frac{T{R}_{DFM}}{T{R}_{CFM}}$, as a comparison measure of the two models. Values above 1 imply that the dynamic model has superior performance to the nondynamic version, with respect to how close the estimated and the true factors are.
Table 1: Performance of factor estimators from 1000 simulations. Trace statistics for various time lengths and sample sizes.
View Table 1
Table 1 reports the results of the trace statistics from the simulation experiment. The numbers in the table refer to the average across 1000 replicates. As an example, we use the case of T = 3, n = 300, and p = 5; 86% of the variability of the true, simulated factor is explained by the factor estimated by the DFM. The explained variability using DFM is 1% higher than the explained variability using CFM. The goodness of fit of the estimated factors, as measured with TR_{DFM}, increases with the size of individuals n in the sample, and the number of repeated observations per individual T. TR_{DFM} varies from 0.77 for a smalln, and smallT sample to 0.97 for a moderate nmoderate T. The goodness of fit of the estimators does not improve as the number of observed tests p per individual increases for a given size n and T. Moreover, the dynamic factor model always performs at least as good as the nondynamic version. The relative performance of the dynamic factor model increases with T. Based on $\frac{T{R}_{DFM}}{T{R}_{CFM}}$ , the relative performance of the dynamic model ranges from 1% better goodness of fit compared to the nondynamic version when n = 10, p = 5 and T = 3 to 9% when n = 300, p = 15 and T = 15.
Application
Alzheimers Disease (AD), the most common form of dementia, is a significant cause of disability and mortality among the elderly. The latest figures show that 5.2 million people in the US, approximately 14% of the population over age 70, are afflicted by AD [37]. As the population ages over the next several decades, this number is expected to increase [38]. The only definitive way to diagnose AD is postmortem, but neuropsychiatrists reach a premortem diagnosis by reviewing and discussing the subject's clinical history, as well as scores from a variety of neuropsychological evaluation tests [39]. The results of the neuropsychological tests which are part of the batteries can exhibit high withinsubject variability [40] and may make diagnosis difficult. Moreover, the emphasis in Alzheimers disease clinical research has shifted to developing interventions before symptoms onset. In order to address this need, researchers are required to develop cognitive measures which discriminate between cognitively healthy subjects and individuals with small cognitive changes who will convert to mild cognitive impairment (MCI).
We used the NACC dataset with visits from September 2005 to June 2013 for testing and evaluation. NACC serves as a repository for data collected at 34 past and present Alzheimer's Disease Centers (ADCs) throughout the United States. The ADCs conduct clinical and biomedical research on Alzheimer's disease and related disorders. Centers enroll their study subjects in various ways, including referral from clinicians, selfreferral by patients themselves or concerned family members, active recruitment through community organizations, and volunteers who wish to contribute to research. Most centers also enroll volunteer control subjects. Study subjects at each center are best regarded as a case series, not necessarily representing all cases of disease in a defined population. More information on the study can be found in Morris et al. (2006) [41].
We focus on a study subsample which includes cognitively healthy participants at initial visit. For all subjects we only considered their neuropsychological test results while cognitively healthy, even though some converted to MCI state at a later visit. For those participants who did not convert to MCI during our observation period, we only considered those with at least 4 visits. To avoid the risk of healthy participants converting to MCI at a future visit, we restricted our analytic period by excluding from further analyses the last two measurement occasions for that group. Consequently, normal controls remain cognitively normal for at least 2 years after the end of the analytic period used in this study. For those participants who converted to MCI we considered those with at least 1 followup with normal cognition. We also excluded nonEnglish speakers as well as subjects with a number of comorbidities: history of stroke, history of transient ischemic attack, history of other cerebrovascular diseases, Parkinson's disease or other Parkinsonism disorder, history of seizures, history of any brain trauma or other neurologic conditions, history of depression or psychiatric disorder. We then created two balanced groups, with n = 149 each, matched by age, sex and education which differ only in their future cognitive state: one group converts to MCI at the following visit after the end of the analytic period (converters), while the other group remains cognitively normal for at least the next two subsequent visits (nonconverters) beyond the end of the restricted observational period. The description of the sample is given in figure 1. The mean (SD) age at initial visit is 75.7 (7.5) with 15.4 (2.5) average years of education. There are 170 (57.1%) women in the sample with 3.0 (1.2) visits on average, and 2.3 (1.3) years of followup since the initial visit.
.
Figure 1: Description of analytic sample
View Figure 1
We considered four factor models using different neuropsychological measures according to their relation to a specific domain: i) memory, ii) attentionpsychomotor speed, iii) language and iv) general cognition. For each factor model, we run both a dynamic and a nondynamic version. We estimated all models using a MacBook Pro with a Intel Core i7 2.3GHz processor on OS X Yosemite. Estimation times vary from 12.01 second for the language domain, to 46.97 seconds for the attentionpsychomotor speed domain. For the general cognition factor, which includes all 11 neuropsychological scores, estimation time was 25.31 minutes. For both versions of the model, the one stepahead prediction errors were tested for normality and residual autocorrelation. Even though both the dynamic and the nondynamic version of all four factor models indicated nonnormal errors (e.g. for general cognition, pvalue<0.0001 for BowmanShenton test for normality and for BoxLjung portmanteau test for autocorrelation at lag 1), further investigation showed that this is caused by outliers from seven participants. These participants have significantly lower estimated factors at the last visit, which may indicate misdiagnosis or untimely diagnosis of MCI. For each neuropsychological test, we run mixed effects regressions using PROC MIXED in SAS 9.3 with random intercepts and random slopes for time to test the hypothesis that there are significant differences in the rate of change by group (converters Vs nonconverters). We also used mixed effects regression on the factor scores estimated by the dynamic as well as the nondynamic factor models. Table 2 shows the factor loadings for each domain and table 3 shows the estimated annual rate of change for each of the neuropsychological tests and for the nondynamic and dynamic factors. For ease of comparison, all outcomes have been standardized, using the mean and standard deviation of all cognitively healthy NACC participants. We note that there is no significant annual change for the group of nonconverters for all neuropsychological measures, with the exception of logical memory: delayed, and for the estimated factors. For the converters, only MMSE, Trails B and Verbal Fluency Test: vegetables show significant decrease at the 5% level, while the factors from the simple (nondynamic) factor model for attention and language show significant decrease over time. For the dynamic factor model estimates, all three domains and total cognition estimates show significant decreases over time for the group which progressed to MCI at the next followup period. For the nondynamic factor model estimates however, the total cognition factor and the memory factor do not show any difference. Given that an important feature that leads to an MCI diagnosis is manifestation of significant cognitive decline, it is important to note that the dynamic factor model estimates show evidence of decline even before conversion to MCI. We also note that both the dynamic and the nondynamic version of factor models show significant differences in the annual rate of change between groups. In general, the estimates of difference of the Factor models are larger and have lower pvalues than the estimates of input variables. Furthermore, the estimates of difference from the dynamic factor model are at least as high with larger pvalues as the equivalent estimates of the simple factor models. This difference is due to the fact that DFM incorporates the longitudinal aspect of the psychometric results of every patient. This may be an indication of increased power for the dynamic factor model, which we explore in the next subsection.
Table 2: Factor loadings for each domain
View Table 2
Table 3: Parameter estimates for annual rate of change for all neuropsychological tests
View Table 3
Power analysis: We also investigated the performance of the observed indicators and the estimated factors with respect to power. Our main aim remains to detect differences by group in the rate of change. In order to assess the power of each outcome, we follow a bootstrapping scheme using the data described in the previous section. We first assume that there is, indeed, a difference in the annual rate of change between normal controls and MCI while they are both cognitively normal. For a given sample size n, we perform the following steps:
Simulation scheme
1. Select n/2 matched pairs with replacement.
2. Estimate factor scores for all domains and for total cognition using simple and dynamic factor models.
3. Run a mixed effect regression on the estimated factor using time since first visit, group (converters Vs nonconverters) and time×group interaction, along with age at initial visit as covariates.
4. Is the estimate of time×group interaction significant at the 5% error level?
5. Repeat for 1000 times.
Table 4 shows the power of detecting significant differences for different sample sizes for all outcomes. We note that power of the dynamic factor model estimates is higher than the power of the nondynamic version. For the total summary index, the power for the dynamic version varies from 49.9% for n = 120 to 96.7% for n = 240. The power for the nondynamic version is much lower and goes up to 76.9% for n = 240. These results indicate that a smaller sample size is required for a given power in order to find significant differences in the rate of change by groups. Using the results from table 4, we can calculate the required sample size for both the dynamic and the Nondynamic factor models for an 80% power at α = 5%. For the total cognition index, the DFM model requires a sample size of 187 while the nondynamic version (CFM) requires a sample size of 252. We get similar results for the other domains: memory (^{n}DFM = 370, nCFM = 485), attention (nDFM = 334, nCFM = 546), language (nDF M = 414, nCFM = 1419).
Table 4: Power analysis for group differences
View Table 4
We also note that the power of the factor models is always higher than the power of the individual neuropsychological tests. This indicates that using factor models increases the power of detecting significant differences in the rate of change. One notable exception is the Boston Naming Test (BNT) in the language domain. BNT has a higher power for all sample sizes considered compared to the nondynamic factor estimates. It also has a higher power than the dynamic factor model when n = 120. For larger sample sizes, the dynamic factor model estimates have higher power than BNT.
Conclusion
In this article, we developed an algorithm to estimate a dynamic factor model for data typical in large epidemiological studies and apply it on latent cognitive variables. We compared it with equivalent factor models which do not use temporal information in the estimation, and showed that the dynamic factor model estimates are more accurate as reflected by comparison of fit statistics in simulation experiments. They are also more precise than the nondynamic version estimates as shown by improved power to detect differences in the rate of decline. Since the estimated latent index is a weighted average of the concurrent observed values, the reason for the improved performance of the dynamic factor model is due to the fact that the weighting scheme of DFM takes into account any withinsubject variability over time and any crosscorrelation of tests. In the nondynamic version, weights depend on the correlation between tests as well as on betweensubject variability. Measures that are highly correlated or have increased betweensubject variability will receive higher weight. The main limitations of the nondynamic approach are that we do not use any information from the withinsubject variability over time. If we do not account for variability over time we may over (under) inflate the weights. In the dynamic factor model, the estimated latent variable is a weighted average of observed values from all time points. Concurrent values are weighted higher than observations further back into the past which will be discounted exponentially. The rate of discount will depend on the variability of each observed measure over time. For example, in the dynamic factor model, past observations of measures that are stable over time will be discounted less. Koopman and Harvey (2003) provide a general description of the weighting schemes for the model defined by equations (2.1) and (2.2) [42].
The dynamic factor model can be extended to allow for observed variables loading to multiple factors or for studies where participants may be clustered due to familial or other relationship. The current model is applicable to data with short temporal component and unequally spaced observations. This is a particular strength of the estimation algorithm, since most of the observational studies on cognition have these specific characteristics. A limitation of the current study is that the estimated factor is not validated with changes in biomarkers, such as volumetric data from MRI scans. Additionally, even though NACC battery is well validated and we consider the tests which load to specific domains as known, this may not be true in other applications. Another limitation of this current study is the use of aggregate scores for each test rather than the scores of each specific item used in each test. Crane et al. (2008) show that it is advantageous to use the item scores to derive latent factors in longitudinal studies [43]. Unfortunately, in the NACC study, as in many large studies, the data for items are not readily available for all participants. The methodology presented in this paper can be easily applied in most large studies where only the aggregate scores for each test are available. The dynamic factor model is particularly useful when we are interested in finding differences in the rate of cognitive change between groups. This advantage can be used in future observational studies researching the heterogeneity in rates of progression of MCI and AD patients. The DFM model can be easily used in searching specific thresholds above which the risk of conversion to MCI increases. This is especially important for designing future clinical trials that need to identify healthy participants at high risk of significant decline in cognition.
Acknowledgement
This study was supported by U01 AG016976 through a Junior Investigator Grant and by a SASIIF Grant award.
References

Snyder PJ, Jackson CE, Petersen RC, Khachaturian AS, Kaye J, et al. (2011) Assessment of cognition in mild cognitive impairment: A comparative study. Alzheimers Dementia 7: 338355.

Hayden KM, Jones RNS, Zimmer C, Plassman BL, Browndyke JN, et al. (2011) Factor structure of the national alzheimer's coordinating centers uniform dataset neuropsychological battery: An evaluation of invariance between and within groups over time. Alzheimer Dis Assoc Disord 25: 128137.

Johnson JK, Gross AL, Pa J, McLaren DG, Park LQ, et al. (2012) Longitudinal change in neuropsychological performance using latent growth models: a study of mild cognitive impairment. Brain Imaging and Behavior 6: 540550.

Proust C, JacqminGadda H, Taylor JM G, Ganiayre J, Commenges D (2006) A nonlinear model with latent process for cognitive evolution using multivariate longitudinal data. Biometrics 62: 10141024.

Locascio JJ, Atri A (2011) An overview of longitudinal data analysis methods for neurological research. Dement Geriatr Cogn Dis Extra 1: 330357.

Park LQ, Gross AL, McLaren D, Pa J, Johnson JK, et al. (2012) Confirmatory factor analysis of the ADNI neuropsychological battery. Brain Imaging and Behavior 6: 528539.

Engle R, Watson M (1981) A onefactor multivariate time series model of metropolitan wage rates. Journal of the American Statistical Association 76: 774781.

Molenaar P (1985) A dynamic factor model for the analysis of multivariate time series. Psychometrika 50: 181202.

Nesselroade JR, McArdle JJ, Aggen SH, Meyers J M (2002) Dynamic factor analysis models for representing process in multivariate timeseries. In: DS Moskowitz, SL Hershberger. Modeling intraindividual variability with repeated measures data: Methods and applications. Lawrence Erlbaum Associates Publishers: 235265.

Nesselroade JR, Molenaar PC (2004) Applying dynamic factor analysis in behavioral and social science research. In: David Kalpan, The Sage handbook of quantitative methodology for the social sciences. Sage publications, USA 33544.

Chow SM, Nesselroade JR, Shifren K, McArdle JJ (2004) Dynamic structure of emotions among individuals with parkinson's disease. Structural Equation Modeling 11: 560582.

Gu F, Preacher K, Ferrer E (2014) A state space modeling approach to mediation analysis. Journal of Educational and Behavioral Statistics 39: 117143

Chow SM, Hamaker EL, Fujita F, Boker SM (2009) Representing timevarying cyclic dynamics using multiplesubject statespace models. The British Psychological Society 62: 683716.

Molenaar P, Gooijer JGD, Schmitz B (1992) Dynamic factor analysis of nonstationary multivariate time series. Psychometrika 57: 333349.

Aguilar O, Huerta G, Prado R, West M (1998) Bayesian inference on latent structure in time series. In: JM Bernardo, JO Berger, AP Dawid, AFM Smith, Bayesian statistics 6. Oxford University press USA 116.

Rue H, Martino S, Chopin N (2009) Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71: 319392.

Zuur AF, Fryer RJ, Jolliffe IT, Dekker R, Beukema JJ (2003) Estimating common trends in multivariate time series using dynamic factor analysis. Environmetrics 14: 665685.

Shumway RH, Stoffer DS (1982) An approach to time series smoothing and forecasting using the em algorithm. Journal of Time Series Analysis 3: 253264.

Liu C, Rubin D (1998) Maximum likelihood estimation of factor analysis using the ecme alogorithm with complete and incomplete data. Statistica Sinica 8: 729747.

Koopman S (1993) Disturbance smoother for state space models. Biometrika 80: 117126.

Durbin J, Koopman S (2001) Time series analysis by state space methods (No. 24). Oxford, UK. Oxford University Press.

Kohn, R, Ansley C (1989) A fast algorithm for signal extraction, influence and crossvalidation in state space models. Biometrika 76: 6579.

deJong P (1991) The diffuse Kalman fiter. The Annals of Statistics 19: 10731083.

Zirogiannis N, Tripodis Y (2014) Dynamic factor analysis for short panels: Estimating performance trajectories for water utilities. In Agricultural and Applied Economics Association Annual Meeting Proceedings, Minneapolis, MN.

Molenaar P, Campbell C (2009) The new personspecific paradigm in psychology. Current Directions in Psychological Science 18: 112117.

Hekler EB, Buman MP, Poothakandiyil N, Rivera DE, Dzierzewski JM, et al. (2013) Exploring behavioral markers of longterm physical activity maintenance a case study of system identi cation modeling within a behavioral intervention. Health Education & Behavior 40: 51S62S.

Dempster A, Laird N, Rubin D (1977) Maximum likelihood from incomplete data via EM algorithm. Journal of the Royal Statistical Society 39: 138.

Meng X, Dyk DV (1997) The EM algorithman old folksong sung to a fast new tune. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 59: 511567.

McLachlan G, Peel D (2000) Finite mixture models. WileyInterscience.

Rubin D, Thayer D (1982) EM algorithms for ML factor analysis. Psychometrika 47: 6976.

Kalman R (1960) A new approach to linear filtering and prediction problems. Journal of Basic Engineering 82: 3545.

McLachlan G, Krishnan T (2008) The EM algorithm and extensions (2 Edition). A John Wiley & Sons, Inc., Publication, New Jersey.

Harvey A (1990) The econometric analysis of time series (2 Edition). The MIT Press, USA.

Doornik, Jurgen A (2009) An objectoriented matrix programming language ox 6. Timberlake Consultants Ltd.

Doz C, Giannone D, Reichlin L (2011) A quasimaximum likelihood approach for large, approximate dynamic factor models. Review of Economics and Statistics 94: 10141024.

Stock JH, Watson MW (2002) Forecasting using principal components from a large number of predictors. Journal of the American Statistical Association 97: 11671179.

Alzheimer's Association (2014) 2014 Alzheimer's disease facts and figures. Alzheimer's & Dementia 10.

Brookmeyer R, Johnson E, ZieglerGraham K, Arrighi HM (2007) Forecasting the global burden of alzheimer's disease. Alzheimers Dement 3: 186191.

Duara R, Loewenstein DA, Greig M, Acevedo A, Potter E, et al. (2010) Reliability and validity of an algorithm for the diagnosis of normal cognition, MCI and dementia: Implications for multicenter research studies. Am J Geriatr Psychiatry 18: 363370.

Behl P, Stefurak TL, Black SE (2005) Progress in clinical neurosciences: cognitive markers of progression in alzheimer's disease. Can J Neurol Sci 32: 140151.

Morris JC, Weintraub S, Chui HC, Cummings J, Decarli C, et al. (2006) The uniform data set (UDS): Clinical and cognitive variables and descriptive data from alzheimer disease centers. Alzheimer Dis Assoc Disord 20: 210216.

Koopman SJ, Harvey AC (2003) Computing observation weights for signal extraction and filtering. Journal of Economic Dynamics and Control 27: 13171333.

Crane PK, Narasimhalu K, Gibbons LE, Mungas DM, Haneuse S, et al. (2008) Item response theory facilitated cocalibrating cognitive tests and reduced bias in estimated rates of decline. Journal of Clinical Epidemiology 61: 10181027.