[S] ANOVAs and MANOVAs for repeated measures (solved examples)

Gabriel Baud-Bovy (baudbovy@fpshp1.unige.ch)
Tue, 17 Feb 98 22:29:56 +0100


Dear all,

First, I would like to than those of you who responded to my previous
postings (thompson@HSB-PPC.baylor.edu, Rick Bilonick and Joseph Lucke).
In particular, Joseph Lucke gave me useful pointers to analyze the
doubly multivariate design (see also the last example below):

> Your response matrix Y=cbind(wtloss,esteem) should be a n x (p1+p2) matrix,
> where n is the number of subjects and p1 is the number of measurement
> episodies (months, I suspect) for wtloss and p2 is the number of esteem
> measuremets. I hope p1=p2 as this will make life a little easier.
>
> We're going to be forced to construct the design matrix on Y, call it T.
> This will be a little complicated as T will be a block diagonal matrix
> (doubly multivariate) and depends on the details of your
> repeated-measrues design
>
> Then we will use manova(Y %*%T ~ 1) for the analysis. The 1 on the right
> means only the intercept is a between subjects "factor". All the important
> information is in T.

In this summary, I replicated some examples of univariate and multivariate
analysis for repeated measures that I found in various statistical books.

Disclaimer: The comments may be wrong (I am not an expert in statistics)
but the results are those that I found in the books.

Now, I am working with my own data and I still find difficult to setup
the data frame and the design matrices T on Y to test all the desired
effects (I have a doubly multivariate deisgn with 2 dependent variables
and 2 within factors). Two responders suggested that I should do these
analyses by matrix computation. If anyone has some function for it
or a comment on this posting, I would be very interested.

Finally, regarding my question on the F adjustements in univariate
analysis of reapeated measure, I received the answer of John James
from MathSoft support:

>I have asked our
>technical specialist in this area, and he confirms that s-plus can not
>currently compute the Huynh-Feldt or the Greenhouse-Geisser adjustments in the
>case of repeated measures analysis with an ANOVA. Until the theory of and
>software for linear mixed effects models was available, the H-F and G-G
>adjustments were standard.
>
>However, we have implemented a more appropriate functionality - lme(), which
>deals directly with the case. I expect the s-news will compare the various
>options much more succinctly than I can.

and thompson@HSB-PPC.baylor.edu indicated to me that:

> For step-by-step illustrations of computations of some "F-adjustments"
> you mention, see Kirk, R.E. (1995) Experimental Design, Brooks/Cole
> (e.g., chap. 7 p.280ff). I also recall from earlier editions of
> Tabachnik and Fidell that they give you the matrix computations for
> all analyses.

I am planning to write a function for the H-F F-adjustement
in some simple cases but if anyone has already a function for it,
I would be interested.

Gabriel Baud-Bovy

#
# Examples of univariate and multivariate analysis of variance
# for repeated measure designs:
#
# 1) Repeated Measure Design (Hays, Table 13.21.2, p. 518)
# (1 dependent variable, 2 independent variables: 0 between, 2 within)
# 2) Repeated measures design (Stevens, chapt 13.2, p. 442)
# (1 dependent variable, 1 independent variable: 0 between, 1 within)
# 3) Repeated measure design (Stevens, chapt 13.12, p. 468)
# (1 dependent variable, 3 independent variables: 1 beteween, 2 within).
# 4) Doubly multivariate design (Tabachnick and Fidell, chapt 10.5.3, p. 476)
# (2 dependent variables, 2 independent variables: 1 between, 1 within)
#
# References:
#
# Hays (1988, 4th ed.) "Statistics"
# Stevens (1992, 2nd ed) "Applied Multivariate Statistics for the Social
Sciences"
# Tabachnick and Fidell (1996, 3rd ed.) "Using Multivariate Statistics"

