Re: [R] post hoc in repeated measures of anova

From: Dieter Menne <dieter.menne_at_menne-biomed.de>
Date: Fri, 21 Dec 2007 20:59:05 +0000 (UTC)

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
}



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Fri 21 Dec 2007 - 21:03:06 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 21 Dec 2007 - 22:30:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.