10.1 Repeated Measures and Longitudinal Data

The term repeated measures refers to experimental designs (or observational studies) in which each experimental unit (or subject) is measured at several points in time. The term longitudinal data is also used for this type of data.

Typical Design

Experimental units are randomly allocated to one of g treatments. A short time series is observed for each observation. An example in which there are 3 treatment groups with 3 units per treatment, and each unit is measured at four times is as follows:

In the example design just given, we have observed nine short time series, one for each experimental unit. The focus of the analysis would be to determine treatment differences, both in mean level and patterns across time. We may also want to characterize the overall time pattern (assuming that pattern doesn’t differ substantially for different treatments).

One approach to analyzing these data is the ANOVA model

Y = Treatment + Animal(Treatment) + Time + Treatment*Time + Error

• Animal is a “random” factor assumed to have mean 0 and an unknown constant variance. “Random” means that the animals are considered to be a random sample from a larger population of animals.
• The notation Animal(Treatment) specifies that animals were nested in treatment meaning that different animals were in the different treatment groups.
• The Treatment factor measures whether the mean response differs for the three treatments when we average over all animals and all times.
• The Time factor measures whether the mean response differs over time when we average over all animals and all treatments.
• The Time*Treatment interaction which is sensitive to whether the pattern across time depends upon the specific treatment used.
• The errors are assumed to be independently normally distributed with mean 0 and constant variance.

An analysis of variance for this example is given below. All factors are statistically significant. The Treatment factor is significant with p-value = .0001737. The Time factor is significant with p-value = .0001117 and the interaction between treatment and time is significant with p-value = .0207002. The significant interaction confirms that the pattern of change over time differs for the treatments.

The assumptions for this model lead to another assumption – that, in theory, the correlation is equal for all time gaps between observations. For instance, the correlation between data at times 1 and 2 is the same as the correlation between data at times 1 and 3, and is also the same as the correlation between times 1 and 4. This is called the compound symmetry assumption.

For the phlebitis example in which each animal is measured at four times, the assumed correlation matrix for the measurements at the four times has the structure

\(\left(\begin{array}{cccc}1 & \rho & \rho & \rho \\ \rho & 1 & \rho & \rho \\ \rho & \rho & 1 & \rho \\ \rho & \rho & \rho & 1\end{array}\right)\)

More flexible models for the correlation structure of the observations across time are available in R and SAS, but not in Minitab.

In R, the function gls within the nlme library can be used to specify several different structures for the correlations among measurements. For any specified assumption, maximum likelihood estimation is used to estimate the model parameters (including parameters in the correlation matrix). For each model, R calculates AIC and BIC statistics that can be used to compare models. An “adjusted” analysis of variance can also be calculated.

To use gls, you must first use the command

library(nlme)

Among the possibilities for assumptions about the correlations among observations across time are:

• Compound symmetry with the equal variances at the different time points (the assumption above)
• A completely unstructured set of correlations, meaning that the correlations are different for the various time gaps and have no ARIMA structure. This also allows for either equal or unequal variances at the different time points.
• An AR(1) structure, either with or without unequal variances
• A general ARIMA structure

Find a complete list of possibilities at Correlation Structure Classes.

The gls command is described here: Fit Linear Model Using Generalized Least Squares.

Code for the phlebitis example will be given (and explained) at the end of this lesson.

For the phlebitis data, we examined four possible models for the correlation structure of the repeated measurements:

• Compound symmetry
• Completely unstructured with possibly unequal variances
• AR(1) with equal variances
• AR(1) with unequal variances

The other components of each model were the same as before – the treatment factor, the time factor, the interaction between time and treatment, and a random animal factor.

A comparative summary of results, produced in R, is the following:

The AR(1) model for correlations among repeated measures gives the lowest AIC and BIC statistics, although not by much. The original compound symmetry model is a close second. The completely unstructured model (for correlations) is the worst based on AIC and BIC. This model has the maximum log likelihood value (a good thing), but the gain in log likelihood was outweighed by the increased number of parameters that had to be estimated.

