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

**|**HHS Author Manuscripts**|**PMC3621129

Formats

Article sections

Authors

Related links

Neuroimage. Author manuscript; available in PMC 2014 May 15.

Published in final edited form as:

Published online 2013 January 26. doi: 10.1016/j.neuroimage.2013.01.034

PMCID: PMC3621129

NIHMSID: NIHMS440022

Yimei Li,^{a} John H. Gilmore,^{b} Dinggang Shen,^{c,}^{e} Martin Styner,^{b} Weili Lin,^{c,}^{e} and Hongtu Zhu^{d,}^{e}

The readers are welcome to request reprints from Dr. Hongtu Zhu. Email: ude.cnu.soib@uhzh; Phone: 919-966-7272

The publisher's final edited version of this article is available at Neuroimage

See other articles in PMC that cite the published article.

Many large-scale longitudinal imaging studies have been or are being widely conducted to better understand the progress of neuropsychiatric and neurodegenerative disorders and normal brain development. The goal of this article is to develop a multiscale adaptive generalized estimation equation (MAGEE) method for spatial and adaptive analysis of neuroimaging data from longitudinal studies. MAGEE is applicable to making statistical inference on regression coefficients in both balanced and unbalanced longitudinal designs and even twin and familial studies, whereas standard software platforms have several major limitations in handling these complex studies. Specifically, conventional voxel-based analyses in these software platforms involve Gaussian smoothing imaging data and then independently fitting a statistical model at each voxel. However, the conventional smoothing methods suffer from the lack of spatial adaptivity to the shape and spatial extent of region of interest and the arbitrary choice of smoothing extent, while independently fitting statistical models across voxels does not account for the spatial properties of imaging observations and noise distribution. To address such drawbacks, we adapt a powerful propagation**–**separation (PS) procedure to sequentially incorporate the neighboring information of each voxel and develop a new novel strategy to solely update a set of parameters of interest, while fixing other nuisance parameters at their initial estimators. Simulation studies and real data analysis show that MAGEE significantly outperforms voxel-based analysis.

Many large-scale longitudinal neuroimaging studies including the Alzeimer’s disease neuroimaging initiative and the NIH magnetic resonance imaging study of normal brain have been or are being widely conducted to better understand the progress of neuropsychiatric and neurodegenerative diseases or normal brain development (Evans and Group., 2006; Almli et al., 2007; Petersen et al., 2010; Skup et al., 2011; Meltzer et al., 2009; Kim et al., 2010). The primary goal of longitudinal neuroimaging studies is to characterize individual change in neuroimaging measurements (e.g., volumetric and morphometric measurements) over time, and the effect of some covariates (or predictors) of interest, such as diagnostic status and gender, on the individual change (Petersen et al., 2010; Evans and Group., 2006). A distinctive feature of longitudinal neuroimaging data is that neuroimaging data have a temporal order. Imaging measurements of the same individual usually exhibit positive correlation and the strength of the correlation decreases with the time separation. Ignoring temporal correlation structure in imaging measures would likely influence subsequent statistical inference, such as increase in false positive and negative errors, which may lead to misleading scientific inference (Diggle et al., 2002; Fitzmaurice et al., 2004). However, the analysis of longitudinal imaging data has been hindered by the lack of advanced tools, which effectively integrate advanced image processing and statistical tools for analyzing complex and correlated imaging data along with behavioral and clinical data.

Standard software platforms have several major limitations. Standard neuroimaging software platforms including statistical parametric mapping (SPM) (www.fil.ion.ucl.ac.uk/spm/) and FMRIB Software Library (FSL) (www.fmrib.ox.ac.uk/fsl/), among many others, cannot accurately model longitudinal data when there are more than two visits (repeated measurements) (Nichols and Waldorp, 2010). Specifically, FSL can only accommodate a univariate measure at the second level (e.g., comparing visit 2 - visit 1) and SPM, even though it models the correlation among repeated measures, unrealistically assumes that the correlation is equal over the whole brain. In contrast, proper longitudinal modeling is available in standard statistical software platforms including proc MIXED and proc GEE in SAS and *lme4* and *nlme* in R. Recently, analysis of functional neuroImages (AFNI) (afni.nimh.nih.gov/afni/) adopts the linear mixed effects modeling packages *nlme* (Pinheiro et al., 2011) and *lme4* (Bates et al., 2011) in R for longitudinal functional magnetic resonance imaging data. Moreover, the Freesurfer implements the linear mixed effects modeling in the Freesurfer’s LME Matlab tool-box (http://surfer.nmr.mgh.harvard.edu/fswiki/LinearMixedEffectsModels) (Bernal-Rusiel et al., 2013). The conventional analyses of longitudinal neuroimaging data, referred to as voxel-based analysis, may be carried out in two major steps: Gaussian smoothing the imaging data and subsequently fitting a statistical model at each voxel by using either SAS or R. As discussed below, the voxel-based analysis is generally not optimal in power and the use of Gaussian smoothing may introduce substantial bias in statistical results.

The voxel-based analysis has several major limitations. First, it is common to apply a single Gaussian kernel with the full width half maximum in the range of 8–16mm to imaging data in order to account for registration errors, to Gaussianize the data, and to integrate imaging signals from a region, rather than from a single voxel. As pointed out in (Jones et al., 2005; Zhang and Davatzikos, 2011; Zhao et al., 2012; T.Ball et al., 2012), such Gaussian smoothing method can suffer from several major drawbacks including the arbitrary choice of smoothing extent and the lack of spatial adaptivity to the shape and spatial extent of the region of interest. Thus, it is suboptimal in power. In addition, as discussed in (Li et al., 2012), directly smoothing imaging data from twin and familial studies can introduce substantial bias in estimating these factors and lead to a dramatic increase of the numbers of false positives and false negatives. Second, as pointed out in (Worsley et al., 2004; Li et al., 2011), the voxel-based analysis essentially treats all voxels as independent units in the estimation stage, and thus it does not explicitly account for the spatial properties (e.g., location and smoothness) of imaging observations.

There are several attempts to address the limitations of voxel-based analysis. In (Zhang and Davatzikos, 2011), an optimally-discriminative voxel-based analysis was proposed to determine the spatially adaptive smoothing of images, followed by applying voxel-wise group analysis. The key drawback of the optimally-discriminative voxel-based analysis is that it uses the imaging data twice for both optimal weights determination and group analysis, and thus the test statistics calculated for the group analysis do not have a simple asymptotic null distribution, such as the *t* distribution. Thus, the optimally-discriminative voxel-based analysis has to resort to permutation test to calculate the *p*-values of test statistics. However, the permutation methods are not only computational intensive, but also require the so-called complete exchangeability. Such complete exchangeability is in fact a very strong assumption, and thus the optimally-discriminative voxel-based analysis is limited to both univariate imaging measure and two-group comparisons and cannot control for other continuous covariates of interest, such as age. Moreover, the optimally-discriminative voxel-based analysis has not been extended to analyze longitudinal neuroimaging data. In (Tabelow et al., 2006, 2008; Polzehl et al., 2010), the authors generalized a powerful propagation**–**separation (PS) approach (Polzehl and Spokoiny, 2000, 2006) to develop a multiscale adaptive linear model to adaptively and spatially de-noise functional magnetic resonance images and diffusion tensor images from a single subject and analyze neuroimaging data from cross-sectional studies. Recently, in (Zhu et al., 2009; Li et al., 2011; Skup et al., 2012), a multiscale adaptive regression model and a multiscale adaptive generalized method of moments approach were developed to integrate the PS approach (Polzehl and Spokoiny, 2000, 2006) with statistical modeling at each voxel for spatial and adaptive analysis of neuroimaging data from multiple subjects. All these PS related methods, however, only allow simultaneously smoothing all parameters.

This article has two major aims. The first one is to review a class of statistical methods called generalized estimating equation (GEE) for general neuroimaging researchers. We illustrate that GEE is a powerful tool for making statistical inference on regression coefficients in both balanced and unbalanced longitudinal designs and even twin and familial studies. The second aim is to develop a multiscale adaptive generalized estimating equation (MAGEE) for the spatial and adaptive analysis of longitudinal neuroimaging data. Compared with the existing literature including (Zhu et al., 2009; Li et al., 2011; Skup et al., 2012; Polzehl et al., 2010), we make several novel contributions. (i) MAGEE integrates the PS approach with GEE, which is a semiparametric model, into a simultaneous smoothing and estimation framework, allowing adaptively smoothing images while accounting for the spatial pattern of activation regions. (ii) We develop a new novel strategy of estimation and testing hypothesis of interest in MAGEE. Specifically, the new strategy allows solely smoothing the images of a set of parameters of interest, while fixing other parameters at their initial estimates. For instance, the scientific interest of many neuroimaging studies typically focuses on the comparison of imaging measures across diagnostic groups, while controlling for age, gender, and other covariates. MAGEE allows solely smoothing the images for parameter estimates of the diagnostic effect without smoothing the images of other parameter estimates, such as age and gender. (iii) We use simulated data sets to show that the new strategy can dramatically gain statistical power in some scenarios. (iv) Theoretically, in the appendix, the adaptive estimates and test statistics of MAGEE are shown to have appropriate statistical properties. We will validate companion software for MAGEE and release it to the public.

In a typical longitudinal study, one collects a fixed number of repeated measurements on all study participants at a set of common time points. When all individuals have the same number of repeated measurements on a common set of occasions, the study is “balanced” over time. Many of the early statistical methods, such as repeated-measures analysis of variance, have been developed specifically for balanced longitudinal designs. However, in most longitudinal studies over a relatively long duration in the health sciences, some individuals almost always miss their scheduled visit or date of observation. Consequently, the sequence of observation times may vary across individuals. In that case, we call the data “unbalanced” over time.

Missing data, a ubiquitous problem in longitudinal studies, can be caused by various reasons, such as skipped assessments, bad MRI scans, or study dropout. Therefore, in practice, the longitudinal data are necessarily unbalanced and they are often called “incomplete” to emphasize the fact that an intended measurement for an individual could not be obtained. Complete case analysis, a common and simple method for handling incomplete data, focuses on all individuals with complete measurements from the analysis. This approach, however, can be highly inefficient when a large proportion of the subjects are excluded. Moreover, when the individuals with complete data are not a random sample from the target population, this approach can also seriously bias estimates of longitudinal change. Fortunately, most statistical methods for longitudinal analysis, such as GEE discussed below, accommodate incomplete data under less stringent assumptions, such as missing at random (Diggle et al., 2002; Fitzmaurice et al., 2004). A good longitudinal analysis should include serious assessment of these assumptions for the data at hand and consideration of the effects of their violation on the results of the analysis, which is beyond the score of this paper. See, for example, (Ibrahim and Molenberghs, 2009) for an exhaustive review of missing data methods in longitudinal studies.

In a typical longitudinal neuroimaging study, we observe repeated imaging and clinical measures from *n* subjects. Let *m _{i}* denote the total number of time points and

The GEE method for longitudinal data has 2 main components: (i) a mean model for the mean response and its dependence on covariates; and (ii) a working covariance model for the covariance among repeated measures. GEE allows characterization and comparison of changes in the imaging measure of interest over time, complex models for the covariance, and accommodation of incomplete data. GEE can also handle unbalanced data, accommodate continuous and discrete covariates, and model the covariance in a parsimonious way (Liang and Zeger, 1986; Diggle et al., 2002). Moreover, GEE is also free of distributional assumption. For notational simplicity, we temporarily drop voxel *d* from our notation.

The mean model of GEE is to characterize the trajectories of imaging measure over time and their association with inter-individual differences in selected covariates (e.g., diagnostic group, and gender). Such trajectories can be linear or nonlinear. Specifically, the mean model of GEE is given by

$$E\{{\mathbf{y}}_{i}({t}_{ij})\mid {\mathbf{x}}_{ij}\}={\mu}_{ij}=\mu ({\mathbf{x}}_{ij},\mathit{\beta})\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i=1,\dots ,n;\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}j=1,\dots ,{m}_{i},$$

