### Statistical Computing Seminars Repeated Measures Analysis with R

There are a number of situations that can arise when the analysis includes between groups effects as well as within subject effects. We start by showing 4 example analyses using measurements of depression over 3 time points broken down by 2 treatment groups. In the first example we see that the two groups differ in depression but neither group changes over time. In the second example the two groups grow in depression but at the same rate over time.  In the third example, the two groups start off being quite different in depression but end up being rather close in depression. The fourth example shows the groups starting off at the same level of depression, and one group group increases over time whereas the other group decreases over time.

Note that in the interest of making learning the concepts easier we have taken the liberty of using only a very small portion of the output that R provides and we have inserted the graphs as needed to facilitate understanding the concepts. The code needed to actually create the graphs in R has been included.

#### Demo Analysis #1

The between groups test indicates that the variable group is significant, consequently in the graph we see that the lines for the two groups are rather far apart. The within subject test indicate that there is not a significant time effect, in other words, the groups do not change in depression over time. In the graph we see that the groups have lines that are flat, i.e. the slopes of the lines are approximately equal to zero. Also, since the lines are parallel, we are not surprised that the interaction between time and group is not significant.


## Convert variables to factor
demo1 <- within(demo1, {
group <- factor(group)
time <- factor(time)
id <- factor(id)
})

par(cex = .6)

with(demo1, interaction.plot(time, group, pulse,
ylim = c(5, 20), lty= c(1, 12), lwd = 3,
ylab = "mean of pulse", xlab = "time", trace.label = "group"))

demo1.aov <- aov(pulse ~ group * time + Error(id), data = demo1)
summary(demo1.aov)

Error: id
Df Sum Sq Mean Sq F value    Pr(>F)
group      1 155.04 155.042    3721 1.305e-09 ***
Residuals  6   0.25   0.042
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
Df  Sum Sq  Mean Sq F value Pr(>F)
time        2 0.08333 0.041667       1 0.3966
group:time  2 0.08333 0.041667       1 0.3966
Residuals  12 0.50000 0.041667


#### Demo Analysis #2

The between groups test indicates that the variable group is not significant, consequently in the graph we see that the lines for the two groups are rather close together. The within subject test indicate that there is a significant time effect, in other words, the groups do change in depression over time. In the graph we see that the groups have lines that increase over time. Again, the lines are parallel consistent with the finding that the interaction is not significant.


## Convert variables to factor
demo2 <- within(demo2, {
group <- factor(group)
time <- factor(time)
id <- factor(id)
})
par(cex = .6)

with(demo2, interaction.plot(time, group, pulse,
ylim = c(10, 40), lty = c(1, 12), lwd = 3,
ylab = "mean of pulse", xlab = "time", trace.label = "group"))

demo2.aov <- aov(pulse ~ group * time + Error(id), data = demo2)
summary(demo2.aov)

Error: id
Df  Sum Sq Mean Sq F value Pr(>F)
group      1  15.042  15.042  0.8363 0.3957
Residuals  6 107.917  17.986

Error: Within
Df Sum Sq Mean Sq F value    Pr(>F)
time        2 978.25  489.12 53.6845 1.032e-06 ***
group:time  2   1.08    0.54  0.0595    0.9426
Residuals  12 109.33    9.11
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#### Demo Analysis #3

The between groups test indicates that the variable group is significant, consequently in the graph we see that the lines for the two groups are rather far apart. The within subject test indicate that there is a significant time effect, in other words, the groups do change over time, both groups are getting less depressed over time. Moreover, the interaction of time and group is significant which means that the groups are changing over time but are changing in different ways, which means that in the graph the lines will not be parallel. In the graph we see that the groups have non-parallel lines that decrease over time and are getting progressively closer together over time.


## Convert variables to factor
demo3 <- within(demo3, {
group <- factor(group)
time <- factor(time)
id <- factor(id)
})

par(cex = .6)

with(demo3, interaction.plot(time, group, pulse,
ylim = c(10, 60), lty = c(1, 12), lwd = 3,
ylab = "mean of pulse", xlab = "time", trace.label = "group"))

demo3.aov <- aov(pulse ~ group * time + Error(id), data = demo3)
summary(demo3.aov)

Error: id
Df  Sum Sq Mean Sq F value    Pr(>F)
group      1 2035.04 2035.04  343.15 1.596e-06 ***
Residuals  6   35.58    5.93
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
Df  Sum Sq Mean Sq F value    Pr(>F)
time        2 2830.33 1415.17 553.761 1.517e-12 ***
group:time  2  200.33  100.17  39.196 5.474e-06 ***
Residuals  12   30.67    2.56
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#### Demo Analysis #4

The within subject test indicate that the interaction of time and group is significant. The main effect of time is not significant. However, the significant interaction indicates that the groups are changing over time and they are changing in different ways, in other words, in the graph the lines of the groups will not be parallel. The between groups test indicates that there the variable group is significant. In the graph for this particular case we see that one group is increasing in depression over time and the other group is decreasing in depression over time.


## Convert variables to factor
demo4 <- within(demo4, {
group <- factor(group)
time <- factor(time)
id <- factor(id)
})
par(cex = .6)

