Covariance Structure Selection in Longitudinal Models

Introduction

Longitudinal clinical trials collect repeated measurements from the same participant over time. A patient may have a baseline measurement followed by assessments at Week 2, Week 4, Week 8, Week 12, and later visits.

These observations are not independent.

Measurements from the same patient tend to be correlated because they share the patient's underlying biological characteristics, treatment response, disease trajectory, and measurement environment.

A longitudinal mixed model therefore needs to describe not only the expected mean trajectory but also the covariance among repeated observations within a subject.

Covariance structure selection is therefore not a cosmetic modeling decision. It determines how the model represents within-subject dependence and can affect standard errors, degrees of freedom, confidence intervals, test statistics, and conclusions about treatment effects.

The General Longitudinal Mixed Model

A common linear mixed model can be written as:

where:

• \(\mathbf{y}_i\) is the vector of repeated observations for subject \(i\).
• \(\mathbf{X}_i\boldsymbol{\beta}\) represents the fixed-effects mean structure.
• \(\mathbf{Z}_i\mathbf{b}_i\) represents subject-specific random effects, if included.
• \(\boldsymbol{\epsilon}_i\) represents residual variation.

The marginal covariance matrix of the repeated observations is:

The matrix \(\mathbf{R}_i\) describes the residual covariance structure after accounting for the random effects.

Why Covariance Structure Matters

Suppose a patient has measurements:

The four observations provide four pieces of information, but they do not represent four independent patients.

If the correlation between measurements is ignored, the model can incorrectly estimate the amount of independent information contained in the data.

Covariance structure selection attempts to answer questions such as:

• Are all visits equally correlated?
• Does correlation decline as the time between visits increases?
• Does residual variability increase over time?
• Should each pair of visits have its own covariance?
• Are observations closer in time more strongly correlated?

The Covariance Matrix

For four repeated measurements, the covariance matrix might be written as:

The diagonal elements are variances.

The off-diagonal elements are covariances.

A correlation matrix instead standardizes the covariance terms:

Different covariance structures impose different patterns on these variances and correlations.

The Most Common Covariance Structures

The most commonly encountered structures in clinical longitudinal analyses include:

Here \(m\) denotes the number of repeated measurements.

Compound Symmetry

Compound symmetry, or CS, assumes that every measurement has the same variance and every pair of measurements has the same correlation.

For four visits:

There are only two covariance parameters:

• \(\sigma^2\): common variance
• \(\rho\): common correlation

When CS Makes Sense

CS can be reasonable when measurements are approximately equally variable and the correlation does not appear to depend strongly on how far apart the visits are.

It is particularly attractive when the number of repeated measurements is large relative to the available sample size because it uses very few covariance parameters.

Autoregressive AR(1)

An AR(1) covariance structure assumes that correlation declines as the lag between measurements increases.

For equally spaced visits:

For four visits, the correlation matrix is:

If \(\rho=0.70\), for example:

Thus adjacent observations are more strongly correlated than observations separated by several visits.

When AR(1) Makes Sense

• Measurements are approximately equally spaced.
• Correlation decreases with increasing lag.
• Nearby visits are expected to be more similar than distant visits.

Heterogeneous Compound Symmetry

A heterogeneous compound symmetry structure, often denoted CSH, relaxes the equal-variance assumption of CS.

For four visits:

The correlations remain common, but the variances can differ by visit.

This can be useful when variability increases as disease progression or treatment exposure increases.

Heterogeneous AR(1)

Heterogeneous AR(1) allows both the variances and the correlations to vary according to an AR(1)-type pattern.

The covariance between visits \(j\) and \(k\) can be represented as:

This structure is substantially more flexible than ordinary AR(1), but it also requires substantially more covariance parameters.

Toeplitz Covariance

A Toeplitz structure assumes that covariance depends on the lag between measurements, but unlike AR(1), it does not require a geometric decay pattern.

For four visits:

The first off-diagonal represents lag 1, the second lag 2, and so on.

Toeplitz therefore provides more flexibility than AR(1) while still imposing a useful lag-based pattern.

Unstructured Covariance

The unstructured covariance structure, abbreviated UN, estimates every variance and covariance separately.

For four visits:

For \(m\) visits, UN requires:

covariance parameters.

For example:

Spatial Covariance Structures

Longitudinal data are sometimes collected at irregularly spaced times.