The table above also includes likelihood ratio tests that compare the models sequentially, based on the order of models given in the command that created the table. The L. Ratio column = −2 *Difference in log likelihood values. This statistic has an approximate chi-square distribution with df = difference in DF values. The DF value given for a model is the total number of parameters that had to be estimated for all effects in the model. The chi-square distribution is used to determine a p-value for a test of whether the additional parameters = 0 or not. Failure to reject means that the models are not significantly different. That’s the case here.

The following table compares analysis of variance results for the best two models of the models that we tried. There’s not much difference! For this example, we get essentially the same results for the important factors of interest regardless of the model used for correlations across time.

We might be interested in whether the pattern across time is linear, quadratic, and so on, and whether treatments differ with regard to linear trend, quadratic trend, and so on. This can be done using orthogonal polynomial components for Time. The phrase orthogonal polynomial refers to a coding of time (or any other variable) such that polynomial components (linear, quadratic, cubic, and so on) are statistically independent of one another.

It might help to look back at the graph above in this lesson. We’re trying to characterize the trend pattern across time and compare the treatments with regard to the pattern(s).

In our example, there are four time points so we can fit a cubic polynomial (at most) to the mean responses at the four times. This is an extension of “two points determine a line.” In this case, four points determine a cubic. We then examine the significance of the various polynomial components in order to describe the trend across time, and we can also look at the interactions between the polynomial components and treatments.

It’s important to note that the use of polynomial trends does not change the overall results of F tests for the factors Time and Time*Treatment. We are simply considering the four time points in a more structured way.

The following output was produced in R by specifying an orthogonal polynomial of degree 3 (a cubic) for time and using an AR(1) specification for the correlations among repeated measures.

• The first two lines under the headings give comparison of treatment 2 to treatment 1 and treatment 3 to treatment 1. We can tell this by the numbers at the end of the coefficient descriptions and the fact that the there isn’t a line for treatment 1.
• The next three lines give significance results for the linear, quadratic and cubic components (in that order) of a cubic polynomial for the time trend for treatment 1 only.
• The final six lines give significance results for interactions between treatment and the polynomial components. The following two lines have to do with the linear component:

These results indicate that the linear trend (slope) differs for the treatments. More specifically the first line indicates that the slope differs for treatments 1 and 2. The second line indicates that the linear slope differs for the treatments 1 and 3.

The remaining interaction terms are not statistically significant.

In maximum likelihood estimation, the AIC and BIC values change due to the use of an orthogonal polynomial for the time trend. This might lead you to think that the fit has been improved. That’s not the case. The overall fit and analysis of variance is still the same, but the scaling of the likelihood is affected by the scaling of the coding for the factors.

Here’s what happened when we tried to compare the fits of the original model with no structure for time trend and the model with an orthogonal polynomial for time:

The data are in a “csv” file named phlebitis.csv. The dataset includes variables named Treatment, Time, Animal, and Y.

The first line reads the data, indicating that the first row gives text variable names, and a “,” separates data value. The second line makes the variable names available to R commands and the last line removes them when finished working with phlebitisdata. The next two lines have explanations given within the command lines.

The interaction.plot command created the Phlebitis Interaction Plot. The “lty” parameter assigned line types to the 3 treatment. If we have four treatments, for example, switch that to lty = c(1:4). The lwd = 2 parameter is optional. It increases the line widths. We did this to be sure that the line patterns (dashes, etc, could be seen).

The command that starts “nestinginfo” defines the nesting property of animals within treatments. A “new” dataset is created that actually is the same as the original, but information about nesting is added. (Nesting means different subjects are used within differing levels of another factor.)

The gls commands define a basic model involving Treatment, Time and their interaction, along with the nesting and a correlation model. The “factor” specification orders R to treat Time and Treatment as categorical factors rather than regression type variables.

Two commands used to examine results are summary and anova.

• The command summary(object name), where object name is a name attached to a gls analysis will give AIC and BIC values along with coefficient estimates. In the code this was used only for the last fit (with polynomial trend), but it could have been used for any fit.
• The command anova(object name) will give an analysis of variance for the gls with that object name.
• The command anova(list of object names) – for example anova(fit.compsym, fit.nostruct, fit.ar1, fit.ar1het) will compare AIC and BIC values for the various models.

Alternative references for this lesson are

PSU World Campus STAT 502 Lesson 10

UCLA-idre: Repeated Measures Analysis with R