with(demo4, interaction.plot(time, group, pulse,
ylim = c(10, 60), lty = c(1, 12), lwd = 3,
ylab = "mean of pulse", xlab = "time", trace.label = "group"))

demo4.aov <- aov(pulse ~ group * time + Error(id), data = demo4)
summary(demo4.aov)

Error: id
Df  Sum Sq Mean Sq F value    Pr(>F)
group      1 2542.04 2542.04  628.96 2.646e-07 ***
Residuals  6   24.25    4.04
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
Df Sum Sq Mean Sq  F value    Pr(>F)
time        2    1.0    0.50   0.0789    0.9246
group:time  2 1736.3  868.17 137.0789 5.438e-09 ***
Residuals  12   76.0    6.33
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#### Exercise data

The data called exer, consists of people who were randomly assigned to two different diets: low-fat and not low-fat and three different types of exercise: at rest, walking leisurely and running. Their pulse rate was measured at three different time points during their assigned exercise: at 1 minute, 15 minutes and 30 minutes.


## Convert variables to factor
exer <- within(exer, {
diet <- factor(diet)
exertype <- factor(exertype)
time <- factor(time)
id <- factor(id)
})
print(exer)
id diet exertype pulse time
1   1    1        1    85    1
2   1    1        1    85    2
3   1    1        1    88    3
4   2    1        1    90    1
5   2    1        1    92    2
6   2    1        1    93    3
7   3    1        1    97    1
8   3    1        1    97    2
9   3    1        1    94    3
10  4    1        1    80    1
11  4    1        1    82    2
12  4    1        1    83    3
13  5    1        1    91    1
14  5    1        1    92    2
15  5    1        1    91    3
16  6    2        1    83    1
17  6    2        1    83    2
18  6    2        1    84    3
19  7    2        1    87    1
20  7    2        1    88    2
21  7    2        1    90    3
22  8    2        1    92    1
23  8    2        1    94    2
...


#### Exercise example, model 1 (time and diet)

Let us first consider the model including diet as the group variable. The graph would indicate that the pulse rate of both diet types increase over time but for the non-low fat group (diet=2) the pulse rate is increasing more over time than for the low fat group (diet=1)


par(cex=.6)

with(exer, interaction.plot(time, diet, pulse,
ylim = c(90, 110), lty = c(1, 12), lwd = 3,
ylab = "mean of pulse", xlab = "time", trace.label = "group"))



Looking at the results we conclude that the effect of time is significant but the interaction of time and diet is not significant. The between subject test of the effect of diet is also not significant. Consequently, in the graph we have lines that are not flat, in fact, they are actually increasing over time, which was expected since the effect of time was significant. Furthermore, the lines are approximately parallel which was anticipated since the interaction was not significant.


diet.aov <- aov(pulse ~ diet * time + Error(id),
data = exer)
summary(diet.aov)

Error: id
Df  Sum Sq Mean Sq F value  Pr(>F)
diet       1  1261.9 1261.88  3.1471 0.08694 .
Residuals 28 11227.0  400.97
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
Df Sum Sq Mean Sq F value    Pr(>F)
time       2 2066.6 1033.30 11.8078 5.264e-05 ***
diet:time  2  192.8   96.41  1.1017    0.3394
Residuals 56 4900.6   87.51
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#### Exercise example, model 2 (time and exercise type)

Next, let us consider the model including exertype as the group variable.


with(exer, interaction.plot(time, exertype, pulse,
ylim = c(80, 130), lty = c(1, 2, 4), lwd = 2,
ylab = "mean of pulse", xlab = "time"))



The interaction of time and exertype is significant as is the effect of time. The between subject test of the effect of exertype is also significant. Consequently, in the graph we have lines that are not parallel which we expected since the interaction was significant. Furthermore, we see that some of the lines that are rather far apart and at least one line is not horizontal which was anticipated since exertype and time were both significant.


exertype.aov <- aov(pulse ~ exertype * time + Error(id), data = exer)
summary(exertype.aov)

Error: id
Df Sum Sq Mean Sq F value   Pr(>F)
exertype   2 8326.1  4163.0  27.001 3.62e-07 ***
Residuals 27 4162.8   154.2
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
Df Sum Sq Mean Sq F value    Pr(>F)
time           2 2066.6 1033.30  23.543 4.446e-08 ***
exertype:time  4 2723.3  680.83  15.512 1.651e-08 ***
Residuals     54 2370.1   43.89
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#### Further Issues

Missing Data

• Compare aov and lme functions handling of  missing data (under construction).

#### Variance-Covariance Structures

Independence

As though analyzed using between subjects analysis.

s2
0 s2
0 0 s2

Compound Symmetry

Assumes that the variance-covariance structure has a single variance (represented by s2) for all 3 of the time points and a single covariance (represented by s1) for each of the pairs of trials.  This structure is illustrated by the half matrix below.

s2
s1 s2
s1 s1 s2

Unstructured

Assumes that each variance and covariance is unique.   Each trial has its own variance (e.g. s12 is the variance of trial 1) and each pair of trials has its own covariance (e.g. s21 is the covariance of trial 1 and trial2).  This structure is illustrated by the half matrix below.

s12
s21 s22
s31 s32 s32