For example, suppose measurements occur at:

• Day 0
• Day 7
• Day 21
• Day 60
• Day 120

The gaps are clearly not equal.

A discrete AR(1) structure based simply on visit number treats the difference between visits 1 and 2 similarly to the difference between visits 4 and 5, even though the actual time intervals differ dramatically.

A continuous-time spatial covariance model can instead use actual time distance.

A common exponential correlation model is:

where \(t_j\) and \(t_k\) are the actual measurement times.

Covariance Structure vs. Random Effects

One of the most important concepts in longitudinal modeling is distinguishing the covariance structure from the random-effects structure.

Consider a random-intercept model:

where:

and:

Two observations from the same subject have covariance:

This induces a compound-symmetry-like marginal correlation pattern.

Adding a random slope produces a richer covariance pattern:

The covariance between observations can then depend on their time values.

Marginal and Conditional Perspectives

Mixed models can be interpreted conditionally on subject-specific random effects or marginally over those random effects.

This distinction matters when thinking about covariance structure.

The conditional residual covariance is represented by \(\mathbf{R}\), whereas the marginal covariance incorporates the random-effects contribution:

Two models that use different random-effects and residual structures can sometimes produce similar marginal covariance patterns.

Therefore, covariance selection should be considered as part of the complete model specification rather than as an isolated technical step.

How Do We Select a Covariance Structure?

There is no single universal rule.

A practical selection process combines:

• Scientific knowledge
• Study design
• Visit spacing
• Observed variance patterns
• Observed correlation patterns
• Model convergence
• AIC or BIC
• Likelihood-based comparisons when appropriate
• Interpretability
• Stability of treatment-effect estimates

Maximum Likelihood vs. Restricted Maximum Likelihood

This is one of the most important technical issues in covariance structure selection.

Maximum Likelihood

Maximum likelihood, or ML, estimates both fixed effects and variance parameters using the full likelihood.

ML is appropriate when comparing models with different fixed-effects structures because the fixed-effects contribution is part of the likelihood comparison.

Restricted Maximum Likelihood

Restricted maximum likelihood, or REML, estimates variance components after accounting for the fixed-effects estimation.

REML generally provides less biased variance-component estimates in finite samples.

However, REML likelihoods from models with different fixed-effects structures are not directly comparable for ordinary likelihood-ratio model selection.

AIC

The Akaike Information Criterion is:

where:

• \(\ell\) is the maximized log-likelihood.
• \(k\) is the number of estimated parameters.

Lower AIC is preferred.

AIC rewards better likelihood fit but penalizes model complexity.

BIC

The Bayesian Information Criterion is:

where \(n\) is the effective sample size used in the likelihood.

BIC applies a stronger complexity penalty than AIC when the sample size is greater than approximately \(e^2\).

Consequently, BIC will often favor a simpler covariance structure than AIC.

Do Not Select Covariance Structure by AIC Alone

Suppose the following results are obtained:

AIC favors UN, while BIC favors CSH.

This is not a contradiction.

It means the criteria place different weights on model flexibility.

Likelihood Ratio Tests

For nested models, a likelihood ratio statistic is:

Under standard regularity conditions, this can be compared with a chi-square distribution with degrees of freedom equal to the difference in the number of parameters.

However, covariance parameters can lie on the boundary of their parameter space. This means the usual chi-square approximation can require special attention in some variance-component testing problems.

Empirical Clues from the Data

Before fitting highly parameterized structures, examine the observed data.

Useful exploratory quantities include:

• Variance by visit
• Standard deviation by visit
• Pairwise correlations by visit
• Mean trajectory
• Number of observations at each visit
• Actual elapsed time between measurements

Suppose the observed standard deviations look approximately like:

Equal-variance structures such as CS or AR(1) may be reasonable candidates.

Now suppose:

A heterogeneous variance structure deserves consideration.

Correlation Patterns Can Be Informative

Suppose the estimated correlations appear approximately:

This pattern is consistent with decreasing correlation as time separation increases.

An AR-type covariance structure may therefore be scientifically reasonable.

A Critical Issue: Visit Number vs. Actual Time

Consider two patients.

Patient A is measured at:

Patient B is measured at:

Both may have four nominal visits, but the actual time separations are very different.

A standard AR(1) based on visit index treats the observation sequence according to lags of 1, 2, and 3 visits.

A continuous-time covariance model instead uses actual time separation.