===================================
# Example: Mixed effects model (Hays, Table 13.21.2, p. 518)
#
# Repeated Mesure Design:
# - one random variable: subject
# - 2 within-subject fixed effects: shape (2 levels), color (2 levels)
# - 1 replication
#
fname<-list(shape=1:2,color=1:2,su=1:12)
data<-fac.design(c(2,2,12),factor.names=fname,replications=1)
data<-data[order(data$shape,data$color,data$su),]
data$obs<-c(
49,47,46,47,48,47,41,46,43,47,46,45,
48,46,47,45,49,44,44,45,42,45,45,40,
49,46,47,45,49,45,41,43,44,46,45,40,
45,43,44,45,48,46,40,45,40,45,47,40)
data<-data[order(data$su,data$shape,data$color),]
---------------------
# repeated measure design treated as a randomized design
summary(aov(obs~(shape+color)*su,data))
Df Sum of Sq Mean Sq F Value Pr(F)
shape 1 12.0 12.00000 4.721311 0.0505358
color 1 12.0 12.00000 4.721311 0.0505358
su 11 226.5 20.59091 8.101341 0.0005455
shape:su 11 9.5 0.86364 0.339791 0.9580180
color:su 11 17.5 1.59091 0.625931 0.7769581
Residuals 12 30.5 2.54167
Note: F are not correct for a repeated design
-------------------
summary(aov(obs~shape+color+Error(su/(shape+color)),data))
Error: su
Df Sum of Sq Mean Sq F Value Pr(F)
Residuals 11 226.5 20.59091

Error: shape %in% su
Df Sum of Sq Mean Sq F Value Pr(F)
shape 1 12.0 12.00000 13.89474 0.003337568
Residuals 11 9.5 0.86364

Error: color %in% su
Df Sum of Sq Mean Sq F Value Pr(F)
color 1 12.0 12.00000 7.542857 0.01901175
Residuals 11 17.5 1.59091

Error: Within
Df Sum of Sq Mean Sq F Value Pr(F)
Residuals 12 30.5 2.541667

note: F(shape %in% su) = MS(shape)/MS(shape:su) = 13.89474
F(color %in% su) = MS(color)/MS(color:su) = 7.542857
--------------------
summary(aov(obs~shape:color+Error(su/(shape:color)),data))
Error: su
Df Sum of Sq Mean Sq F Value Pr(F)
Residuals 11 226.5 20.59091

Error: shape:color %in% su
Df Sum of Sq Mean Sq F Value Pr(F)
shape:color 3 24.0 8.000000 4.591304 0.008570969
Residuals 33 57.5 1.742424

note: F = MS(shape:color)/MS(shape:color:su) = 4.59 is correct

-------------------
# Warning: the documentation for SPLUS suggests that the following
# analzsis is correct (see Guide to Statistical Analysis, S+ 3.2,
# chapter 11.5) but it is not the case :
summary(aov(obs~shape+color+Error(su),data))
Error: su
Df Sum of Sq Mean Sq F Value Pr(F)
Residuals 11 226.5 20.59091
Error: Within
Df Sum of Sq Mean Sq F Value Pr(F)
shape 1 12.0 12.00000 7.095652 0.01172588
color 1 12.0 12.00000 7.095652 0.01172588
Residuals 34 57.5 1.69118
Note: the F are not correct for a repetead measures design
==================================
# Example: Single group repeated measures (Stevens, 13.2, p.442)
#
# Design:
# - one random variable (subject)
# - 1 fixed effect (drug)
# - 1 replication
-------------
# Data frame
fname<-list(drug=1:4,su=1:5)
data<-fac.design(c(4,5),factor.names=fname,replications=1)
data<-data[order(data$drug,data$su),]
data$obs<-c(
30,14,24,38,26,
28,18,20,34,28,
16,10,18,20,14,
34,22,30,44,30)
data<-data[order(data$su,data$drug),]
--------------
# completely randomized analysis:
summary(aov(obs~drug,data=data))
Df Sum of Sq Mean Sq F Value Pr(F)
drug 3 698.2 232.7333 4.692204 0.01552821
Residuals 16 793.6 49.6000
--------------
# repeated measures analysis:
summary(aov(obs~drug + Error(su/drug),data=data))
Error: su
Df Sum of Sq Mean Sq F Value Pr(F)
Residuals 4 680.8 170.2