Autoregressive

Another common covariance structure which is frequently observed in repeated measures data is an autoregressive structure, which recognizes that observations which are more proximate are more correlated than measures that are more distant.  This structure is illustrated by the half matrix below.

s2
sr s2
sr2 sr s2

Autoregressive Heterogeneous Variances

If the variances change over time, then the covariance would look like this.

s12
sr     s22
sr2   sr     s32

However, we cannot use this kind of covariance structure in a traditional repeated measures analysis (using the aov function), but we can use it in the gls function.

Let's look at the correlations, variances and covariances for the exercise data.


mat <- with(exer, matrix(c(pulse[time==1], pulse[time==2], pulse[time==3]), ncol = 3))
var(mat)

[,1]      [,2]      [,3]
[1,] 37.84368  48.78851  60.28506
[2,] 48.78851 212.11954 233.76092
[3,] 60.28506 233.76092 356.32299
cor(mat)
[,1]      [,2]      [,3]
[1,] 1.0000000 0.5445409 0.5191479
[2,] 0.5445409 1.0000000 0.8502755
[3,] 0.5191479 0.8502755 1.0000000


#### Exercise example, model 2 using the gls function

Even though we are very impressed with our results so far, we are not completely convinced that the variance-covariance structure really has compound symmetry. In order to compare models with different variance-covariance structures we have to use the gls function (gls = generalized least squares) and try the different structures that we think our data might have.

Compound Symmetry

The first model we will look at is one using compound symmetry for the variance-covariance structure. This model should confirm the results of the results of the tests that we obtained through the aov function and we will be able to obtain fit statistics which we will use for comparisons with our models that assume other variance-covariance structures.

In order to use the gls function we need to include the repeated structure in our data set object. We do this by using the groupedData function and the id variable following the bar notation indicates that observations are repeated within id. We then fit the model using the gls function and we use the corCompSymm function in the corr argument because we want to use compound symmetry.  We obtain the 95% confidence intervals for the parameter estimates, the estimate of rho and the estimated of the standard error of the residuals by using the intervals function.


library(nlme)
longg <- groupedData(pulse ~ as.numeric(exertype) * as.numeric(time) | id, data = exer)
fit.cs <- gls(pulse ~ exertype * time, data = longg,
corr = corCompSymm(, form= ~ 1 | id) )
summary(fit.cs)

Generalized least squares fit by REML
Model: pulse ~ exertype * time
Data: longg
AIC      BIC    logLik
612.8316 639.1706 -295.4158

Correlation Structure: Compound symmetry
Formula: ~1 | id
Parameter estimate(s):
Rho
0.455816


Unstructured

We now try an unstructured covariance matrix. Option "corr = corSymm" specifies that the correlation structure is unstructured. Option "weights = varident(form = ~ 1 | time)" specifies that the variance at each time point can be different.


fit.un <- gls(pulse ~ exertype * time, data = longg,
corr=corSymm(form = ~ 1 | id),
weights = varIdent(form = ~ 1 | time))
summary(fit.un)

Generalized least squares fit by REML
Model: pulse ~ exertype * time
Data: longg
AIC      BIC    logLik
607.7365 643.6532 -288.8682

Correlation Structure: General
Formula: ~1 | id
Parameter estimate(s):
Correlation:
1     2
2 0.434
3 0.417 0.583
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | time
Parameter estimates:
1        2        3
1.000000 1.596720 1.877599

anova(fit.un)
Denom. DF: 81
numDF  F-value p-value
(Intercept)       1 8184.123  <.0001
exertype          2    6.426  0.0026
time              2   22.324  <.0001
exertype:time     4   14.387  <.0001


Autoregressive

From previous studies we suspect that our data might actually have an auto-regressive variance-covariance structure so this is the model we will look at next. However, for our data the auto-regressive variance-covariance structure does not fit our data much better than the compound symmetry does.


fit.ar1 <- gls(pulse ~ exertype * time, data = longg,
corr = corAR1(, form= ~ 1 | id))
summary(fit.ar1)

Generalized least squares fit by REML
Model: pulse ~ exertype * time
Data: longg
AIC      BIC    logLik
612.1163 638.4553 -295.0582

Correlation Structure: AR(1)
Formula: ~1 | id
Parameter estimate(s):
Phi
0.4992423

anova(fit.ar1)
Denom. DF: 81
numDF  F-value p-value
(Intercept)       1 6167.352  <.0001
exertype          2   26.990  <.0001
time              2   18.196  <.0001
exertype:time     4   11.733  <.0001


Autoregressive with heterogeneous variances

Now we suspect that what is actually going on is that the we have auto-regressive covariances and heterogeneous variances.


fit.arh1 <- gls(pulse ~ exertype * time, data = longg,
corr = corAR1(, form = ~ 1 | id), weight = varIdent(form = ~ 1 | time))
summary(fit.arh1)

Generalized least squares fit by REML
Model: pulse ~ exertype * time
Data: longg
AIC      BIC    logLik
605.7693 636.8971 -289.8846

Correlation Structure: AR(1)
Formula: ~1 | id
Parameter estimate(s):
Phi
0.5100781
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | time
Parameter estimates:
1        2        3
1.000000 1.561315 1.796993