Missing Data and Covariance Structure

Longitudinal clinical trials frequently have incomplete repeated measurements.

For example:

Mixed models can accommodate unbalanced longitudinal data under appropriate missing-data assumptions without requiring every subject to have observations at every visit.

However, missingness can make covariance estimation more difficult, particularly for highly parameterized structures such as UN.

Covariance Structure and MAR

Standard likelihood-based longitudinal mixed models are often interpreted under a missing-at-random framework.

The likelihood uses the observed measurements and does not require simple mean imputation for missing outcomes.

But the validity of the analysis depends on the missingness assumptions and the correctly specified model.

Changing the covariance structure does not automatically solve informative dropout.

A Worked Clinical Example

Suppose a randomized clinical trial compares an active treatment with placebo.

The primary continuous endpoint is measured at:

• Baseline
• Week 2
• Week 4
• Week 8
• Week 12
• Week 24

The planned longitudinal model contains:

• Treatment
• Visit
• Treatment-by-visit interaction
• Baseline covariate
• Subject-level repeated observations

The scientific question is whether the treatment groups differ in mean change over time.

Step 1: Define the Mean Model

A categorical-visit model might be:

The treatment-by-visit interaction allows treatment differences to vary across visits.

Step 2: Fit Candidate Covariance Structures

Suppose four candidate structures are considered:

Step 3: Compare the Models

Suppose the models produce:

UN has the best raw likelihood, but its large number of covariance parameters produces a much larger BIC.

CSH has the smallest AIC in this hypothetical comparison.

The selection should therefore not be reduced to the single statement "UN fits best."

Step 4: Examine Covariance Estimates

Suppose CSH produces estimated standard deviations:

These estimates suggest increasing variability over time.

If that pattern is clinically plausible, CSH may be preferable to CS.

Step 5: Examine the Treatment Effect

The final model should not be selected solely because its AIC is smallest.

Compare the treatment effect estimates across reasonable covariance structures.

The treatment estimate is stable across structures.

That provides reassurance that the scientific conclusion is not being driven by a fragile covariance assumption.

What If the Treatment Effect Changes Dramatically?

Suppose instead the results were:

The UN model may be too flexible for the available information.

Its treatment estimate is similar, but its uncertainty is much larger.

This is one reason a highly flexible covariance model should not automatically be selected simply because it minimizes the unpenalized likelihood.

Convergence Is Not Optional

A covariance model that does not converge properly should not be selected.

Warning signs include:

• Failure to converge
• Non-positive-definite covariance matrix
• Variance estimates approaching zero
• Correlation estimates approaching \(\pm1\)
• Extremely large standard errors
• Huge condition numbers
• Different answers from different starting values
• Optimization warnings

Boundary Estimates

Suppose an estimated variance is:

or a correlation is:

This can indicate that the data do not provide enough information to estimate the covariance parameters reliably.

A simpler structure may be more appropriate.

Positive Definiteness

A valid covariance matrix must be positive semidefinite.

For a covariance matrix \(\mathbf{R}\), we require:

for every vector \(\mathbf{a}\).

This requirement is why covariance structures cannot simply be constructed by independently estimating arbitrary pairwise correlations.

The model-fitting software imposes the mathematical constraints required for a valid covariance matrix.

R Example Using lme4

The lme4 package is widely used for linear mixed models, particularly for random-effects structures.

library(lme4)

fit_ri <- lmer(
  outcome ~ treatment * visit + baseline +
    (1 | subject),
  data = dat,
  REML = TRUE
)

summary(fit_ri)

This model contains a random intercept but does not provide the same flexible residual covariance structures available in specialized mixed-model software.

R Example Using nlme: Compound Symmetry

library(nlme)

fit_cs <- lme(
  outcome ~ treatment * visit + baseline,
  random = ~ 1 | subject,
  correlation = corCompSymm(
    form = ~ visit | subject
  ),
  data = dat,
  method = "REML"
)

summary(fit_cs)

R Example Using AR(1)

fit_ar1 <- lme(
  outcome ~ treatment * visit + baseline,
  random = ~ 1 | subject,
  correlation = corAR1(
    form = ~ visit | subject
  ),
  data = dat,
  method = "REML"
)

summary(fit_ar1)

The AR(1) model assumes correlation depends on the lag in the ordering of the measurement occasions.

R Example with Continuous-Time Correlation

