Chapter 24 R Code for the Repeated Measures
#Read the data
d=read.delim("http://dnett.github.io/S510/RepeatedMeasures.txt")
#Create Factors
d$Program=as.factor(d$Program)
d$Subj=as.factor(d$Subj)
d$Timef=as.factor(d$Time)
#Load the nlme package
library(nlme)
# compound symmetry
lme(Strength ~ Program * Timef, data = d,random = ~1|Subj)## Linear mixed-effects model fit by REML
## Data: d
## Log-restricted-likelihood: -710.4101
## Fixed: Strength ~ Program * Timef
## (Intercept) Program2 Program3 Timef4 Timef6 Timef8
## 79.687500000 1.360119048 0.062500000 0.875000000 1.125000000 1.312500000
## Timef10 Timef12 Timef14 Program2:Timef4 Program3:Timef4 Program2:Timef6
## 1.562500000 1.437500000 1.437500000 -0.255952381 -0.675000000 -0.267857143
## Program3:Timef6 Program2:Timef8 Program3:Timef8 Program2:Timef10 Program3:Timef10 Program2:Timef12
## -0.875000000 0.163690476 -1.012500000 0.008928571 -1.512500000 0.229166667
## Program3:Timef12 Program2:Timef14 Program3:Timef14
## -1.587500000 0.610119048 -1.587500000
##
## Random effects:
## Formula: ~1 | Subj
## (Intercept) Residual
## StdDev: 3.098924 1.094017
##
## Number of Observations: 399
## Number of Groups: 57gls(Strength ~ Program * Timef, data = d,
correlation = corCompSymm(form = ~1|Subj))## Generalized least squares fit by REML
## Model: Strength ~ Program * Timef
## Data: d
## Log-restricted-likelihood: -710.4101
##
## Coefficients:
## (Intercept) Program2 Program3 Timef4 Timef6 Timef8
## 79.687500000 1.360119048 0.062500000 0.875000000 1.125000000 1.312500000
## Timef10 Timef12 Timef14 Program2:Timef4 Program3:Timef4 Program2:Timef6
## 1.562500000 1.437500000 1.437500000 -0.255952381 -0.675000000 -0.267857143
## Program3:Timef6 Program2:Timef8 Program3:Timef8 Program2:Timef10 Program3:Timef10 Program2:Timef12
## -0.875000000 0.163690476 -1.012500000 0.008928571 -1.512500000 0.229166667
## Program3:Timef12 Program2:Timef14 Program3:Timef14
## -1.587500000 0.610119048 -1.587500000
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | Subj
## Parameter estimate(s):
## Rho
## 0.8891805
## Degrees of freedom: 399 total; 378 residual
## Residual standard error: 3.286366# AR(1)
gls(Strength ~ Program * Timef, data = d,
correlation = corAR1(form = ~1|Subj))## Generalized least squares fit by REML
## Model: Strength ~ Program * Timef
## Data: d
## Log-restricted-likelihood: -633.4018
##
## Coefficients:
## (Intercept) Program2 Program3 Timef4 Timef6 Timef8
## 79.687500000 1.360119048 0.062500000 0.875000000 1.125000000 1.312500000
## Timef10 Timef12 Timef14 Program2:Timef4 Program3:Timef4 Program2:Timef6
## 1.562500000 1.437500000 1.437500000 -0.255952381 -0.675000000 -0.267857143
## Program3:Timef6 Program2:Timef8 Program3:Timef8 Program2:Timef10 Program3:Timef10 Program2:Timef12
## -0.875000000 0.163690476 -1.012500000 0.008928571 -1.512500000 0.229166667
## Program3:Timef12 Program2:Timef14 Program3:Timef14
## -1.587500000 0.610119048 -1.587500000
##
## Correlation Structure: AR(1)
## Formula: ~1 | Subj
## Parameter estimate(s):
## Phi
## 0.9517769
## Degrees of freedom: 399 total; 378 residual
## Residual standard error: 3.280242# general structure
gls(Strength ~ Program * Timef, data = d,
correlation = corSymm(form = ~1|Subj),
weight = varIdent(form = ~1|Timef))## Generalized least squares fit by REML
## Model: Strength ~ Program * Timef
## Data: d
## Log-restricted-likelihood: -617.4479
##
## Coefficients:
## (Intercept) Program2 Program3 Timef4 Timef6 Timef8
## 79.687500000 1.360119048 0.062500000 0.875000000 1.125000000 1.312500000
## Timef10 Timef12 Timef14 Program2:Timef4 Program3:Timef4 Program2:Timef6
## 1.562500000 1.437500000 1.437500000 -0.255952381 -0.675000000 -0.267857143
## Program3:Timef6 Program2:Timef8 Program3:Timef8 Program2:Timef10 Program3:Timef10 Program2:Timef12
## -0.875000000 0.163690476 -1.012500000 0.008928571 -1.512500000 0.229166667
## Program3:Timef12 Program2:Timef14 Program3:Timef14
## -1.587500000 0.610119048 -1.587500000
##
## Correlation Structure: General
## Formula: ~1 | Subj
## Parameter estimate(s):
## Correlation:
## 1 2 3 4 5 6
## 2 0.960
## 3 0.925 0.940
## 4 0.872 0.877 0.956
## 5 0.842 0.860 0.937 0.960
## 6 0.809 0.827 0.898 0.909 0.951
## 7 0.797 0.792 0.876 0.887 0.917 0.953
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Timef
## Parameter estimates:
## 2 4 6 8 10 12 14
## 1.000000 1.038707 1.104340 1.071323 1.173689 1.157118 1.203159
## Degrees of freedom: 399 total; 378 residual
## Residual standard error: 2.963148AIC and BIC for repeated measures in R:
- \(AIC = -2\ell(\hat\theta) + 2k\)
- \(BIC = -2\ell(\hat\theta) + k\ln(n)\)
- k = number of mean parameters (rank of X) + number of variance parameters
- n = total number of observations - rank of X