anova(fit.arh1)
Denom. DF: 81
numDF  F-value p-value
(Intercept)       1 8284.813  <.0001
exertype          2    9.134   3e-04
time              2   21.918  <.0001
exertype:time     4   13.805  <.0001


Model comparison (using the anova function)

We can use the anova function to compare competing models to see which model fits the data best.

anova(fit.cs, fit.un)
Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit.cs     1 11 612.8316 639.1706 -295.4158
fit.un     2 15 607.7365 643.6532 -288.8682 1 vs 2 13.09512  0.0108

anova(fit.cs, fit.ar1)
Model df      AIC      BIC    logLik
fit.cs      1 11 612.8316 639.1706 -295.4158
fit.ar1     2 11 612.1163 638.4553 -295.0582

anova(fit.cs, fit.arh1)
Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit.cs       1 11 612.8316 639.1706 -295.4158
fit.arh1     2 13 605.7693 636.8971 -289.8846 1 vs 2 11.06236   0.004


The two most promising structures are Autoregressive Heterogeneous Variances and Unstructured since these two models have the smallest AIC values and the -2 Log Likelihood scores are significantly smaller than the -2 Log Likelihood scores of other models.

#### Exercise example, model 3 (time, diet and exertype)---using the aov function

Looking at models including only the main effects of diet or exertype separately does not answer all our questions. We would also like to know if the people on the low-fat diet who engage in running have lower pulse rates than the people participating in the not low-fat diet who are not running. In order to address these types of questions we need to look at a model that includes the interaction of diet and exertype. After all the analysis involving the variance-covariance structures we will look at this model using both functions aov and gls.

In the graph of exertype by diet we see that for the low-fat diet (diet=1) group the pulse rate for the two exercise types: at rest and walking, are very close together, indeed they are almost flat, whereas the running group has a higher pulse rate that increases over time. For the not low-fat diet (diet=2) group the same two exercise types: at rest and walking, are also very close together and almost flat. For this group, however, the pulse rate for the running group increases greatly over time and the rate of increase is much steeper than the increase of the running group in the low-fat diet group.
The within subject tests indicate that there is a three-way interaction between diet, exertype and time. In other words, the pulse rate will depend on which diet you follow, the exercise type you engage in and at what time during the the exercise that you measure the pulse. The interactions of time and exertype and diet and exertype are also significant as are the main effects of diet and exertype.


par(cex = .6)
with(exer, interaction.plot(time[diet==1], exertype[diet==1], pulse[diet==1],
ylim = c(80, 150), lty = c(1, 12, 8),
trace.label = "exertype", ylab = "mean of pulse", xlab = "time"))
title("Diet = 1")

with(exer, interaction.plot(time[diet==2], exertype[diet==2], pulse[diet==2],
ylim = c(80, 150), lty = c(1, 12, 8),
trace.label = "exertype", ylab = "mean of pulse", xlab = "time"))
title("Diet = 2")



Looking at the graphs of exertype by diet.


both.aov <- aov(pulse ~ exertype * diet * time + Error(id), data = exer)
summary(both.aov)

Error: id
Df Sum Sq Mean Sq F value    Pr(>F)
exertype       2 8326.1  4163.0 47.9152 4.166e-09 ***
diet           1 1261.9  1261.9 14.5238 0.0008483 ***
exertype:diet  2  815.8   407.9  4.6945 0.0190230 *
Residuals     24 2085.2    86.9
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
Df  Sum Sq Mean Sq F value    Pr(>F)
time                2 2066.60 1033.30 31.7206 1.662e-09 ***
exertype:time       4 2723.33  680.83 20.9005 4.992e-10 ***
diet:time           2  192.82   96.41  2.9597   0.06137 .
exertype:diet:time  4  613.64  153.41  4.7095   0.00275 **
Residuals          48 1563.60   32.57
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#### Exercise example, model 3 (time, diet and exertype)--using the gls fuction

For the gls model we will use the autoregressive heterogeneous variance-covariance structure since we previously observed that this is the structure that appears to fit the data the best (see discussion of variance-covariance structures). We do not expect to find a great change in which factors will be significant but we do expect to have a model that has a better fit than the anova model.
The graphs are exactly the same as the anova model and we find that the same factors are significant. However, since the model has a better fit we can be more confident in the estimate of the standard errors and therefore we can be more confident in the tests and in the findings of significant factors. The model has a better fit than the model only including exertype and time because both the -2Log Likelihood and the AIC has decrease dramatically. The -2 Log Likelihood decreased from 579.8 for the model including only exertype and time to 505.3 for the current model.


longa <- groupedData(pulse ~ as.numeric(exertype) * as.numeric(diet) * as.numeric(time) | id,
data = exer)
both.arh1 <- gls(pulse ~ exertype * diet * time, data = longa,
corr = corAR1(, form = ~ 1 | id), weight = varIdent(form = ~ 1 | time))
summary(both.arh1)

Generalized least squares fit by REML
Model: pulse ~ exertype * diet * time
Data: longa
AIC      BIC    logLik
549.2788 599.3654 -252.6394