For irregularly spaced measurements, a continuous-time correlation structure may be more appropriate.

fit_ct_ar1 <- lme(
  outcome ~ treatment * visit + baseline,
  random = ~ 1 | subject,
  correlation = corCAR1(
    form = ~ time_days | subject
  ),
  data = dat,
  method = "REML"
)

summary(fit_ct_ar1)

Here the correlation is based on the actual continuous measurement time rather than simply the visit number.

R Example with Heterogeneous Variance

A variance function can allow the residual variance to differ across visits.

fit_csh <- lme(
  outcome ~ treatment * visit + baseline,
  random = ~ 1 | subject,
  correlation = corCompSymm(
    form = ~ visit | subject
  ),
  weights = varIdent(
    form = ~ 1 | visit
  ),
  data = dat,
  method = "REML"
)

summary(fit_csh)

The correlation is constrained to a common value while the residual variances can differ across visits.

Comparing AIC in R

AIC(fit_cs, fit_ar1, fit_csh)
BIC(fit_cs, fit_ar1, fit_csh)

Lower values indicate a preferred model according to the corresponding information criterion.

Using ML When Comparing Fixed Effects

Suppose you want to compare:

outcome ~ treatment + visit + baseline

with:

outcome ~ treatment * visit + baseline

The models have different fixed effects.

They should generally be compared using ML rather than REML.

fit_main <- lme(
  outcome ~ treatment + visit + baseline,
  random = ~ 1 | subject,
  correlation = corAR1(form = ~ visit | subject),
  data = dat,
  method = "ML"
)

fit_interaction <- lme(
  outcome ~ treatment * visit + baseline,
  random = ~ 1 | subject,
  correlation = corAR1(form = ~ visit | subject),
  data = dat,
  method = "ML"
)

anova(fit_main, fit_interaction)

SAS Example: Compound Symmetry

In SAS PROC MIXED, compound symmetry can be specified with TYPE=CS.

proc mixed data=dat method=reml;
    class subject treatment visit;

    model outcome =
        treatment
        visit
        treatment*visit
        baseline
        / ddfm=kr;

    repeated visit /
        subject=subject
        type=cs;

run;

The Kenward-Roger degrees-of-freedom method is commonly used in clinical longitudinal analyses when appropriate.

SAS Example: AR(1)

proc mixed data=dat method=reml;
    class subject treatment visit;

    model outcome =
        treatment
        visit
        treatment*visit
        baseline
        / ddfm=kr;

    repeated visit /
        subject=subject
        type=ar(1);

run;

SAS Example: Unstructured Covariance

proc mixed data=dat method=reml;
    class subject treatment visit;

    model outcome =
        treatment
        visit
        treatment*visit
        baseline
        / ddfm=kr;

    repeated visit /
        subject=subject
        type=un;

run;

For \(m\) visits, the UN structure estimates:

covariance parameters.

SAS Example: Heterogeneous Compound Symmetry

proc mixed data=dat method=reml;
    class subject treatment visit;

    model outcome =
        treatment
        visit
        treatment*visit
        baseline
        / ddfm=kr;

    repeated visit /
        subject=subject
        type=csh;

run;

SAS Example: Compare Candidate Structures

A practical workflow is to fit candidate structures under the same fixed-effects model and collect fit statistics.

ods output FitStatistics=fit_cs;

proc mixed data=dat method=reml;
    class subject treatment visit;

    model outcome =
        treatment
        visit
        treatment*visit
        baseline
        / ddfm=kr;

    repeated visit /
        subject=subject
        type=cs;
run;

ods output FitStatistics=fit_ar1;

proc mixed data=dat method=reml;
    class subject treatment visit;

    model outcome =
        treatment
        visit
        treatment*visit
        baseline
        / ddfm=kr;

    repeated visit /
        subject=subject
        type=ar(1);
run;

The resulting fit statistics can then be reviewed alongside convergence, covariance estimates, and treatment-effect results.

Why Unstructured Is Not Always the Gold Standard

UN is often described as the most flexible covariance structure.

That does not mean it is automatically the best structure.

Suppose a study has 10 visits.

UN requires:

covariance parameters.

That can be a substantial burden, particularly if the study has modest sample size or considerable dropout.

A simpler structure such as AR(1), CSH, or Toeplitz may provide a better balance between flexibility and stability.

The Parameter-to-Information Problem