Error: drug %in% su
Df Sum of Sq Mean Sq F Value Pr(F)
drug 3 698.2 232.7333 24.75887 1.992501e-05
Residuals 12 112.8 9.4000
--------------
# reformat data frame for multivariate analysis
data.maov<-as.data.frame(cbind(
unique(data$su),
t(sapply(split(data$obs,data$su),rbind))))
dimnames(data.maov)[[2]] <-
c("su",
outer("drug",levels(data$drug),paste,sep=""))
-----------------
# test drug effect
tmp<-manova(cbind(drug1,drug2,drug3,drug4)%*%contr.poly(4)~1,data=data.maov)
summary(tmp,intercept=T,test="pillai")
Df Pillai Trace approx. F num df den df P-value
(Intercept) 1 0.97707 28.41231 3 2 0.03419
Residuals 4
summary(tmp,intercept=T,test="hotelling-lawley")
Df Hotelling-Lawley approx. F num df den df P-value
(Intercept) 1 42.61846 28.41231 3 2 0.03419
Residuals 4
summary(tmp,intercept=T,test="wilk")
Df Wilks Lambda approx. F num df den df P-value
(Intercept) 1 0.02293 28.41231 3 2 0.03419
Residuals 4
--------------
# note: instead of contro.poly(4), one can use for example
contr<-matrix(c(
1, 0, 0,
0, 1, 0,
0, 0, 1,
-1,-1,-1),nrow=4,byrow=T)
#======================================
# Example: Repeated measure design (one beteween, two within).
# Stevens (1992, 2nd ed) chapt 13.12 p. 468
# Experiment on the effectiveness of drug treatement
# Design:
# one dependent variable
# one between-subject factor: group (two groups of 8 subjects)
# two within-subject factors: drug (2 levels), dose (3 levels)
-----------
# data frame (dij: effect of drug i at dose j)
data<-matrix(c(
19,22,28,16,26,22,
11,19,30,12,18,28,
20,24,24,24,22,29,
21,25,25,15,10,26,
18,24,29,19,26,28,
17,23,28,15,23,22,
20,23,23,26,21,28,
14,20,29,25,29,29,
16,20,24,30,34,36,
26,26,26,24,30,32,
22,27,23,33,36,45,
16,18,29,27,26,34,
19,21,20,22,22,21,
20,25,25,29,29,33,
21,22,23,27,26,35,
17,20,22,23,26,28),nrow=16,byrow=T)
data<-cbind.data.frame(su=1:16,gp=rep(1:2,rep(8,2)), data)
data$gp<-factor(data$gp)
dimnames(data)[[2]]<-c("su","gp","d11","d12","d13","d21","d22","d23")
------------
# test group main effect
contr<-matrix(c(
1,
1,
1,
1,
1,
1),nrow=6,byrow=T)
tmp<-aov(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp , data) # aov
summary(tmp,intercept=T)
Df Sum of Sq Mean Sq F Value Pr(F)
(Intercept) 1 333795.1 333795.1 1461.324 0.00000000
gp 1 1620.1 1620.1 7.092 0.01854578 # group main effect
Residuals 14 3197.9 228.4
------------
# test main effect drug and drug:group interaction
contr<-matrix(c(
1,
1,
1,
-1,
-1,
-1),nrow=6,byrow=T)
tmp<-aov(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp , data) # aov
summary(tmp,intercept=T)
Df Sum of Sq Mean Sq F Value Pr(F)
(Intercept) 1 2093.062 2093.062 13.00111 0.002865858 # drug main effect
gp 1 1958.063 1958.063 12.16255 0.003624241 # drug:group interaction
Residuals 14 2253.875 160.991
-------------
# test main effect dose and dose:group interaction
contr<-matrix(c(
1, 0,
0, 1,
-1,-1,
1, 0,
0, 1,
-1,-1),nrow=6,byrow=T)
tmp<-manova(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp , data) # manova
summary(tmp,test="wilk",intercept=T)
Df Wilks Lambda approx. F num df den df P-value
(Intercept) 1 0.20466 25.26075 2 13 3e-05 # dose
main effect
gp 1 0.81738 1.45223 2 13 0.26962 #
dose:group main effect
Residuals 14
------------
# test drug:dose and drug:dose:group interaction
contr<-matrix(c(
1,-1,
0, 2,
-1,-1,
-1, 1,
0,-2,
1, 1,),nrow=6,byrow=T)
tmp<-manova(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp , data) # manova
summary(tmp,test="wilk",intercept=T)
Df Wilks Lambda approx. F num df den df P-value
(Intercept) 1 0.87396 0.93739 2 13 0.41658 #
drug:dose interaction
gp 1 0.85686 1.08583 2 13 0.36637 #
drug:dose:group
Residuals 14
==================================================
# Example: Doubly multivariate design (Tabachnick and Fidell, Chapt 10.5.3,
p. 476)
#
# Design:
# two dependent variables: wtloss, esteem
# one between variables: cbw
# two within variables: month
-----------------------------------------
# data frame
d1_ c(4,4,4,3,5,6,6,5,5,3,4,5, 6,5,7,6,3,5,4,4,6,7,4,7,
8,3,7,4,9,2,3,6,6,9,7,8)
d2_ c(3,4,3,2,3,5,5,4,4,3,2,2, 3,4,6,4,2,5,3,2,5,6,3,4,
4,6,7,7,7,4,5,5,6,5,9,6)
d3_ c(3,3,1,1,2,4,4,1,1,2,2,1, 2,1,3,2,1,4,1,1,3,4,2,3,
2,3,4,1,3,1,1,2,3,2,4,1)
d4_10+c(4,3,7,1,6,7,7,3,4,4,6,5, 2,3,7,6,6,3,2,2,7,9,5,6,
6,9,5,6,3,6,3,5,5,6,6,7)
d5_10+c(3,4,2,1,5,8,6,5,4,5,6,3, 1,4,1,5,7,1,1,1,6,9,5,4,
2,9,1,2,2,3,3,2,3,4,6,7)
d6_10+c(5,7,6,2,4,8,9,5,5,3,1,6, 4,5,8,8,5,5,4,1,9,9,5,8,
6,6,9,8,7,7,6,8,8,7,9,7)
tabach_cbind.data.frame( su=rep( paste("s",as.character(1:36), sep=""), 3),
cbw=rep( rep( c("ctl","diet","diex"),rep(12,3) ), 3 ),
month=rep(1:3,rep(36,3)),
wtloss=c(d1,d2,d3), esteem=c(d4,d5,d6) )
tabach$month_factor(tabach$month)
-------------------------
# Singly multivariate analysis
# test month main effect and month:cbw interaction (assume homogeneity of
covariance)
tmp<-manova(cbind(wtloss,esteem)~month*cbw+ Error(su/month),data=tabach)
summary(tmp,test="wilk",intercept=T)
Error: su
Df Wilks Lambda approx. F num df den df P-value
cbw 2 0.77301 2.19821 4 64 0.07907
Residuals 33
Error: month %in% su
Df Wilks Lambda approx. F num df den df P-value
month 2 0.1988 40.3992 4 130 0 # month effect
month:cbw 4 0.6791 3.4687 8 130 0.0012 # month:cbw
interaction
Residuals 66
----------------------
# Data Format for doubly multivariate analysis
# Between-Subject Dependent variables (esteem and wtloss) for each
level (3)
# Variable (cbw) of the within-Subject variable (month)
# cbw esteem1 esteem2 esteem3 wtloss1 wtloss2 wtloss3
# S1 . . . .
# S2
# ..
# S36
#
t1<-sapply(split(tabach$wtloss,tabach$month),rbind)
t2<-sapply(split(tabach$esteem,tabach$month),rbind)
tabach2<-cbind.data.frame(su=paste("s",as.character(1:36), sep=""),
cbw=rep( c("ctl","diet","diex"),rep(12,3)),
t1,t2)
names(tabach2)<-c("su","cbw","wl1","wl2","wl3","e1","e2","e3")
-------------------
# doubly multivariate analysis
# test month main effect and month:cbw interaction
# design matrix on the response matrix (block diagonal matrix)
contr<-matrix(0,nrow=6,ncol=4)
contr[1:3,1:2]<-contr.poly(3)
contr[4:6,3:4]<-contr.poly(3)
tmp<-manova(cbind(wl1,wl2,wl3,e1,e2,e3) %*% contr ~ cbw , data=tabach2,)
summary(tmp,test="wilk",intercept=T)
Df Wilks Lambda approx. F num df den df P-value
(Intercept) 1 0.11239 59.2305 4 30 0 # month
main effect
cbw 2 0.39542 4.42706 8 60 0.0003 # month:cbw
interaction
Residuals 33

-------------------------------------------------------------
Gabriel Baud-Bovy baudbovy@fpshp1.unige.ch
Université de Genève, FAPSE tel. +41 22 705 97 67
9, route de Drize fax +41 22 300 14 82
1227 Carouge, Switzerland home tel. +41 22 320 21 38
-------------------------------------------------------------

-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news