Correlation Structure: AR(1)
Formula: ~1 | id
Parameter estimate(s):
Phi
0.360999
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | time
Parameter estimates:
1        2        3
1.000000 1.490617 1.171196
anova(both.arh1)
Denom. DF: 72
numDF   F-value p-value
(Intercept)            1 13391.413  <.0001
exertype               2    42.121  <.0001
diet                   1    17.330  0.0001
time                   2    30.822  <.0001
exertype:diet          2     4.738  0.0117
exertype:time          4    20.248  <.0001
diet:time              2     2.797  0.0676
exertype:diet:time     4     4.452  0.0029


#### Contrasts and interaction contrasts for model 3

From the graphs in the above analysis we see that the runners (exertype level 3) have a pulse rate that is increases much quicker than the pulse rates of the two other groups. We would like to know if there is a statistically significant difference between the changes over time in the pulse rate of the runners versus the change over time in the pulse rate of the walkers and the people at rest across diet groups and across time. Furthermore, we suspect that there might be a difference in pulse rate over time and across exercise type between the two diet groups. But to make matters even more complicated we would like to test if the runners in the low fat diet group are statistically significantly different from all the other groups (i.e. the runners in the non-low fat diet, the walkers and the people at rest in both diet groups). Since we are being ambitious we also want to test if the runners in the low fat diet group (diet=1) are different from the runners in the non-low fat diet group (diet=2).

In order to implement contrasts coding for diet and exertype we will make copies of the variables. If they were not already factors, we would need to convert them to factors first. Note that we are still using the data frame longa which has the hierarchy characteristic that we need for the gls function.


longa[, c("ef", "df", "tf")] <- longa[, c("exertype", "diet", "time")]


Now we can attach the contrasts to the factor variables using the contrasts function. We need to use the contrast coding for regression which is discussed in the chapter 6 in our regression web book (note that the coding system is not package specific so we arbitrarily choose to link to the SAS web book.) For the contrast coding of ef and tf we first create the matrix containing the contrasts and then we assign the contrasts to them. The contrasts coding for df is simpler since there are just two levels and we can therefore assign the contrasts directly without having to create a matrix of contrasts.


m <- matrix( c( c(-1/2, 1/2, 0), c(-1/3, -1/3, 2/3) ), ncol=2)
contrasts(longa$ef) <- m (contrasts(longa$tf) <- m)

[,1]       [,2]
[1,] -0.5 -0.3333333
[2,]  0.5 -0.3333333
[3,]  0.0  0.6666667
(contrasts(longa$df) <- c(-1/2, 1/2)) [1] -0.5 0.5  Now that we have all the contrast coding we can finally run the model. Looking at the results the variable ef1 corresponds to the contrast of exertype=1 versus exertype=2 and it is not significant indicating that there is no difference between the pulse rate of the people at rest and the people who walk leisurely. The variable ef2 corresponds to the contrast of exertype=3 versus the average of exertype=1 and exertype=2. This contrast is significant indicating that there is a difference between the mean pulse rate of the runners compared to the walkers and the people at rest. The variable df1 corresponds to the contrast of the two diets and it is significant indicating that the mean pulse rate of the people on the low-fat diet is different from that of the people on a non-low fat diet. The interaction ef2:df1 corresponds to the contrast of the runners on a low fat diet (people who are in the group exertype=3 and diet=1) versus everyone else. This contrast is significant indicating the the mean pulse rate of the runners on a low fat diet is different from everyone else's mean pulse rate.  model.cs <- gls(pulse ~ ef * df * tf, data = longa, corr = corCompSymm(, form = ~ 1 | id) ) summary(model.cs) Generalized least squares fit by REML Model: pulse ~ ef * df * tf Data: longa AIC BIC logLik 547.6568 593.1901 -253.8284 Correlation Structure: Compound symmetry Formula: ~1 | id Parameter estimate(s): Rho 0.3572133 Coefficients: Value Std.Error t-value p-value (Intercept) 99.70000 0.982533 101.47246 0.0000 ef1 4.36667 2.406704 1.81438 0.0738 ef2 20.05000 2.084266 9.61969 0.0000 df1 7.48889 1.965065 3.81101 0.0003 tf1 8.40000 1.473658 5.70010 0.0000 tf2 7.10000 1.276225 5.56328 0.0000 ef1:df1 0.46667 4.813407 0.09695 0.9230 ef2:df1 12.76667 4.168533 3.06263 0.0031 ef1:tf1 2.80000 3.609709 0.77569 0.4405 ef2:tf1 18.90000 3.126100 6.04587 0.0000 ef1:tf2 0.20000 3.126100 0.06398 0.9492 ef2:tf2 18.45000 2.707282 6.81495 0.0000 df1:tf1 2.93333 2.947315 0.99526 0.3229 df1:tf2 5.66667 2.552450 2.22009 0.0296 ef1:df1:tf1 -0.40000 7.219418 -0.05541 0.9560 ef2:df1:tf1 11.20000 6.252200 1.79137 0.0774 ef1:df1:tf2 -3.40000 6.252200 -0.54381 0.5883 ef2:df1:tf2 21.20000 5.414564 3.91537 0.0002  The contrasts that we were not able to obtain in the previous code were the tests of the simple effects, i.e. testing for difference between the two diets at exertype=3. We would like to test the difference in mean pulse rate of the people following the two diets at a specific level of exertype. We would like to know if there is a difference in the mean pulse rate for runners (exertype=3) in the lowfat diet (diet=1) versus the runners in the non-low fat diet (diet=2). In order to obtain this specific contrasts we need to code the contrasts for diet at each level of exertype and include these in the model. For more explanation of why this is the case we strongly urge you to read chapter 5 in our web book that we mentioned before. Looking at the results the variable e3d12 corresponds to the contrasts of the runners on the low fat diet versus the runners on the non-low fat diet. This contrast is significant indicating that the mean pulse rate of runners on the low fat diet is different from that of the runners on a non-low fat diet.  longa$e1d12 <- (-1/2*(longa$exertype==1 & longa$diet==1))
longa$e1d12[longa$exertype==1 & longa$diet==2] <- 1/2 longa$e2d12 <- (-1/2*(longa$exertype==1)) longa$e2d12[longa$exertype==2 & longa$diet==2] <- 1/2