A covariance structure becomes increasingly difficult to estimate as the number of parameters grows relative to the number of subjects contributing information.

This is especially important because covariance parameters are estimated from within-subject repeated observations.

A trial with 1,000 observations may still have only 100 subjects.

The effective information for subject-level covariance estimation is therefore very different from treating all 1,000 observations as independent.

Small Samples and Overparameterization

Suppose there are only 60 subjects and 8 planned visits.

UN requires:

covariance parameters.

This may be ambitious depending on the amount of missingness, treatment allocation, and distribution of observations across visits.

A lower-dimensional structure may be much more stable.

Covariance Structure Selection Is Not the Same as Mean Model Selection

Suppose the treatment-by-visit interaction is significant under CS but not under UN.

This does not necessarily mean that one model is "correct" and the other is "wrong."

The estimated standard errors depend on the covariance model.

The appropriate question is whether the covariance model is adequately supported by the data and scientific context.

Pre-Specification in Clinical Trials

Covariance structure selection should ideally be addressed before the unblinded primary analysis.

The protocol or statistical analysis plan may specify:

• Primary mean model
• Primary covariance structure
• Alternative covariance structures
• Model-selection criteria
• Rules for handling convergence failure
• Degrees-of-freedom method
• Estimation method
• Missing-data assumptions

For example:

Should the Best-Fitting Structure Always Be the Primary Structure?

Not necessarily.

Suppose two structures have nearly identical AIC values:

The AIC difference is only:

Such a small difference does not provide strong evidence that CSH is meaningfully better.

If AR(1) has a clear scientific rationale and substantially simpler interpretation, it may be preferable.

AIC Differences

A useful quantity is:

The model with the smallest AIC has \(\Delta=0\).

Models with small differences in AIC have similar empirical support relative to the candidate set.

This is more informative than simply ranking models 1, 2, 3, and 4.

AICc for Smaller Samples

When sample size is relatively small relative to the number of estimated parameters, a finite-sample correction can be useful.

A common form is:

where \(k\) is the number of parameters and \(n\) is the relevant sample size.

The correction becomes increasingly important as \(k\) approaches \(n\).

Common Covariance Selection Scenarios

Common Mistakes

1. Assuming CS is always the default answer. CS is simple and useful, but it may be unrealistic when correlation changes with time.
2. Using AR(1) with highly irregular visit spacing without considering actual time. The discrete AR(1) model is based on visit lag rather than elapsed time.
3. Choosing UN automatically because it makes fewer assumptions. UN may require many parameters and can be unstable.
4. Comparing REML likelihoods across different fixed-effects models. Use ML when comparing fixed-effects structures.
5. Choosing the structure solely by p-value. Covariance selection should not be driven by whether the treatment effect becomes statistically significant.
6. Ignoring convergence warnings. An unstable covariance model is not a valid preferred model simply because its AIC is low.
7. Ignoring variance heterogeneity. Residual variability may increase or decrease substantially over time.
8. Ignoring dropout. Late-visit covariance parameters can become difficult to estimate when many subjects discontinue.
9. Confusing random effects with residual covariance. Both contribute to within-subject dependence but play different roles in the model.
10. Changing the covariance structure after seeing treatment results. Post hoc selection based on the desired inferential result can undermine the credibility of the analysis.

A Practical Model Selection Strategy

A robust approach is to start with a scientifically reasonable set of candidate structures rather than fitting every possible structure.

What If Candidate Models Give Similar Results?

This is often a good outcome.

Suppose CS, AR(1), and CSH all produce similar treatment estimates and similar confidence intervals.

Then the scientific conclusion is relatively robust to covariance specification.

In such a situation, selection can reasonably emphasize:

• Scientific interpretability
• Parsimony
• Convergence
• Pre-specification
• Consistency with prior studies

What If Candidate Models Give Different Results?

This deserves more attention.

Large differences in standard errors or treatment conclusions may indicate that the data do not strongly identify the within-subject covariance pattern.

Investigate:

• Whether the most flexible model is overparameterized
• Whether there is substantial dropout
• Whether visit spacing is irregular
• Whether variances differ substantially by visit
• Whether the mean model is correctly specified
• Whether the random-effects structure is appropriate

Sensitivity analyses may be particularly valuable in this situation.

Covariance Structure and Degrees of Freedom

In finite samples, uncertainty in estimated covariance parameters can affect the degrees of freedom used for hypothesis tests.