(1)

where ** β** is a

$$\mu ({\mathbf{x}}_{ij},\mathit{\beta})={\beta}_{1}+{\beta}_{2}{t}_{ij}+{\beta}_{3}{t}_{ij}^{2}+{\beta}_{4}{\text{g}}_{i}+{\beta}_{5}{\text{g}}_{i}{t}_{ij}+{\beta}_{6}{\text{g}}_{i}{t}_{ij}^{2},$$

(2)

where ** β** = (

Results of Real Longitudinal Data I: the raw line curves (in black) and fitted growth curves (in red) from MAGEE-A at a selected voxel obtained from the genu (panel (A)), splenium (panel (B)), optic radiation (panel (C)), and cerebral peduncle (panel **...**

The working covariance model of GEE is to characterize the correlation among repeated measurements on an individual and heterogeneous variability, which are two common features of longitudinal data. That is, knowledge of the value of the response on one occasion provides information about the likely value of the response on a future occasion and the variance of the response changes over the duration of the study. Specifically, the covariance matrix of **Y*** _{i}* = (

$$E\{[{\mathbf{Y}}_{i}-E({\mathbf{Y}}_{i})]{[{\mathbf{Y}}_{i}-E({\mathbf{Y}}_{i})]}^{T}\}={V}_{i}(\mathit{\theta})={A}_{i}^{1/2}(\mathit{\beta},\mathit{\gamma}){C}_{i}(\mathit{\alpha}){A}_{i}^{1/2}(\mathit{\beta},\mathit{\gamma})$$

(3)

for *i* = 1, …, *n*, where *θ* = (** β**,

As an illustration, we consider several covariance structures (3) for univariate responses. Some additional parametric models for covariance structure can be found in Appendix A. Without loss of generality, we consider the homogeneous variances of **y*** _{i}*(

- exchangeable: Corr(
**y**(_{i}*t*),_{ij}**y**(_{i}*t*)) =_{ik}*ρ*for*j*≠*k*; - autoregressive: Corr(
**y**(_{i}*t*),_{ij}**y**(_{i}*t*)) =_{ik}*ρ*^{|}^{j}^{−}^{k}^{|}for*j*≠*k*; *m***–**dependent: Corr(**y**(_{i}*t*),_{ij}**y**(_{i}*t*)) =_{ik}*ρ*_{|}_{j}_{−}_{k}_{|}for |*j*−*k*| ≤*m*;- unstructured: Corr(
**y**(_{i}*t*),_{ij}**y**(_{i}*t*)) =_{ik}*ρ*for_{jk}*j*≠*k*.

Based on a specific covariance structure of **y*** _{i}*(

$${A}_{i}(\mathit{\beta},\mathit{\gamma})={\sigma}_{y}^{2}{\mathbf{I}}_{{m}_{i}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{C}_{i}({\alpha}_{1})=(1-{\alpha}_{1}){\mathbf{I}}_{{m}_{i}}+{\alpha}_{1}{\mathbf{1}}_{{m}_{i}}{\mathbf{1}}_{{m}_{i}}^{T},$$

(4)

where **1*** _{mi}* is an

There are several advantages and drawbacks associated with the four correlation structures. The *exchangeable* working correlation structure is well known as *compound symmetry* in the longitudinal literature. The exchangeable correlation has simple interpretation, but it may be not reasonable for longitudinal studies with more than three measurements. The *auto-regressive* AR(1) correlation structure assumes the decreasing correlation as the distance between two measures increases. The *m***–***dependent* correlation structure assumes the zero correlation when two measures are *m* steps away. Both the auto-regressive and *m***–**dependent correlation structures are appealing for equally spaced data, less so for unequally spaced data. The *unstructured* working correlation structure leaves the correlation matrix completely unspecified and has *m _{i}*(

In most longitudinal studies, our primary interest focuses on making inference on ** β** or subcomponents of

- Given estimates of
^{(}^{s}^{)}and^{(}^{s}^{)}, one calculates^{(}^{s}^{+1)}as the solution of the following GEE given by$$\sum _{i=1}^{n}{U}_{i}{(\mathit{\beta})}^{T}{V}_{i}{(\mathit{\beta},{\widehat{\mathit{\alpha}}}^{(s)},{\widehat{\mathit{\gamma}}}^{(s)})}^{-1}\{{\mathbf{Y}}_{i}-{\mu}_{i}(\mathit{\beta})\}=\mathbf{0},$$(5)where*μ*(_{i}) = (*β**μ*(**x**_{i}_{1},)*β*, …,^{T}*μ*(**x**,_{imi})*β*)^{T}and^{T}*U*(_{i}) =*β**μ*(_{i})/*β**β*. Specifically, one sets^{T}^{(0)}=^{(}^{s}^{)}and then updates^{(}^{k}^{)}according to a Newton-type algorithm until ||^{(}^{k}^{+1)}−^{(}^{k}^{)}||_{2}is smaller than a fixed small number*ε*, say 10^{−4}. One sets^{(}^{s}^{+1)}=^{(}^{k}^{+1)}. - Given
^{(}^{s}^{+1)}, one computes ${\widehat{\mathbf{r}}}_{ij}^{(s+1)}={\mathbf{y}}_{i}({t}_{ij})-\mu ({\mathbf{x}}_{ij},{\widehat{\mathit{\beta}}}^{(s+1)})$ for all*i*,*j*and uses them to construct some estimates of (,*α*), denoted by (*γ*^{(}^{s}^{+1)},*γ*^{(}^{s}^{+1)}). See Appendix B for some estimation methods for (,*α*).*γ* - Iterate the previous two steps until convergence.
- Obtain the final estimator
**= (****,****,****) and then calculate***V*(_{i}),*θ*=_{ij}**y**(_{i}*t*) −_{ij}*μ*(**x**,_{ij}**) and the sandwich estimator of the covariance matrix of****, denoted by Σ**(_{n}**), which is given by**$${(\sum _{i=1}^{n}{\widehat{U}}_{i}^{T}{\widehat{V}}_{i}^{-1}{\widehat{U}}_{i})}^{-1}(\sum _{i=1}^{n}{\widehat{U}}_{i}^{T}{\widehat{V}}_{i}^{-1}{\widehat{\mathbf{r}}}_{i}{\widehat{\mathbf{r}}}_{i}^{T}{\widehat{V}}_{i}^{-1}{\widehat{U}}_{i}){(\sum _{i=1}^{n}{\widehat{U}}_{i}^{T}{\widehat{V}}_{i}^{-1}{\widehat{U}}_{i})}^{-1},$$(6)where=_{i}**Y**−_{i}*μ*(_{i}**),***Û*=_{i}*U*(_{i}**), and**=_{i}*V*(_{i}**).**

The ** based on GEE (5) has three attractive properties. Firstly, **** can be almost as efficient as the maximum likelihood estimates of **** β** in many practical applications provided that

In longitudinal studies, ** β**(

$${\mathit{\beta}}_{I}(d)={({\beta}_{5}(d),{\beta}_{6}(d))}^{T},\phantom{\rule{0.38889em}{0ex}}{\mathit{\beta}}_{N}(d)={({\beta}_{1}(d),\dots ,{\beta}_{4}(d))}^{T},\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.16667em}{0ex}}({q}_{I},{q}_{N})=(2,4).$$

Moreover, the components of (** α**(

We propose a weighted generalized estimating equations method, referred to as weighted GEE-A, to spatially and adaptively update {*β** _{I}* (

$$\sum _{{d}^{\prime}\in B(d,h)}\omega (d,{d}^{\prime};h)=1\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\omega (d,{d}^{\prime};h)\ge 0\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}\text{all}\phantom{\rule{0.38889em}{0ex}}h\ge 0.$$

weighted GEE-A is based on a set of weighted GEEs, denoted by *G _{n}*(

$$\sum _{{d}^{\prime}\in B(d,h)}\omega (d,{d}^{\prime};h)\sum _{i=1}^{n}{U}_{iI}{({\mathit{\beta}}_{I}(d),{\widehat{\mathit{\beta}}}_{N}({d}^{\prime}))}^{T}{V}_{i}{(\widehat{\mathit{\theta}}({d}^{\prime}))}^{-1}{\mathbf{e}}_{i}({\mathit{\beta}}_{I}(d),{\widehat{\mathit{\beta}}}_{N}({d}^{\prime}))=\mathbf{0},$$

(7)

where **e*** _{i}*(

$${G}_{n}({\widehat{\mathit{\beta}}}_{I}(d,h);\omega ,h)=\mathbf{0}.$$

(8)

A good *ω*(*d*, *d*′; *h*) plays a critical role in preventing oversmoothing the estimates of *β** _{I}* (

We present test statistic based on the weighted GEE (7) at each *d*
for a fixed radius *h*. Our choice of which hypotheses to test is motivated by either a comparison of brain structure (or function) across diagnostic groups or the detection of a change in brain structure (or function) across time (Skup et al., 2011; Meltzer et al., 2009; Kim et al., 2010). These questions usually can be formulated as testing hypotheses about *β** _{I}* (

$${H}_{0,\beta}:R{\mathit{\beta}}_{I}(d)={\mathbf{b}}_{0}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{vs}.\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{H}_{1,\beta}:R{\mathit{\beta}}_{I}(d)\ne {\mathbf{b}}_{0},$$

(9)

where *R* is a *r* × *q _{I}* matrix and

$$R=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\mathbf{b}}_{0}={(0,0)}^{T}.$$

We test the null hypothesis *H*_{0,}* _{β}* using a Wald test statistic given by

$${W}_{\beta}(d,h)={[R{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathbf{b}}_{0}]}^{T}{[R{\mathrm{\sum}}_{n}({\widehat{\mathit{\beta}}}_{I}(d,h)){R}^{T}]}^{-1}[R{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathbf{b}}_{0}],$$

(10)

where Σ* _{n}*(

We develop a PS procedure based on multiscale adaptive generalized estimating equations, referred as to MAGEE-A, by integrating weighted GEE-A and the PS procedure proposed in (Polzehl and Spokoiny, 2000, 2006). Since PS and the choice of its associated parameters have been discussed in details in (Polzehl et al., 2010; Li et al., 2011), we briefly mention them here.

- At each
*d*, the PS procedure evolves along a sequence of nested spheres with increasing radii*h*as follows:_{s}$$\begin{array}{lllllll}{h}_{0}=0\hfill & <\hfill & \hfill {h}_{1}\hfill & <\hfill & \dots \hfill & <\hfill & {h}_{S}={r}_{0};\hfill \\ B(d,{h}_{0})=\{d\}\hfill & \subset \hfill & \hfill B(d,{h}_{1})\hfill & \subset \hfill & \dots \hfill & \subset \hfill & B(d,{h}_{S}).\hfill \end{array}$$(11) - At the scale
*h*_{0}= 0, we just use the voxel-wise GEE estimator**(***d*,*h*_{0}) =**(***d*) and then we fix ((*β*_{N}*d*),(*α**d*),(*γ**d*)) at ((_{N}*d*),**(***d*),**(***d*)). - We combine all information contained in {
(_{I}*d*) :*d*} to calculate weights*ω*(*d*,*d*′;*h*_{1}) at scale*h*_{1}for all*d*. Subsequently, we utilize all data in {*B*(*d*,*h*_{1}) :*d*}, all weights {*ω*(*d*,*d*′;*h*_{1}) :*d*,*d*′ }, and weighted GEEs (7) to estimate**(***d; h*_{1}) for all*d*. In this way, we can sequentially determine*ω*(*d*,*d*′;*h*) and adaptively update_{s}**(***d*,*h*) from_{s}*h*_{0}= 0 to*h*=_{S}*r*_{0}. At the end of PS, we calculate the Wald test statistics*W*(_{β}*d*,*h*) across all voxels_{S}*d*.

A path diagram of the PS procedure for MAGEE-A is given as follows:

$$\begin{array}{ccccccc}B(d,{h}_{0})=\{d\}& \subset & B(d,{h}_{1})& \subset & \dots & \subset & B(d,{h}_{S})\\ \Downarrow & \phantom{\rule{0.16667em}{0ex}}& \Downarrow & \phantom{\rule{0.16667em}{0ex}}& \dots & \phantom{\rule{0.16667em}{0ex}}& \Downarrow \\ \widehat{\mathit{\theta}}(d)& \phantom{\rule{0.16667em}{0ex}}& \omega (d,{d}^{\prime};{h}_{1})& \phantom{\rule{0.16667em}{0ex}}& \dots & \phantom{\rule{0.16667em}{0ex}}& \omega (d,{d}^{\prime};{h}_{S}={r}_{0})\\ \Downarrow & \nearrow & \Downarrow & \nearrow & \dots & \nearrow & \Downarrow \\ {\widehat{\mathit{\beta}}}_{I}(d,{h}_{0})& \phantom{\rule{0.16667em}{0ex}}& {\widehat{\mathit{\beta}}}_{I}(d,{h}_{1})& \phantom{\rule{0.16667em}{0ex}}& \dots & \phantom{\rule{0.16667em}{0ex}}& {\widehat{\mathit{\beta}}}_{I}(d,{h}_{S})\\ \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \Downarrow \\ \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& {W}_{\beta}(d,{h}_{S})\end{array}.$$

(12)

The detailed steps of MAGEE-A are given Appendix D.

We simulated data at all *m* = 23, 232 voxels on a 88 × 88 × 5 phantom image for *n* = 80 subjects. Each slice contains the same activation region. Specifically, at each voxel *d* in
, we simulated *y _{ij}*(

$${y}_{ij}(d)={\mathbf{x}}_{ij}^{T}\mathit{\beta}(d)+{\epsilon}_{ij}(d)\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}j=1,\dots ,{m}_{i}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i=1,\dots ,n,$$

where ** β**(

We fitted GEE with the AR(1) working correlation structure and homogeneous variances and set

$${\mu}_{ij}({\mathbf{x}}_{ij},\mathit{\beta}(d))={\mathbf{x}}_{ij}^{T}\mathit{\beta}(d)={\beta}_{1}(d)+{x}_{ij2}{\beta}_{2}(d)+{x}_{ij3}{\beta}_{3}(d).$$

We used MAGEE-A to spatially and adaptively calculate the parameter estimates of ** β**(

For simplicity, we present some selected results for _{2}(*d*, *h _{s}*). Both MAGEE-A and MAGEE-B show better accuracy of parameter estimates compared with the voxel-wise GEE estimates based on the effective smoothness scales of 4 and 8mm in terms of bias (Fig. 2 (B)–(E)). Moreover, MAGEE-B slightly outperforms MAGEE-A in terms of bias and the root-mean-square error (Table 1), since

Simulation results for scenario I: Panel (A) is a selected slice of the ground truth image for simulation study of the six regions of interest with different gray levels representing *β*_{2}(*d*)=0, 0.05, 0.1, 0.2, 0.3, 0.4 respectively; Panel (B) represents **...**

Simulation results: Average Bias, RMS, SD, and RE of *β*_{2}(*d*) parameters in the six regions of interest at scale *h*_{10} for MAGEE-A, MAGEE-B, and the voxel-wise GEE methods with the heat kernel smoothing with *k* iterations, which are denoted as SGEE- **...**

We tested the hypotheses *H*_{0} : *β*_{2}(*d*) = 0 and *H*_{1} : *β*_{2}(*d*) ≠ 0 across all voxels by using *W _{β}*(

We calculated three statistics based on the two groups of clusters including the dice overlap ratio between the clusters in the ‘good’ group and the true activation regions, the cluster number, and the spatial extent (or the number of voxels) of all clusters in the ‘bad’ group. The value of dice overlap ratio is between 0 and 1 and a larger dice overlap ratio value represents better performance in detecting true positive voxels. Moreover, if the spatial extent of a cluster is smaller than a threshold, such as 10, then the cluster is not detected to be significant by the cluster based thresholding. Table 2 presents the three statistics of each method.

Simulation results: Dice Overlap Ratio (DOR), the non overlap cluster number and cluster size for testing *β*_{2}(*d*) = 0 in the six regions of interest for MAGEE-A, MAGEE-B, and the voxel-wise GEE methods with the heat kernel smoothing with *k* iterations, **...**

Inspecting Table 2 reveals several facts. (i) Both MAGEE-A and MAGEE-B have larger dice overlap ratio compared with voxel-wise GEEs with relatively large smoothing scale (≥4mm). See Fig. 4 (a) for details. (ii) Although MAGEE-A and MAGEE-B have many more non-overlap small clusters than voxel-wise GEEs methods, their spatial extent is much smaller than the standard spatial extent threshold (e.g., 10). (iii) MAGEE-B outperforms MAGEE-A in terms of dice overlap ratio and non-overlap cluster number. (iv) For voxel-wise GEEs, small smoothing scales (≤3mm) outperform large smoothing scales (≤4mm) in terms of dice overlap ratio, and vice versa in terms of non-overlap cluster number. However, as expected, the larger is smoothing scale, the larger is non-overlap cluster. This can lead to false detection of non-overlap clusters.

Simulation results for scenarios I and II for comparing MAGEE-A, MAGEE-B, and the voxel-wise GEE methods with the heat kernel smoothing with *k* iterations, which are denoted as SGEE-*k* for *k* = 1, 4, 9, 16, 25, 36, 49, and 64. Panel (a) presents Dice Overlap **...**

Figure 4 (b) and (c) presents the results for testing *β*_{2}(*d*) = 0 in all voxels in the six regions of interest for all methods by using the false-discovery rate corrected threshold, *p _{c}* = 0.05, based on 1,000 simulations. For each method, we calculated the type I error rate and power based on voxels in all regions of interest for each simulated data set. Figure 4 (b) also presents their average type I error rates, whereas Figure 4 (c) presents the detection power of voxels in each of active regions of interest. MAGEE-A and MAGEE-B have higher detection power in most regions of interest, while their type I error rates are quite small. As smoothing scale is larger than 3mm, the type I error rate of the voxel-wise GEE becomes more liberal and not well controlled. In contrast, although the type I error rate of voxel-wise GEE is small, their detection powers are poor compared with the MAGEE-A and MAGEE-B methods. In summary, the MAGEE methods outperform the voxel-wise GEEs with different smoothing scales in terms of both type I and II error rates.

We used the same setup as scenario I except that we randomly generated *β*_{3}(*d*) from *U*[−5, 5] in the prefixed seven regions of interest with *β*_{2}(*d*) ≠ 0. In this case, *β*_{3}(*d*) varies in these six regions of interest, and thus the *β*_{2} image and the *β*_{3} image do not share the same imaging pattern. We are interested in examining the effect of different patterns in the *β*_{3} image on _{2}(*d*, *h _{s}*) across all

For _{2}(*d*, *h _{S}*), MAGEE-A, which smooths

A wealth of cross-sectional diffusion tensor imaging studies has been conducted on characterizing white matter development (prenatal to adolescent stages) using various diffusion parameters, such as fractional anisotropy and radial diffusivity, in the past decade. Current diffusion tensor imaging studies involving neonates have revealed three phases in the early postnatal brain development including the rapid changes within first 12 months, the slow maturation from 12 to 24 months, and the steady state afterwards. Particularly, in white matter, neonates have significantly lower fractional anisotropy values and significantly higher mean diffusivity values compared to adults (Neil et al., 1998). These diffusion tensor imaging studies also reveal the temporal non-linearity and spatial inhomogeneity of the apparent changes in diffusion tensor imaging parameters within brain (Mukherjee and McKinstry, 2006; Schneider et al., 2004).

We used 38 subjects from a longitudinal study designed to investigate early brain development led by Dr. John Gilmore at the University of North Carolina at Chapel Hill. For each subject, diffusion-weighted images were acquired at 2 weeks, year 1, and year 2. Diffusion-weighted imaging acquisition scheme includes 18 repeated measures of six non-collinear directions, (1,0,1), (−1,0,1), (0,1,1), (0,1,−1), (1,1,0), and (−1,1,0), at a *b***–**value of 1000 s/mm^{2} and a *b* = 0 reference scan. Forty-six contiguous slices with a slice thickness of 2 mm covered a field of view of 256×256 mm^{2} with an isotropic voxel size of 2 × 2 × 2 mm^{3}. High resolution T1 weighted images were acquired using a 3D MP-RAGE sequence. We then calculated a weighted least squares estimation method to construct diffusion tensors (Basser et al., 1994; Zhu et al., 2007). All images were visually inspected before analysis to ensure no bulk motion. All diffusion tensor images (38 subjects, 3 time points each) were registered onto a template, which is a randomly selected brain diffusion tensor image of a 2-year-old subject (Yap et al., 2009).

We use fractional anisotropy images to identify the spatial patterns of white matter maturation. We fitted a GEE with AR(1) correlation structure and homogeneous variance structure and

$${\mu}_{ij}({\mathbf{x}}_{ij},\mathit{\beta}(d))={\beta}_{1}(d)+{\beta}_{2}(d){t}_{ij}+{\beta}_{3}(d){t}_{ij}^{2}+{\beta}_{4}(d){g}_{i}+{\beta}_{5}(d){t}_{ij}{g}_{i}+{\beta}_{6}(d){t}_{ij}^{2}{g}_{i}$$

for *i* = 1, …, *n* and *j* = 1, 2, 3 at each voxel of the template. We are interested in testing two sets of hypothesis including the gender and time interaction effect and then the time effect if the gender and time interaction effect is not significant. We used MAGEE-A and MAGEE-B to address such hypotheses. We also smoothed fractional anisotropy imaging data using an isotropic Gaussian kernel with full width half maximum 8mm and used voxel-wise GEE to address those hypotheses as well. Moreover, although fractional anisotropy is between 0 and 1, *μ*(**x**, ** β**(

We run the PS procedure with *c _{h}* = 1.15 and

$${\mu}_{ij}({\mathbf{x}}_{ij},\mathit{\beta}(d))={\beta}_{1}(d)+{\beta}_{2}(d){t}_{ij}+{\beta}_{3}(d){t}_{ij}^{2}+{\beta}_{4}(d){g}_{i}$$

for *i* = 1, …, *n* and *j* = 1, 2, 3 at each voxel of the template. Similarly, the findings for MAGEE-A are similar to those for MAGEE-B, whereas MAGEE-A detects more significant clusters and voxels compared with MAGEE-B at *S* = 10 (Fig. 5 (B), (C), (B′), and (C′)). As expected, both MAGEE-A and MAGEE-B outperform voxel-wise GEE based on unsmoothed and smoothed FA images in terms of the number of both significant clusters and voxels (Fig. 5 (A) and (A′)). Specifically, the voxel-wise GEE results based on the smoothed fractional anisotropy images show the obvious oversmoothing in cerebrospinal fluid and the gray matter areas, such as the ventricle (Fig. 5 (D) and (D′)).

Results of Real Longitudinal Data I: panels (A), (B), (C), and (D) respectively, show the – log_{10}(*p*) values of *W*_{β}(*d*, *h*_{s}) for testing the age effect for voxel-wise GEE for unsmoothed FA images, MAGEE-A at the 10-th iteration, and MAGEE-B **...**

We also plotted the line plots and fitted growth curves in four selected voxels in the genu, splenium, optic radiation and cerebral peduncle (Fig. 1). Different growth patterns were observed for the genu, splenium, optic radiation, and cerebral peduncle (Fig. 1). The genu and cerebral peduncle have a similar and small fractional anisotropy value at birth, whereas the genu’s fractional anisotropy increases very fast and the cerebral peduncle’s fractional anisotropy increases slowly. The splenium and optic radiation have a similar and moderate fractional anisotropy value at birth, whereas the fractional anisotropy value of the splenium and optic radiation increases relatively slow. There is a substantial variability across different anatomic regions, even though fractional anisotropy is highly significantly different between neonates and 1-year-olds for all regions of interest.

We considered the same early postnatal brain development project led by Dr. Gilmore at the University of North Carolina at Chapel Hill (Gilmore et al., 2010). A total of 49 paired twins (36 males and 62 females) were selected for an illustration. All 49 pairs were scanned as neonates within a few weeks after birth at term. Written consent was obtained from their parents before imaging acquisition. The mean gestational age at magnetic resonance scanning was 246 ± 18.3 days (range: 192 to 270 days). All infants were fed and calmed to sleep on a warm blanket with proper ear protection and they slept comfortably inside the MR scanner. None of infants was sedated during the imaging session.

We then employed a nonlinear fluid deformation based high-dimensional, unbiased atlas computation method to process all 98 diffusion tensor imaging datasets (Goodlett et al., 2009). The atlas building procedure started with an affine registration and was followed by a nonlinear registration of a set of feature images for all subjects. The feature images are the maximum eigenvalue of the Hessian of the fractional anisotropy image, which are sensitive to the geometry of white matter. Using the computed deformation fields, we warped all tensor images into the unbiased atlas space via log-euclidean based interpolation (Arsigny et al., 2006). We also averaged all the warped tensor images to create a study specific diffusion tensor imaging atlas.

Fractional anisotropy has been widely used as a measurement to assess directional organization of the brain, which is greatly influenced by the magnitude and orientation of white matter tracts. Here, fractional anisotropy images are employed to identify the spatial patterns of white matter maturation. We considered a structural equation model given by

$${y}_{ij}(d)={\beta}_{1}(d)+{\beta}_{2}(d){g}_{ij}+{\beta}_{3}(d){z}_{ij}+{a}_{ij}(d)+{c}_{i}(d)+{\epsilon}_{ij}(d)$$

(13)

for *i* = 1, …, 49 and *j* = 1, 2 at each voxel of the template, where *g _{ij}* and

Our objective is to show that MAGEE is applicable to making statistical inference on regression coefficients including *β*_{1}(*d*), *β*_{2}(*d*), and *β*_{3}(*d*) in twin studies. As an illustration, we consider the gender effect, that is *H*_{0} : *β*_{2}(*d*) = 0. We applied MAGEE-A and MAGEE-B and compared their results with those obtained from three other methods including an maximum likelihood estimation method, TwinMARM, and voxel-wise GEE based on Gaussian smoothed images with full width half maximum 8mm. For MAGEE-A and MAGEE-B, we fitted weighted GEE with the compound-symmetry correlation and *μ*(**x*** _{ij}*,

Inspecting Fig. 6 (A)–(J) reveals that the images of parameter estimator obtained from all five methods have similar pattern. All multiscale adaptive methods including TwinMARM, MAGEE-A, and MAGEE-B significantly outperform the maximum likelihood estimation method in terms of the spatial smoothness of significant clusters and the magnitude of standard deviations (Fig. 6). Comparing Fig. 6 (E) with Fig. 6 (A)–(D) reveals some subtle differences between voxel-wise GEE based on the smoothed images and all other methods. This may confirm the biases introduced by the use of Gaussian kernel as discussed in Li et al. (2012).

We have reviewed the method of GEE for analyses of repeated-measures data. GEE is powerful for handling unbalanced designs and missing data obtained from longitudinal, twin, and familial studies, accommodating continuous and discrete covariates and their interactions, and releasing distributional assumptions and complex models for the covariance. Specifically, we discuss how to set up the mean model and the working covariance model of GEE. The mean model is used to directly characterize individual trajectory in repeated measurements as linear and nonlinear functions of time, while adjusting for the effect of other predictors, such as diagnostic status and gender, on the individual trajectory. Although the linear mean models are generally satisfactory approximations for most neuroimaging applications, there are many cases, such as growth curve and dose-response relationships, when an empirical indicated or a theoretically justified nonlinear mean model is more appropriate. For instance, since growth from birth to maturity in human subjects typically is nonlinear in nature, a nonlinear mean model should be used. The logistic and Gompertz nonlinear models described in Appendix A have been widely used to characterize rapid growth shortly after birth, pronounced growth during puberty, and a leveling off sometime before adulthood. In many applications, if a simple nonlinear mean model is unavailable at the beginning, it is common to fit a nonparametric mean model to longitudinal data and then use a simple mean model to approximate the fitted nonparametric mean model (Wang, 2003; Wu and Zhang, 2006; Wang, 1998).

The working covariance model of GEE is used to model and understand the likely sources of random variation in longitudinal data. As discussed in Diggle et al. (2002), a useful working covariance model should include at least three qualitatively different sources of random variation: (i) random effects, (ii) serial correlation, and (iii) measurement error. Various strategies have been developed to incorporate them into specific models (Diggle et al., 2002). For instance, we have discussed four commonly used correlation structures, such as exchangeable, in Section 2.4, and several parametric models for covariance structure in Appendix A. Although these parametric models may be sufficient for longitudinal data with few repeated measures, most parametric models discussed here are not appealing for sparsely and irregularly longitudinal data. Recently, several nonparametric covariance models have been proposed to deal with such issue (Yao et al., 2005; Ramsay and Silverman, 2005). Similar to the mean model, one may fit a nonparametric covariance model to longitudinal data and then use a simple covariance model to approximate the nonparametric model.

We have developed MAGEE for the spatial and adaptive analyses of longitudinal neuroimaging data. MAGEE is essentially a locally adaptive and spatial smoothing method and is adaptive to the spatial pattern and extent of each activation region for each regression coefficient map. Such adaptive and spatial property is very appealing from theoretical and practical perspectives. Specifically, according to the matched filter theorem, the size of the optimal filter should match the size of target signal, while accounting for noise distribution. When there are multiple effect regions with different spatial patterns and extent, multiple filters with different shapes and kernel sizes should be used. Moreover, MAGEE as a weighted GEE method explicitly incorporate the mean model and the working covariance model of GEE. In contrast, the single-filter and multi-filter methods are independent of the model assumptions of GEE (Jones et al., 2005; Poline and Mazoyer, 1994; Zhao et al., 2012; T.Ball et al., 2012; Siegmund and Worsley, 1995). Thus, the smoothed imaging data may not follow the assumed model assumptions of GEE, which can lead to biases for the analysis of longitudinal neuroimaging data. Such biases can be substantial for all kinds of neuroimaging data from longitudinal, twin, and familial studies. See Li et al. (2012) for detailed discussions on biases introduced by directly smoothing twin imaging data. Actually, the same discussions are valid for longitudinal neuroimaging data when either the nonlinear mean model is valid or the covariance model is the primary problem of interest. Therefore, one should interpret the statistical findings obtained from the voxel-based analysis based on directly smoothed longitudinal neuroimaging data with great caution.

We have used both simulations and real imaging data to demonstrate that MAGEE outperforms the voxel-wise GEE method coupled with a single Gaussian kernel. As shown in our simulations, the kernel size in the single Gaussian kernel can have a substantial impact on brain mapping results. Specifically, as the kernel width varies from 1- to 8mm, the dice overlap ratio and the type I and II error rates of voxel-wise GEE change substantially. As shown in the first real imaging data, the use of the single smoothing kernel can lead to the false positive results of significant diffusion property changes in cerebrospinal fluid and the gray matter areas. Our results are also consistent with previous results on single-filter and multi-filter analyses for various neuroimaging data including positron emission tomography, diffusion tensor imaging, functional magnetic resonance imaging, and cortical morphometry (Jones et al., 2005; Poline and Mazoyer, 1994; Zhao et al., 2012; T.Ball et al., 2012; Siegmund and Worsley, 1995). As shown in the second real imaging data, the commonly used Gaussian kernel for smoothing imaging data can introduce biases for the analysis of twin imaging data.

Several important issues need to be addressed in future research. Firstly, the MAGEE procedure is solely powerful for detecting relatively large effect regions, which are smooth interiorly and consistent across subjects. In practice, however, the extent and location of effect regions may vary dramatically across subjects due to both registration error and population heterogeneity. Therefore, it is important and interesting to model population heterogeneity, while accounting for registration error. Secondly, although we focus on parametric growth curves, it is interesting to develop more flexible nonparametric growth curve models, which are important for sparsely and irregularly longitudinal studies. Developing multiscale adaptive methods for such non-parametric models faces up with many new challenges both computationally and theoretically. Thirdly, more research is needed for optimizing the choices of the parameters in MAGEE and for incorporating other edge-preserving local smoothing methods into MAGEE (Qiu, 2005; Qiu and Mukherjee, 2010; Mukherjee and Qiu, 2011). Fourthly, we will extend MAGEE from simple longitudinal studies to more complex longitudinal twin/familial studies. Fifthly, it is interesting and important to treat the whole image as a single piecewisely smoothed function, instead of a collection of isolated voxels, and then develop new statistical models to directly model such functions from cross-sectional and longitudinal studies (Zhu et al., 2011; Ramsay and Silverman, 2005; Yao et al., 2005; Greven et al., 2010).

- Develop MAGEE for longitudinal neuroimaging data.
- Characterize the development of white matter diffusivities.
- Evaluate the finite sample performance of MAGEE.
- MAGEE significantly outperforms VBA.

We may consider other parametric models for mean structure. The *μ*(**x*** _{ij}*,

$$\mu ({\mathbf{x}}_{ij},\mathit{\beta})=\frac{{\beta}_{1}}{1+{\beta}_{2}exp(-{\beta}_{3}{t}_{ij})}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\mu ({\mathbf{x}}_{ij},\mathit{\beta})={\beta}_{1}exp(-{\beta}_{2}exp(-{\beta}_{3}{t}_{ij})),$$

(14)

where **x*** _{ij}* = (

The growth curve pattern for different models. Red: Gomportz growth model with *μ*(**x**, *β*) = 3.6 exp(−2 exp(−2*t*)); Blue: Quadratic growth model for female with *μ*(**x**, *β*) = 0.1 + 0.8*t* − 0.05*t*^{2}; Magenta: **...**

We may consider other parametric models for covariance structure. See (Diggle et al., 2002; Pourahmadi, 2011; Zimmerman and Nunez-Anton, 2001) for overviews of different parametric models for covariance structure. For instance, one may consider a heteroscedastic model for the variances of **y*** _{i}*(

$$\text{Corr}({\mathbf{y}}_{i}({t}_{ij}),{\mathbf{y}}_{i}({t}_{ik}))=\rho (\mid {t}_{ij}-{t}_{ik}\mid ),$$

(15)

where *ρ*(*u*) is a known function of *u*. Two popular choices of *ρ*(*u*) are the exponential correlation model and the Gaussian correlation model, which are, respectively, given by

$$\rho (u)=exp(-{\gamma}_{1}\mid u\mid )\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\rho (u)=exp(-{\gamma}_{1}{u}^{2})$$

(16)

for some *γ*_{1} > 0. Based on these more general forms of covariance structure, one can derive the explicit form of **Y*** _{i}* in (3).

To estimate (** α**,

$$\sum _{i=1}^{n}\frac{\partial \text{Vecs}({V}_{i}({\widehat{\mathit{\beta}}}^{(s+1)},\mathit{\alpha},\mathit{\gamma}))}{\partial (\mathit{\alpha},\mathit{\gamma})}\{{\mathcal{V}}_{Y,i}({\widehat{\mathit{\beta}}}^{(s+1)})-\text{Vecs}({V}_{i}({\widehat{\mathit{\beta}}}^{(s+1)},\mathit{\alpha},\mathit{\gamma}))\}=\mathbf{0}.$$

(17)

In many cases, the dimension of (** α**,

$$\begin{array}{l}{\widehat{\mathit{\alpha}}}^{(s+1)}=\frac{2}{{\sum}_{i=1}^{n}{m}_{i}({m}_{i}-1)\phi}\sum _{i=1}^{n}\sum _{j<k}{\widehat{\mathbf{r}}}_{ij}^{(s+1)}{\widehat{\mathbf{r}}}_{ik}^{(s+1)},\\ {\widehat{\sigma}}_{y}^{(s+1)2}=\frac{1}{{\sum}_{i=1}^{n}{m}_{i}}\sum _{i=1}^{n}\sum _{j=1}^{{m}_{i}}{\widehat{\mathbf{r}}}_{ij}^{(s+1)2}.\end{array}$$

For some other covariance structures, the exact expression of their moment estimates of (** α**,

We derive an approximation of the covariance matrix of * _{I}* (

We need to introduce some notation. Let *θ** _{N}* (

$${G}_{n}({\mathit{\beta}}_{I}(d);\omega ,h)={G}_{n}({\mathit{\beta}}_{I}(d),{\widehat{\mathit{\theta}}}_{N}({d}^{\prime};B(d,h));\omega ,h).$$

We define

$$\begin{array}{l}{U}_{iN\ast}(d,{d}^{\prime})=\frac{\partial {\mu}_{i}({\mathit{\beta}}_{I\ast}(d),{\mathit{\beta}}_{N\ast}({d}^{\prime}))}{\partial {\mathit{\beta}}_{N}(d)},{U}_{iI\ast}(d,{d}^{\prime})=\frac{\partial {\mu}_{i}({\mathit{\beta}}_{I\ast}(d),{\mathit{\beta}}_{N\ast}({d}^{\prime}))}{\partial {\mathit{\beta}}_{I}(d)},\\ {\mathbf{F}}_{\ast}(d)={n}^{-1}\sum _{i=1}^{n}{U}_{i}{({\mathit{\beta}}_{\ast}(d))}^{T}{V}_{i}{({\mathit{\theta}}_{\ast}(d))}^{-1}{U}_{i}({\mathit{\beta}}_{\ast}(d)),\\ {\mathbf{F}}_{\ast}{(d)}^{-1}=\left(\begin{array}{cc}{\mathbf{F}}_{\ast}^{II}(d)& {\mathbf{F}}_{\ast}^{IN}(d)\\ {\mathbf{F}}_{\ast}^{NI}(d)& {\mathbf{F}}_{\ast}^{NN}(d)\end{array}\right),\end{array}$$

(18)

in which we have partitioned **F**_{*}(*d*)^{−1} according to the partition of ** β**(

The derivations consist of five steps as follows.

- Step 1 is to derive an asymptotic expansion of
*G*(_{n}(_{I}*d*,*h*);*ω*,*h*). By using the Taylor’s series expansion, we have$$\mathbf{0}={G}_{n}({\widehat{\mathit{\beta}}}_{I}(d,h);\omega ,h)\approx {G}_{n}({\mathit{\beta}}_{I\ast}(d);\omega ,h)+\frac{\partial {G}_{n}({\mathit{\beta}}_{I\ast}(d);\omega ,h)}{\partial {\mathit{\beta}}_{I}(d)}\{{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathit{\beta}}_{I\ast}(d)\},$$(19)$$\begin{array}{l}{G}_{n}({\beta}_{I\ast}(d);\omega ,h)={G}_{n}({\mathit{\beta}}_{I\ast}(d),{\widehat{\mathit{\theta}}}_{N}({d}^{\prime};B(d,h));\omega ,h)\approx \\ {G}_{n}({\mathit{\beta}}_{I\ast}(d),{\mathit{\theta}}_{N\ast}({d}^{\prime};B(d,h));\omega ,h)+\\ \sum _{{d}^{\prime}\in B(d,h)}\frac{\partial {G}_{n}({\mathit{\beta}}_{I\ast}(d),{\mathit{\theta}}_{N\ast}({d}^{\prime};B(d,h));\omega ,h)}{\partial {\mathit{\theta}}_{N}({d}^{\prime})}\{{\widehat{\mathit{\theta}}}_{N}({d}^{\prime})-{\mathit{\theta}}_{N\ast}({d}^{\prime})\},\end{array}$$(20)where*θ*_{N}_{*}(*d*) and*θ*_{N}_{*}(*d*′;*B*(*d*,*h*)) are, respectively, the true value of*β*(_{N}*d*) and that of*β*(_{N}*d*′;*B*(*d*,*h*)). - Step 2 is to characterize the asymptotic properties of the first-order derivative of
*G*(_{n}**β**_{I}_{*}(*d*),*θ*_{N}_{*}(*d*′;*B*(*d*,*h*));*ω*,*h*) with respect to*θ*(_{N}*d*′). It follows from the law of large number and (1) that$$-{n}^{-1}\frac{\partial {G}_{n}({\mathit{\beta}}_{I\ast}(d),{\mathit{\theta}}_{N\ast}({d}^{\prime};B(d,h));\omega ,h)}{\partial {\mathit{\beta}}_{N}({d}^{\prime})}={\mathbf{E}}_{IN\ast}({d}^{\prime};d,h)\approx \omega (d,{d}^{\prime};h){n}^{-1}\sum _{i=1}^{n}{U}_{iI\ast}{(d,{d}^{\prime})}^{T}{V}_{i}{({\mathit{\theta}}_{\ast}({d}^{\prime}))}^{-1}{U}_{iN}(d,{d}^{\prime}),$$(21)$${n}^{-1}\frac{\partial {G}_{n}({\mathit{\beta}}_{I\ast}(d),{\mathit{\theta}}_{N\ast}({d}^{\prime};B(d,h));\omega ,h)}{\partial (\mathit{\gamma}({d}^{\prime}),\mathit{\alpha}({d}^{\prime}))}\approx \mathbf{0},$$(22)$$-{n}^{-1}\frac{\partial {G}_{n}({\mathit{\beta}}_{I\ast}(d);\omega ,h)}{\partial {\mathit{\beta}}_{I}(d)}={\mathbf{E}}_{II\ast}(d,h)\approx {n}^{-1}\sum _{{d}^{\prime}\in B(d,h)}\omega (d,{d}^{\prime};h)\sum _{i=1}^{n}{U}_{iI}{(d,{d}^{\prime})}^{T}{V}_{i}{({\mathit{\theta}}_{\ast}({d}^{\prime}))}^{-1}{U}_{iI}(d,{d}^{\prime}),$$(23)where ≈ denotes equal except a*o*(1) term._{p} - Step 3 is to characterize the asymptotic expansion of $\sqrt{n}\{{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathit{\beta}}_{I\ast}(d)\}$. Substituting (23) into (19) leads to$$\sqrt{n}\{{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathit{\beta}}_{I\ast}(d)\}\approx {\mathbf{E}}_{II\ast}{(d,h)}^{-1}\frac{1}{\sqrt{n}}{G}_{n}({\mathit{\beta}}_{I\ast}(d);\omega ,h).$$(24)By using (21)–(22) and (24), we have$$\begin{array}{l}\sqrt{n}\{{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathit{\beta}}_{I\ast}(d)\}\approx \\ {\mathbf{E}}_{II\ast}{(d,h)}^{-1}\frac{1}{\sqrt{n}}{G}_{n}({\mathit{\beta}}_{I\ast}(d),{\mathit{\theta}}_{N\ast}({d}^{\prime};B(d,h));\omega ,h)\\ -{\mathbf{E}}_{II\ast}{(d,h)}^{-1}\sum _{{d}^{\prime}\in B(d,h)}{\mathbf{E}}_{IN\ast}({d}^{\prime};d,h)\sqrt{n}\{{\widehat{\mathit{\beta}}}_{N}({d}^{\prime})-{\mathit{\beta}}_{N\ast}({d}^{\prime})\}.\end{array}$$(25)
- Step 4 is to characterize the asymptotic expansion of $\sqrt{n}\{\widehat{\mathit{\beta}}(d)-{\mathit{\beta}}_{\ast}(d)\}$. Since
(_{N}*d*′) is a subvector of**(***d*′) and**(***d*′) is the solution of (5), it follows from the asymptotic results in (Liang and Zeger, 1986) that for any voxel*d*, $\sqrt{n}\{\widehat{\mathit{\beta}}(d)-{\mathit{\beta}}_{\ast}(d)\}$is asymptotically equivalent to$${\mathbf{F}}_{\ast}{(d)}^{-1}\frac{1}{\sqrt{n}}\sum _{i=1}^{n}{U}_{i}(\mathit{\beta}(d)){V}_{i}{({\mathit{\theta}}_{\ast}(d))}^{-1}{\mathbf{e}}_{i}(d;{\mathit{\beta}}_{\ast}(d)),$$(26)where**e**(_{i}*d*′;(*β**d*)) =**Y**(_{i}*d*′) –*μ*(_{i}(*β**d*)) for any*d*′ . - Step 5 is to characterize the asymptotic expansion of $\sqrt{n}\{{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathit{\beta}}_{I\ast}(d)\}$. By substituting (26) into (25), we have$$\sqrt{n}\{{\widehat{\mathit{\beta}}}_{I}(d,h)-{\mathit{\beta}}_{I\ast}(d)\}\approx {\mathbf{E}}_{II\ast}{(d,h)}^{-1}\frac{1}{\sqrt{n}}\sum _{i=1}^{n}\sum _{{d}^{\prime}\in B(d,h)}S({d}^{\prime};d,h),$$(27)where
*S*(*d*′;*d*,*h*) is given by$$\begin{array}{l}\omega (d,{d}^{\prime};h){U}_{iI}{({\mathit{\beta}}_{I\ast}(d),{\mathit{\beta}}_{N\ast}({d}^{\prime}))}^{T}{V}_{i}{({\mathit{\theta}}_{\ast}({d}^{\prime}))}^{-1}{\mathbf{e}}_{i}({d}^{\prime};{\mathit{\beta}}_{I\ast}(d),{\mathit{\beta}}_{N\ast}({d}^{\prime}))\\ -{\mathbf{E}}_{IN\ast}({d}^{\prime};d,h)[{\mathbf{0}}_{{q}_{N}{q}_{I}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\mathbf{I}}_{{q}_{N}}]{\mathbf{F}}_{\ast}{({d}^{\prime})}^{-1}{U}_{i}{({\mathit{\beta}}_{\ast}({d}^{\prime}))}^{T}{V}_{i}{({\mathit{\theta}}_{\ast}({d}^{\prime}))}^{-1}{\mathbf{e}}_{i}({d}^{\prime};{\mathit{\beta}}_{\ast}({d}^{\prime})).\end{array}$$Finally, it follows from (27) that the covariance matrix of $\sqrt{n}{\widehat{\mathit{\beta}}}_{I}(d,h)$, denoted by ${\mathrm{\sum}}_{n}(\sqrt{n}{\widehat{\mathit{\beta}}}_{I}(d,h))$), can be approximated by$${\mathbf{E}}_{II\ast}{(d,h)}^{-1}\frac{1}{n}\sum _{i=1}^{n}[\sum _{{d}^{\prime}\in B(d,h)}S({d}^{\prime};d,h)]{[\sum _{{d}^{\prime}\in B(d,h)}S({d}^{\prime};d,h)]}^{T}{\mathbf{E}}_{II\ast}{(d,h)}^{-1}.$$(28)

The PS procedure for MAGEE-A has five key steps: initialization, weights adaptation, estimation, stop checking, and inference.

- In the initialization step (i), we calculate the voxel-wise GEE estimate
**(***d*) at each voxel*d*. Then, we prefix a geometric series { ${h}_{s}={c}_{h}^{s}$:*s*= 1, …,*S*} of radii with*h*_{0}= 0, where*c*(1, 2), say_{h}*c*= 1.15. We suggest relatively small_{h}*c*, since small_{h}*c*prevents incorporating too many neighboring voxels in the beginning of PS, which improves the robustness of PS and the accuracy of parameter estimation. We then set_{h}*s*= 1 and*h*_{1}=*c*._{h} - In the weights adaptation step (ii), we compute
*D*(_{βI}*d*,*d*′;*h*_{s}_{−1}) and the adaptive weights*ω*(*d*,*d*′;*h*), which are defined as_{s}$$\begin{array}{l}{D}_{{\beta}_{I}}(d,{d}^{\prime};{h}_{s-1})={\mathrm{\Delta}}_{{\beta}_{I}}{(d,{d}^{\prime};{h}_{s-1})}^{T}{\mathrm{\sum}}_{n}{({\widehat{\mathit{\beta}}}_{I}(d;{h}_{s-1}))}^{-1}{\mathrm{\Delta}}_{{\beta}_{I}}(d,{d}^{\prime};{h}_{s-1}),\\ \omega (d,{d}^{\prime};{h}_{s})=\frac{{K}_{\mathit{loc}}({\left|\right|d-{d}^{\prime}\left|\right|}_{2}/{h}_{s}){K}_{st}({D}_{{\beta}_{I}}(d,{d}^{\prime};{h}_{s-1})/{C}_{n})}{{\sum}_{{d}^{\prime}\in B(d,{h}_{s})}{K}_{\mathit{loc}}({\left|\right|d-{d}^{\prime}\left|\right|}_{2}/{h}_{s}){K}_{st}({D}_{{\beta}_{I}}(d,{d}^{\prime};{h}_{s-1})/{C}_{n})},\end{array}$$(29)where Δ(_{βI}*d*,*d*′;*h*_{s}_{−1}) =(_{I}*d*,*h*_{s}_{−1}) −(_{I}*d*′,*h*_{s}_{−1}),*K*(_{loc}*u)*and*K*(_{st}*u*) are two nonnegative kernel functions with compact support,*C*is a number associated with_{n}*n*, and ||·||_{2}denotes the Euclidean norm of a vector (or a matrix). The weights*K*(||_{loc}*d*−*d*′||_{2}/*h*) give less weight to the voxel_{s}*d*′*B*(*d*,*h*), whose location is far from the voxel_{s}*d*. The weights*K*(_{st}*u*) downweight the voxels*d*′ with large differences between(_{I}*d*′,*h*_{s}_{−1}) and(_{I}*d*,*h*_{s}_{−1}). In practice, we set*K*(_{loc}*u*) = (1 −*u*)in order to increase the smoothness of the images of {_{+}(_{I}*d*,*h*) :*d*}, whereas we choose*K*(_{st}*u*) = exp(−*u*) in order to heavily penalize moderate to large differences. We choose*C*=_{n}*n*^{1/3}*χ*^{2}(*q*)_{I}^{0.8}based on our experiments, where χ^{2}(*q*)_{I}is the upper 1^{a}**–***a*-percentile of the*χ*^{2}(*q*) distribution._{I} - In the estimation step (iii), for the radius
*h*, we substitute_{s}*ω*(*d*,*d*′;*h*) into (7) to calculate_{s}(_{I}*d*,*h*) by using the Newton-Raphson algorithm at each voxel_{s}*d*. - In the stop checking step (iv), for
*s > S*_{0}, we compute$${D}_{I}(d;{h}_{s},{h}_{{S}_{0}})={\mathrm{\Delta}}_{{\beta}_{I}}{(d,{h}_{s},{h}_{{S}_{0}})}^{T}{\widehat{\mathrm{\sum}}}_{n}{({\widehat{\mathit{\beta}}}_{I}(d;{h}_{{S}_{0}}))}^{-1}{\mathrm{\Delta}}_{{\beta}_{I}}(d,{h}_{s},{h}_{{S}_{0}}),$$where Δ(_{βI}*d*,*h*,_{s}*h*_{S}_{0}) =(_{I}*d*,*h*_{S}_{0}) –(_{I}*d*,*h*). If_{s}*D*(_{I}*d; h*,_{s}*h*_{S}_{0}) ≥ (*s*), then we set(_{I}*d*,*h*) =_{S}(_{I}*d*,*h*_{s}_{−1}) and*s*=*S*. If*s*=*S*, we go to the inference step (v). If*s*≤*S*_{0}or*D*(_{I}*d; h*,_{s}*h*_{S}_{0}) ≤ for*S*− 1 ≥*s*>*S*_{0}, then we set*h*_{s}_{+1}=*c*, increase_{h}h_{s}*s*by 1 and continue with the weight adaptation step (ii). Throughout the paper, we set*S*_{0}= 3 and (*s*) =*χ*^{2}(*q*)_{I}^{0.8/(}^{s}^{−2)0.9}. - In the inference step (v), when
*s*=*S*, we report the final(_{I}*d*,*h*) and compute the_{S}*p*-values for*W*(_{β}*d*,*h*). In practice, we usually set the maximal step_{S}*S*to be relatively small, say 10, and thus each*B*(*d*,*h*_{10}) only contains a relatively small number of voxels compared to the whole volume. Throughout the paper, we have used the false discovery rate (FDR) method in (Benjamini and Yekutieli, 2001), since the test statistics obtained from PS satisfy the positive dependency condition. We may use other multiple comparison correction methods (e.g., random field theory (RFT), FDR, or permutation methods)(Benjamini and Yekutieli, 2001; Worsley et al., 2004; Chumbley et al., 2010; Salimi-Khorshidi et al., 2011).

^{1}This work was partially supported by NIH grants R01ES17240, MH091645, U54 EB005149, P30 HD03110, RR025747-01, P01CA142538-01, MH086633, AG033387, MH064065, HD053000, and MH070890. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

- Almli CR, Rivkin MJ, McKinstry RC. BDC Group. The nih mri study of normal brain development (objective-2): newborns, infants, toddlers, and preschoolers. IEEE Transactions on Medical Imaging. 2007;35:308–325. [PubMed]
- Arsigny V, Fillard P, Pennec X, Ayache N. Log-euclidean metrics for fast and simple caculus on diffusion tensors. Magnetic Resonance in Medicine. 2006;56:411–421. [PubMed]
- Basser PJ, Mattiello J, LeBihan D. Estimation of the effective self-diffusion tensor from the nmr spin echo. Journal of Magnetic Resonance Ser B. 1994;103:247–254. [PubMed]
- Bates D, Maechler M, Bolker B. lme4: Linear mixed-effects models using s4 classes. R package-42 2011
- Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics. 2001;29:1165–1188.
- Bernal-Rusiel J, Greve D, Reuter M, Fischl B, Sabuncu MR. Statistical analysis of longitudinal neuroimage data with linear mixed effects models. NeuroImage. 2013;66:249–260. [PMC free article] [PubMed]
- Chumbley J, Worsley KJ, Flandin G, Friston KJ. Topological fdr for neuroimaging. Neuroimage. 2010;49 (4):3057–3064. [PMC free article] [PubMed]
- Chung MK, Robbins S, Dalton KM, Davidson RJ, Alexander AL, Evans AC. Cortical thickness analysis in autism via heat kernel smoothing. NeuroImage. 2005;25:1256–1265. [PubMed]
- Diggle P, Heagerty P, Liang KY, Zeger S. Analysis of Longitudinal Data. 2. Oxford University Press; New York: 2002.
- Evans AC. BDC Group. The nih mri study of normal brain development. NeuroImage. 2006;30:184–202. [PubMed]
- Feng R, Zhou G, Zhang M, Zhang H. Analysis of twin data using sas. Biometrics. 2009;65:584–589. [PMC free article] [PubMed]
- Fitzmaurice G, Laird N, Ware J. Applied Longitudinal Analysis. John Wiley and Sons; New York: 2004.
- Gilmore JH, Schmitt JE, Knickmeyer RA, Smithm JK, Lin W, Styner M, Gerig G, Neale MC. Genetic and environmental contributions to neonatal brain structure: a twin study. Human Brain Mapping. 2010;31:1174–1182. [PMC free article] [PubMed]
- Goodlett CB, Fletcher PT, Gilmore JH, Gerig G. Group analysis of dti fiber tract statistics with application to neurodevelopment. NeuroImage. 2009;45:S133–S142. [PMC free article] [PubMed]
- Greven S, Crainiceanu C, Caffo B, Reich B. Longitudinal functional principal components analysis. Electronic Journal of Statistics. 2010;4:1022–1054. [PMC free article] [PubMed]
- Ibrahim JG, Molenberghs G. Missing data methods in longitudinal studies: a review. Test. 2009;18:1–43. [PMC free article] [PubMed]
- Jones DK, Symms MR, Cercignani M, Howard RJ. The effect of filter size on vbm analyses of dt-mri data. NeuroImage. 2005;26:546–554. [PubMed]
- Kim P, Leckman JF, Mayes L, Feldman R, Wang X, Swain JE. The plasticity of human maternal brain: longitudinal changes in brain anatomy during the early postpartum period. Behavioral Neuroscience. 2010;124:695–700. [PMC free article] [PubMed]
- Li Y, Gilmore J, Wang J, Styner M, Lin WL, Zhu H. Twin-marm: two-stage spatial adaptive analysis of twin neuroimaging data. IEEE Transactions on Medical Imaging. 2012;31:1100–1112. [PMC free article] [PubMed]
- Li Y, Zhu H, Shen D, Lin W, Gilmore JH, Ibrahim JG. Multiscale adaptive regression models for neuroimaging data. Journal of the Royal Statistical Society: Series B. 2011;73:559–578. [PMC free article] [PubMed]
- Liang K, Zeger S. Longitudinal data analysis using general linear models. Biometrika. 1986;73:13–22.
- Meltzer JA, Postman-Caucheteux W, McArdle JJ, Braun A. Strategies for longitudinal neuroimaging studies of overt language production. Neuroimage. 2009;47:745–755. [PMC free article] [PubMed]
- Mukherjee P, McKinstry RC. Diffusion tensor imaging and tractography of human brain development. Neuroimaging Clin N Am. 2006;16:19–43. [PubMed]
- Mukherjee PS, Qiu P. 3-d image denoising by local smoothing and nonparametric regression. Technometrics. 2011;53:196–208.
- Neil JJ, Shiran SI, McKinstry RC, Schefft GL, Snyder AZ, Almli CR, Akbudak E, Aronovitz JA, Miller JP, Lee BC, Conturo TE. Normal brain in human newborns: apparent diffusion coefficient and diffusion anisotropy measured by using diffusion tensor mr imaging. Radiology. 1998;209:57–66. [PubMed]
- Nichols T, Waldorp L. Accurate and computationally efficient analysis of longitudinal fmri data. 17th Annual Meeting of the Organization for Human Brain Mapping; Barcelona. 2010.
- Petersen R, Aisen P, Beckett L, Donohue M, Gamst A, Harvey D, Jack C, Jagust W, Shaw L, Toga A, Trojanowski J, Weiner M. Alzheimer’s disease neuroimaging initiative (ADNI): clinical characterization. Neurology. 2010;74:201–209. [PMC free article] [PubMed]
- Pinheiro J, Bates D, DebRoy S, Sarkar D. the R Development Core team. nlme: linear and nonlinear mixed effects models. R package version 3.1 2011
- Poline J, Mazoyer B. Enhanced detection in brain activation maps using a multifiltering approach. J Cereb Blood Flow Metab. 1994;14:639–642. [PubMed]
- Polzehl J, Spokoiny VG. Adaptive weights smoothing with applications to image restoration. J R Statist Soc B. 2000;62:335–354.
- Polzehl J, Spokoiny VG. Propagation-separation approach for local likelihood estimation. Probab Theory Relat Fields. 2006;135:335–362.
- Polzehl J, Voss HU, Tabelow K. Structural adaptive segmentation for statistical parametric mapping. Neuroimage. 2010;52:515–523. [PubMed]
- Pourahmadi M. Modeling covariance matrices: The glm and regularization perspectives. Statistical Science. 2011;26:369–387.
- Qiu P. Image Processing and Jump Regression Analysis. New York: John Wiley & Sons; 2005.
- Qiu P, Mukherjee PS. Edge structure preserving image denoising. Signal Processing. 2010;90:2851–2862.
- Ramsay JO, Silverman BW. Springer Series in Statistics. 2. Springer; New York: 2005. Functional data analysis.
- Ratcliffe S, Shults J. Geeqbox: A matlab toolbox for generalized estimating equations and quasi-least squares. Journal of Statistical Software. 2008;25:1–14.
- Salimi-Khorshidi G, Smith S, Nichols TE. Adjusting the effect of nonstationarity in cluster-based and tfce inference. Neuroimage. 2011;54 (3):2006–2019. [PubMed]
- Schneider J, Il’yasov K, Hennig J, Martin E. Fast quantitative diffusion-tensor imaging of cerebral white matter from the neonatal period to adolescence. Neuroradiology. 2004;46:258–266. [PubMed]
- Siegmund DO, Worsley KJ. Testing for a signal with unknown location and scale in a stationary gaussian random field. Ann Statist. 1995;23:608–639.
- Silver M, Montana G, Nichols TE, ADNI False positives in neuroimaging genetics using voxel-based morphometry data. NeuroImage. 2011;54:992–1000. [PMC free article] [PubMed]
- Skup M, Zhu H, Zhang H. Multiscale adaptive marginal analysis of longitudinal neuroimaging data with time-varying covariates. Biometrics. 2012 in press. [PMC free article] [PubMed]
- Skup M, Zhu HT, Wang Y, Giovanello KS, Lin JA, Shen DG, Shi F, Gao W, Lin W, Fan Y, Zhang HP, ADNI Sex differences in grey matter atrophy patterns among ad and amci patients: Results from adni. NeuroImage. 2011;56:890–906. [PMC free article] [PubMed]
- Tabelow K, Polzehl J, Spokoiny V, Voss HU. Diffusion tensor imaging: structural adaptive smoothing. NeuroImage. 2008;39:1763–1773. [PubMed]
- Tabelow K, Polzehl J, Voss HU, Spokoiny V. Analyzing fmri experiments with structural adaptive smoothing procedures. NeuroImage. 2006;33:55–62. [PubMed]
- Ball T, Breckel TP, Mutschler I, Aertsen A, Schulze-Bonhage A, Hennig J, Speck O. Variability of fmri-response patterns at different spatial observation scales. Human Brain Mapping. 2012;33:1155–1171. [PubMed]
- Wang N. Marginal nonparametric kernel regression accounting for within-subject correlation. Biometrika. 2003;90:43–52.
- Wang Y. Smoothing spline models with correlated random errors. Journal of the American Statistical Association. 1998;93 (441):341–348.
- Worsley KJ, Taylor JE, Tomaiuolo F, Lerch J. Unified univariate and multivariate random field theory. NeuroImage. 2004;23:189–195. [PubMed]
- Wu HL, Zhang JT. Nonparametric Regression Methods for Longitudinal Data Analysis. John Wiley & Sons, Inc; Hoboken, New Jersey: 2006.
- Yao F, Müller H-G, Wang J-L. Functional linear regression analysis for longitudinal data. Ann Statist. 2005;33(6):2873–2903.
- Yap P, Wu G, Zhu H, Lin W, Shen D. Timer: tensor image morphing for elastic registration. NeuroImage. 2009;47:549–563. [PMC free article] [PubMed]
- Zhang T, Davatzikos C. Odvba: optimally-discriminative voxel-based analysis. IEEE Transactions on Medical Imaging. 2011;30:1441–1454. [PMC free article] [PubMed]
- Zhao L, Boucher M, Rosa-Neto P, Evans A. Impact of scale space search on age- and gender-related changes in mri-based cortical morphometry. Human Brain Mapping. 2012 in press. [PubMed]
- Zhu H, Kong L, Li R, Styner M, Gerig G, Lin W, Gilmore JH. Fadtts: functional analysis of diffusion tensor tract statistics. NeuroImage. 2011;56:1412–1425. [PMC free article] [PubMed]
- Zhu HT, Li YM, Ibrahim JG, Lin WL, Shen DG. Marm: Multiscale adaptive regression models for neuroimaging data. Information Processing in Medical Imaging. 2009;5636:314–325. [PMC free article] [PubMed]
- Zhu HT, Zhang HP, Ibrahim JG, Peterson BG. Statistical analysis of diffusion tensors in diffusion-weighted magnetic resonance image data (with discussion) Journal of the American Statistical Association. 2007;102:1085–1102.
- Zimmerman D, Nunez-Anton V. Parametric modeling of growth curve data: an overview (with discussions) Test. 2001;10:1–73.

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. |