longa$e3d12 <- (-1/2*(longa$exertype==3 & longa$diet==1)) longa$e3d12[longa$exertype==3 & longa$diet==2] <- 1/2

modela.cs <- gls(pulse ~ ef + e1d12 + e2d12 + e3d12 , data = longa,
corr = corCompSymm(, form = ~ 1 | id) )
summary(modela.cs)

Coefficients:
Value Std.Error  t-value p-value
(Intercept) 100.27778  1.134531 88.38699  0.0000
ef1          -0.83333  5.644220 -0.14764  0.8830
ef2          19.18333  2.251265  8.52113  0.0000
e1d12         3.00000  3.403593  0.88142  0.3806
e2d12         6.93333  6.807186  1.01853  0.3114
e3d12        16.00000  3.403593  4.70091  0.0000


### Unequally Spaced Time Points

#### Modeling Time as a Linear Predictor of Pulse

We have another study which is very similar to the one previously discussed except that in this new study the pulse measurements were not taken at regular time points.  In this study a baseline pulse measurement was obtained at time = 0 for every individual in the study. However, subsequent pulse measurements were taken at less regular time intervals.  The second pulse measurements were taken at approximately 2 minutes (time = 120 seconds); the pulse measurement was obtained at approximately 5 minutes (time = 300 seconds); and the fourth and final pulse measurement was obtained at approximately 10 minutes (time = 600 seconds). The data for this study is displayed below.


study2 <- within(study2, {
id <- factor(id)
exertype <- factor(exertype)
diet <- factor(diet)
})
study2[1:20, ]

id exertype diet pulse time
1   1        1    1    90    0
2   1        1    1    92  228
3   1        1    1    93  296
4   1        1    1    93  639
5   2        1    1    90    0
6   2        1    1    92   56
7   2        1    1    93  434
8   2        1    1    93  538
9   3        1    1    97    0
10  3        1    1    97  150
11  3        1    1    94  295
12  3        1    1    94  541
13  4        1    1    80    0
14  4        1    1    82  121
15  4        1    1    83  256
16  4        1    1    83  575
17  5        1    1    91    0
18  5        1    1    92  161
19  5        1    1    91  252
20  5        1    1    91  526


In order to get a better understanding of the data we will look at a scatter plot of the data with lines connecting the points for each individual.


library(lattice)

##
par(cex = .6)
xyplot(pulse ~ time, data = study2, groups = id,
type = "o", panel = panel.superpose)

xyplot(pulse ~ time | exertype, data = study2, groups = id,
type = "o", panel = panel.superpose)

xyplot(pulse ~ time | diet, data = study2, groups = id,
type = "o", panel = panel.superpose)



This is a situation where multilevel modeling excels for the analysis of data with irregularly spaced time points.  The multilevel model with time as a linear effect is illustrated in the following equations.

Level 1 (time): Pulse = β0j + β1j (Time) + rij
Level 2 (person): β0j =  γ00  + γ01(Exertype) + u0j
Level 2 (person): β1j =  γ10  + γ11(Exertype) + u1j

Substituting the level 2 model into the level 1 model we get the following single equations. Note: The random components have been placed in square brackets.

Pulse = γ00 + γ01(Exertype) + γ10(Time) + γ11(Exertype*time) + [ u0j + u1j(Time) + rij ]

Since this model contains both fixed and random components, it can be analyzed using the lme function as shown below.


time.linear <- lme(pulse ~ exertype * time,
random = list(id = pdDiag(~ time)), data = study2)
summary(time.linear)

Linear mixed-effects model fit by REML
Data: study2
AIC      BIC    logLik
856.8227 881.4485 -419.4113

Random effects:
Formula: ~time | id
Structure: Diagonal
(Intercept)       time Residual
StdDev:    5.821602 0.01151569 5.692529