Methods such as the Kenward-Roger approximation attempt to improve small-sample inference by accounting for covariance-parameter estimation.

Another commonly used approach is the Satterthwaite approximation.

Covariance Selection in MMRM

The mixed model for repeated measures, or MMRM, is widely used for longitudinal continuous endpoints in randomized clinical trials.

A typical MMRM includes:

• Treatment
• Visit
• Treatment-by-visit interaction
• Baseline covariate
• Repeated observations within subject

Covariance structure selection is therefore a central part of MMRM specification.

UN is commonly considered because it makes relatively few assumptions about the covariance matrix, but alternatives such as CS, AR(1), CSH, or Toeplitz may be more stable depending on the study.

Why MMRM Makes Covariance Selection Especially Important

In an MMRM, the treatment-by-visit interaction may be estimated using within-subject longitudinal information.

Therefore, the covariance matrix directly affects the precision with which these visit-specific treatment differences are estimated.

A poorly specified covariance structure can lead to inefficient or unstable inference.

Baseline Adjustment and Covariance Structure

Baseline adjustment is frequently included in longitudinal models.

For example:

The covariance structure describes dependence among post-baseline outcomes after accounting for the fixed effects.

It is therefore important to distinguish:

• Adjustment for baseline in the mean model
• Within-subject covariance among repeated outcomes

These address different statistical questions.

Covariance Selection and Sensitivity Analysis

A strong longitudinal analysis can evaluate whether conclusions are robust under alternative reasonable covariance structures.

For example:

If the treatment conclusion remains similar across these analyses, confidence in the robustness of the result increases.

Do Not Use Covariance Selection to "Find Significance"

Suppose an analyst fits six covariance structures and reports the one that produces the smallest treatment p-value.

This is poor statistical practice.

The covariance structure should be selected according to a prespecified or scientifically defensible procedure, not according to which model produces the most favorable treatment conclusion.

A Compact Decision Table

A Complete Covariance Selection Checklist

• Confirm the measurement schedule.
• Determine whether visit intervals are approximately equal.
• Examine the number of subjects relative to the number of visits.
• Assess the amount and pattern of dropout.
• Specify the fixed-effects mean model.
• Identify plausible covariance structures.
• Fit models under the same data and fixed-effects specification.
• Use REML for covariance comparison when fixed effects are unchanged.
• Use ML when comparing different fixed-effects structures.
• Review AIC.
• Review BIC.
• Consider likelihood-ratio tests for appropriate nested models.
• Check convergence.
• Check positive definiteness.
• Inspect covariance parameter estimates.
• Compare treatment-effect estimates and standard errors.
• Consider sensitivity analyses.
• Document the final selection rationale.

The Most Important Concept

The most important point is that covariance structure selection is a balance between realism and estimability.

A very simple model such as CS may impose unrealistic assumptions.

A very flexible model such as UN may impose too few structural assumptions but require too many parameters.

AR(1), CSH, Toeplitz, heterogeneous AR-type structures, and continuous-time models occupy different positions along this spectrum.

The best covariance structure is therefore not necessarily the one with the largest number of parameters or the smallest raw likelihood.

Summary

Covariance structure selection is a fundamental component of longitudinal mixed-effects modeling.

The major structures have different assumptions:

The selection process should consider:

• Scientific plausibility
• Visit spacing
• Observed variance patterns
• Observed correlation patterns
• AIC and BIC
• Likelihood comparisons
• Convergence
• Positive definiteness
• Parameter stability
• Treatment-effect robustness

For clinical longitudinal studies, the final covariance structure should be documented clearly in the statistical analysis plan and supported by a reproducible model-selection strategy.

References

Laird, N.M. & Ware, J.H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963–974.

Verbeke, G. & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer.

Fitzmaurice, G.M., Laird, N.M. & Ware, J.H. (2011). Applied Longitudinal Analysis. Wiley.

Little, R.J.A., et al. (2012). The prevention and treatment of missing data in clinical trials. New England Journal of Medicine, 367, 1355–1360.

Kenward, M.G. & Roger, J.H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53, 983–997.

Wolfinger, R. (1993). Covariance structure selection in general mixed models. Communications in Statistics — Simulation and Computation, 22, 1079–1106.

Pinheiro, J.C. & Bates, D.M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.

Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19, 716–723.

Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6, 461–464.