[R] post hoc in repeated measures of anova

[Dieter Menne]

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
}