Fixed effects: pulse ~ exertype * time
Value Std.Error DF  t-value p-value
(Intercept)    91.07179  2.274185 87 40.04591  0.0000
exertype2       3.51075  3.220818 27  1.09002  0.2853
exertype3      12.62517  3.226279 27  3.91323  0.0006
time            0.00158  0.005244 87  0.30185  0.7635
exertype2:time  0.00716  0.007599 87  0.94272  0.3484
exertype3:time  0.05477  0.007531 87  7.27233  0.0000
Correlation:
(Intr) exrty2 exrty3 time   exrt2:
exertype2      -0.706
exertype3      -0.705  0.498
time           -0.311  0.220  0.220
exertype2:time  0.215 -0.320 -0.152 -0.690
exertype3:time  0.217 -0.153 -0.320 -0.696  0.481

Standardized Within-Group Residuals:
Min          Q1         Med          Q3         Max
-2.51744499 -0.31571963 -0.05192003  0.26453221  3.19537209

Number of Observations: 120
Number of Groups: 30

anova(time.linear)
numDF denDF  F-value p-value
(Intercept)       1    87 6373.823  <.0001
exertype          2    27   24.229  <.0001
time              1    87   49.954  <.0001
exertype:time     2    87   30.677  <.0001


Graphs of predicted values. The first graph shows just the lines for the predicted values one for each level of exertype. It is obvious that the straight lines do not approximate the data very well, especially for exertype group 3. The rest of the graphs show the predicted values as well as the observed values. The predicted values are the darker straight lines; the line for exertype group 1 is blue, for exertype group 2 it is red and for exertype group 3 the line is green. In this graph it becomes even more obvious that the model does not fit the data very well.


fitted <- fitted(time.linear, level=0)

with(study2, plot(time[exertype==3], fitted[exertype==3], ylim = c(50, 150),
xlab = "time", ylab = "predicted", type = "b", col = "green"))
with(study2, points(time[exertype==2], fitted[exertype==2],
pch = 4, type = "b", col = "red"))
with(study2, points(time[exertype==1], fitted[exertype==1],
pch = 16, type = "b", col = "blue"))

xyplot(pulse[exertype==1] ~ time[exertype==1], data = study2, groups = id,
type = "o", ylim = c(50, 150), xlim = c(0, 800),
panel = panel.superpose, col = "blue")
with(study2, lines(time[exertype==1], fitted[exertype==1],
ylim = c(50, 150),  xlim = c(0, 800),
type = "b", col = "dark blue", lwd = 4))

xyplot(pulse[exertype==2] ~ time[exertype==2], data = study2, groups=id,
type = "o", ylim = c(50, 150), xlim = c(0, 800),
panel = panel.superpose, col = "red")
with(study2, lines(time[exertype==2], fitted[exertype==2],
ylim = c(50, 150),  xlim = c(0, 800),
type = "b", col = "dark red", lwd = 4))

xyplot(pulse[exertype==3] ~ time[exertype==3], data = study2, groups = id,
type = "o", ylim = c(50, 150), xlim = c(0, 800),
panel = panel.superpose, col = "green")
with(study2, lines(time[exertype==3], fitted[exertype==3],
ylim = c(50, 150), xlim = c(0, 800),
type = "b", col = "dark green", lwd = 4))



#### Modeling Time as a Quadratic Predictor of Pulse

To model the quadratic effect of time, we add time*time to the model. We see that term is significant.


study2$time2 <- study2$time^2
time.quad <- lme(pulse ~ exertype * time + time2,
random = list(id = pdDiag(~ time)), study2)

Linear mixed-effects model fit by REML
Data: study2
AIC      BIC    logLik
859.3578 886.6317 -419.6789

Random effects:
Formula: ~time | id
Structure: Diagonal
(Intercept)      time Residual
StdDev:    5.763908 0.0122891 4.981454

Fixed effects: pulse ~ exertype * time + time2
Value Std.Error DF  t-value p-value
(Intercept)    88.75438 2.2189692 86 39.99802  0.0000
exertype2       3.56744 3.0669444 27  1.16319  0.2549
exertype3      12.92326 3.0722775 27  4.20641  0.0003
time            0.03524 0.0086532 86  4.07252  0.0001
time2          -0.00005 0.0000106 86 -4.82532  0.0000
exertype2:time  0.00563 0.0073822 86  0.76218  0.4480
exertype3:time  0.05253 0.0073320 86  7.16392  0.0000

numDF denDF  F-value p-value
(Intercept)       1    86 6717.113  <.0001
exertype          2    27   22.247  <.0001
time              1    86   53.952  <.0001
time2             1    86   28.249  <.0001
exertype:time     2    86   30.481  <.0001


Graphs of predicted values. The first graph shows just the lines for the predicted values one for each level of exertype. The curved lines approximate the data better than the straight lines of the model with time as a linear predictor. The rest of the graphs show the predicted values as well as the observed values. The predicted values are the very curved darker lines; the line for exertype group 1 is blue, for exertype group 2 it is orange and for exertype group 3 the line is green. This model fits the data better, but it appears that the predicted values for the exertype group 3 have too little curvature and the predicted values for exertype groups 1 and 2 have too much curvature.


fitted2 <- fitted(time.quad, level = 0)
a <- with(study2,
data.frame(time, fitted2, exertype)[order(exertype, time), ])

