[R] post hoc in repeated measures of anova
Dieter Menne
dieter.menne at menne-biomed.de
Fri Dec 21 21:59:05 CET 2007
Kazumi Maniwa <Kazumi.Maniwa <at> uni-konstanz.de> writes:
>
> Hallo, I have this dataset with repeated measures. There are two
> within-subject factors, "formant" (2 levels: 1 and 2) and "f2 Ref" (25
> levels: 670, 729, 788, 846, 905, 1080, 1100, 1120, 1140, 1170, 1480,
> 1470, 1450, 1440, 1430, 1890, 1840, 1790, 1740, 1690, 2290, 2210,
> 2120, 2040, 1950), and one between-subject factor, lang (2 levels:1
> and 2). The response variable is "thresh".
..
>
> I ran a three-way ANOVA with repeated measures:
>
> vThresh<-read.table("rRes.csv",header=T,sep=",")
> attach(vThresh)
>
vThresh.aov<-aov(thresh~factor(lang)*factor(formant)*factor(f2Ref)+
Error(factor(sub)))
> It will be great if you could help me with codes for post hoc in
> repeated measures ANOVA.
First, I suggest that you avoid attach(), and convert the numbers to factors
in the data frame. The formula and the results are much easier to read, and
the risk of making some mistake is reduced. lm and lme will work perfectly
without factor(), but the results can be misleading or simply wrong, because
numbers are treated as what one used to call continuous covariables in the
old days (still persisting in psychology departments).
vThresh$formant = as.factor(vThresh$formant)
... same for the other factors. Better even, give you factors nice names
(ok, with formants, numbers _are_ a natural order).
Then, use lm, which gives you a basic set of contrasts. Make sure that
you understand that (Intercept) stands for the base level of all factors
(e.g. 1), and the others are roughly differences, or differences of
differences (interactions).
vThresh.lm = lm(thresh~lang*formant*f2Ref+Error(factor(sub),data=vThresh)
or, even better, use lme instead of lm when you have repeated measurements:
library(nlme)
vThresh.lme = lme(thresh~lang*formant*f2Ref, data=vThresh, random=~1|subject)
You can use stepAIC from MASS to prune down your model to relevant terms.
Then, use estimable in package gmodels, or glht in package multcomp to
compute the post-hoc tests.
Getting the contrasts right can be tricky; I use an Excel Spreadsheet to
construct the contrast with a bit of conditional formatting (0 white, yellow
1) to easier spot errors, and read the named contrast ranges via ODBC.
Dieter
# some auxillary function to construct contrast in Excel and have the
# variable nicely named in estimable
"readContrasts" =
function(cname,excelfile) {
require(gdata)
require(RODBC) # for trim
channel=odbcConnectExcel(excelfile)
cn=sqlQuery(channel,paste("select * from",cname),as.is=TRUE)
odbcClose(channel)
if (class(cn) != "data.frame") stop(cn[[2]],call.=FALSE)
cn[,1] = trim(cn[,1]) # trim in gdata
cn
}
"getContrasts" =
function(cname,excelfile,rows=NULL) {
cn=readContrasts(cname,excelfile)
if (!is.null(rows)) cn = cn[rows,]
colnames= cn[,1]
rownames = gsub('#','.',colnames(cn))#[-1]
# Use _ as placeholder for an empty field
vars = do.call("rbind", strsplit(rownames,'\\.'))
vars[vars=='_'] = ''
# upper left corner must contain the names of the variables
varnames = vars[1,]
vars = vars[-1,,drop=FALSE]
rownames(vars)=rownames[-1]
colnames(vars)=varnames
cn = t(as.matrix(cn[,-1]))
rownames(cn) = rownames[-1]
colnames(cn) = colnames
attr(cn,"vars") = data.frame(vars)
cn
}
More information about the R-help
mailing list