> #Repeated measures Example 
> 
> library(nlme)

Attaching package: 'nlme'


	The following object(s) are masked _by_ .GlobalEnv :

	 RatPupWeight 

> 
> #help(ergoStool)
> #Devore (2000) cites data from an article in _Ergometrics_ (1993,
> #pp. 519-535) on ``The Effects of a Pneumatic Stool and a
> #One-Legged Stool on Lower Limb Joint Load and Muscular Activity.''
> #This data frame contains the following columns:
> #
> #     effort a numeric vector giving the effort (Borg scale) required to
> #          arise from a stool.
> #
> #     Type a factor with levels 'T1',  'T2',  'T3', and  'T4' giving the
> #          stool type.
> #
> #     Subject an ordered factor giving a unique identifier for the
> #          subject in the experiment.
> 
> ergoStool
Grouped Data: effort ~ Type | Subject
   effort Type Subject
1      12   T1       1
2      15   T2       1
3      12   T3       1
4      10   T4       1
5      10   T1       2
6      14   T2       2
7      13   T3       2
8      12   T4       2
9       7   T1       3
10     14   T2       3
11     13   T3       3
12      9   T4       3
13      7   T1       4
14     11   T2       4
15     10   T3       4
16      9   T4       4
17      8   T1       5
18     11   T2       5
19      8   T3       5
20      7   T4       5
21      9   T1       6
22     11   T2       6
23     11   T3       6
24     10   T4       6
25      8   T1       7
26     12   T2       7
27     12   T3       7
28     11   T4       7
29      7   T1       8
30     11   T2       8
31      8   T3       8
32      7   T4       8
33      9   T1       9
34     13   T2       9
35     10   T3       9
36      8   T4       9
> 
> pdf("ergoplots.pdf")
> par(pty="s")
> plot(ergoStool)
> 
> #visually compare magnitude of the effects of Type and Subject factors
> plot.design(ergoStool)
> 
> fm1 <- lme(effort~Type, data=ergoStool, random=~1 | Subject)
> anova(fm1)
            numDF denDF  F-value p-value
(Intercept)     1    24 455.0075  <.0001
Type            3    24  22.3556  <.0001
> summary(fm1)
Linear mixed-effects model fit by REML
 Data: ergoStool 
       AIC      BIC   logLik
  133.1308 141.9252 -60.5654

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    1.332465 1.100295

Fixed effects: effort ~ Type 
               Value Std.Error DF   t-value p-value
(Intercept) 8.555556 0.5760123 24 14.853079  0.0000
TypeT2      3.888889 0.5186838 24  7.497610  0.0000
TypeT3      2.222222 0.5186838 24  4.284348  0.0003
TypeT4      0.666667 0.5186838 24  1.285304  0.2110
 Correlation: 
       (Intr) TypeT2 TypeT3
TypeT2 -0.45               
TypeT3 -0.45   0.50        
TypeT4 -0.45   0.50   0.50 

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.80200345 -0.64316591  0.05783115  0.70099706  1.63142054 

Number of Observations: 36
Number of Groups: 9 
> 
> #examine the fit and residuals
> intervals(fm1)
Approximate 95% confidence intervals

 Fixed effects:
                 lower      est.    upper
(Intercept)  7.3667247 8.5555556 9.744386
TypeT2       2.8183781 3.8888889 4.959400
TypeT3       1.1517114 2.2222222 3.292733
TypeT4      -0.4038442 0.6666667 1.737177
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Subject 
                   lower     est.    upper
sd((Intercept)) 0.749509 1.332465 2.368835

 Within-group standard error:
    lower      est.     upper 
0.8292494 1.1002946 1.4599324 
> plot(fm1)
> qqnorm(resid(fm1))
> 
> #for residuals vs fitted values by Subject, also see plot.lme
> plot(fm1, resid(.)~fitted(.) | Subject, abline=0, lty=2)
> dev.off()
null device 
          1 
> 
> 
> 
> proc.time()
   user  system elapsed 
  5.940   0.153   6.151