with(a, {
plot(time[exertype==3], fitted2[exertype==3], ylim = c(50, 150),
xlab = "time", ylab = "predicted", col = "green", type = "b")
points(time[exertype==2], fitted2[exertype==2],
pch = 4, col = "red", type = "b")
points(time[exertype==1], fitted2[exertype==1],
pch = 16, col = "blue", type = "b")

xyplot(pulse[exertype==1] ~ time[exertype==1], groups = id, data = study2,
ylim = c(50, 150), xlim = c(0, 800), type = "o",
panel=panel.superpose, col="blue")
with(a, lines(time[exertype==1], fitted2[exertype==1],
ylim = c(50, 150), xlim = c(0, 800),
type = "b", col = "dark blue", lwd = 4))

xyplot(pulse[exertype==2] ~ time[exertype==2], groups = id, data = study2,
ylim=c(50, 150), xlim=c(0, 800), type="o",
panel=panel.superpose, col="red")
with(a, lines(time[exertype==2], fitted2[exertype==2],
ylim = c(50, 150), xlim = c(0, 800),
type = "b", col = "dark red", lwd = 4))

xyplot(pulse[exertype==3] ~time[exertype==3], groups = id, data = study2,
ylim = c(50, 150), xlim = c(0, 800), type = "o",
panel = panel.superpose, col = "green")
with(a, lines(time[exertype==3], fitted2[exertype==3],
ylim = c(50, 150), xlim = c(0, 800),
type = "b", col = "dark green", lwd = 4))



#### Modeling Time as a Quadratic Predictor of Pulse, Interacting by Exertype

We can include an interaction of time*time*exertype to indicate that the different exercises not only show different linear trends over time, but that they also show different quadratic trends over time, as shown below.  The time*time*exertype term is significant.


time.quad2 <- lme(pulse ~ exertype * time + exertype * time2,
random = list(id = pdDiag(~ time)), data = study2)

Linear mixed-effects model fit by REML
Data: study2
AIC      BIC    logLik
864.3194 896.8337 -420.1597

Random effects:
Formula: ~time | id
Structure: Diagonal
(Intercept)       time Residual
StdDev:    5.996351 0.01192684 3.869181

Fixed effects: pulse ~ exertype * time + exertype * time2
Value Std.Error DF  t-value p-value
(Intercept)     90.81505 2.1902570 84 41.46319  0.0000
exertype2        2.66056 3.1027921 27  0.85747  0.3987
exertype3        7.28070 3.0989266 27  2.34943  0.0264
time             0.00554 0.0099971 84  0.55411  0.5810
time2           -0.00001 0.0000136 84 -0.45396  0.6510
exertype2:time   0.01889 0.0142534 84  1.32521  0.1887
exertype3:time   0.13928 0.0146086 84  9.53418  0.0000
exertype2:time2 -0.00002 0.0000198 84 -0.94590  0.3469
exertype3:time2 -0.00014 0.0000208 84 -6.66757  0.0000

numDF denDF  F-value p-value
(Intercept)        1    84 6779.738  <.0001
exertype           2    27   19.074  <.0001
time               1    84   69.223  <.0001
time2              1    84   43.847  <.0001
exertype:time      2    84   38.968  <.0001
exertype:time2     2    84   24.767  <.0001


Graphs of predicted values. The first graph shows just the lines for the predicted values one for each level of exertype. The lines now have different degrees of curvature which approximates the data much better than the other two models. The rest of graphs show the predicted values as well as the observed values. The line for exertype group 1 is blue, for exertype group 2 it is orange and for exertype group 3 the line is green. This model fits the data the best with more curvature for exertype group 3 and less curvature for exertype groups 1 and 2.


fitted3 <- fitted(time.quad2, level = 0)
a <- with(study2,
data.frame(time, fitted3, exertype)[order(exertype, time), ])

with(a, {
plot(time[exertype==3], fitted3[exertype==3], ylim = c(50, 150),
xlab = "time", ylab = "predicted", col = "green", type = "b")
points(time[exertype==2], fitted3[exertype==2],
pch = 4, col = "red", type = "b")
points(time[exertype==1], fitted3[exertype==1],
pch = 16,  col = "blue", type = "b")

xyplot(pulse[exertype==1] ~ time[exertype==1], groups = id, data = study2,
ylim = c(50, 150), xlim = c(0, 800), type = "o",
panel = panel.superpose, col = "blue")
with(a, lines(time[exertype==1], fitted3[exertype==1],
ylim = c(50, 150), xlim = c(0, 800),
type = "b", col = "dark blue", lwd = 4))

xyplot(pulse[exertype==2] ~ time[exertype==2], groups = id, data = study2,
ylim = c(50, 150), xlim = c(0, 800), type = "o",
panel = panel.superpose, col = "red")
with(a, lines(time[exertype==2], fitted3[exertype==2],
ylim = c(50, 150), xlim = c(0, 800),
type = "b", col = "dark red", lwd = 4))

xyplot(pulse[exertype==3] ~ time[exertype==3], groups = id,  data = study2,
ylim = c(50, 150), xlim = c(0, 800), type = "o",
panel = panel.superpose, col = "green")
with(a, lines(time[exertype==3], fitted3[exertype==3],
ylim = c(50, 150), xlim = c(0, 800),
type = "b", col = "dark green", lwd = 4))



The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.