Dear Helios, I've now had a chance to look at your code for the factorltest.mlm() function. I agree that the function makes it easier to test hypotheses in repeated-measures ANOVAs. When I have some more time, I'll make a few suggestions (off list) for improving the user interface to the function.
Best, John -------------------------------- John Fox Senator William McMaster Professor of Social Statistics Department of Sociology McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox > -----Original Message----- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On > Behalf Of Helios de Rosario > Sent: September-22-11 12:41 PM > To: r-help@r-project.org > Subject: [R] Wrapper of linearHypothesis (car) for post-hoc of repeated > measures ANOVA > > For some time I have been looking for a convenient way of performing post-hoc > analysis to Repeated Measures ANOVA, that would be acceptable if sphericity > is violated (i.e. leaving aside post-hoc to lme models). > > The best solution I found was John Fox's proposal to similar requests in R- > help: > http://tolstoy.newcastle.edu.au/R/e2/help/07/09/26518.html > http://tolstoy.newcastle.edu.au/R/e10/help/10/04/1663.html > > However, I think that using linearHypothesis() is not as straightforward as I > would desire for testing specific contrasts between factor levels. The > hypotheses must be defined as linear combinations of the model coefficients > (subject to response transformations according to the intra-subjects design), > which may need some involved calculations if one is thinking on differences > between "this and that" factor levels (either between-subjects or intra- > subjects), and the issue gets worse for post-hoc tests on interaction > effects. > > For that reason, I have spent some time in writing a wrapper to > linearHypothesis() that might be helpful in those situations. I copy the > commented code at the end of this message, because although I have > successfully used it in some cases, I would like more knowledgeable people > to put it to test (and eventually help me create a worthwile contribution for > other people that could find it useful). > > This function (which I have called "factorltest.mlm") needs the multivariate > linear model and the intrasubject-related arguments (idata, > idesign...) that would be passed on to Anova() (from car) for a repeated > measures analysis, or directly the Anova.mlm object returned by Anova() > instead of idata, idesign... (I have tried to explain it clearly in the > commentaries to the code.) > > Moreover, it needs an argument "levelcomb": a list that represents the level > combinations of factors to be tested. There are different ways of > representing those combinations (through names of factor levels, or > coefficient vectors/matrices), and depending on the elements of that list the > test is made for main effects, simple effects, interaction contrasts, etc. > > For instance, let me use an example with the OBrienKaiser data set (as in the > help documentation for Anova() and linearHypothesis()). > > The calculation of the multivariate linear model and Anova is copied from > those help files: > > > phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, > 5)), > + levels=c("pretest", "posttest", "followup")) > > hour <- ordered(rep(1:5, 3)) > > idata <- data.frame(phase, hour) > > mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, > + post.1, post.2, post.3, post.4, post.5, > + fup.1, fup.2, fup.3, fup.4, fup.5) ~ > treatment*gender, > + data=OBrienKaiser) > > av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour) > > Then, let's suppose that I want to test pairwise comparisons for the > significant main effect "treatment" (whose levels are c("control","A","B")). > For the specific contrast between treatment "A" > and the "control" group I can define "levelcomb" in the following > (equivalent) ways: > > > levelcomb <- list(treatment=c("A","control")) levelcomb <- > > list(treatment=c(A=1,control=-1)) levelcomb <- > > list(treatment=c(-1,1,0)) > > Now, let's suppose that I am interested in the (marginally) significant > interaction between treatment and phase. First I test the simple main effect > of phase for different levels of treament (e.g. for the "control" > group). To do this, levelcomb must have one variable for each interacting > factor (treatment and phase): levelcomb$treatment will specify the treatment > that I want to fix for the simple main effects test ("control"), and > levelcomb$phase will have a NA value to represent that I want to test all > orthogonal contrasts within that factor: > > > levelcomb <- list(treatment="control",phase=NA) > > I could also use numeric vectors to define the levels of "treatment" > that I want to fix, as in the previous example, or if I want a more > complicated combination (e.g. if I want to test the phase effect for pooled > treatments "A" and "B"): > > > levelcomb <- list(treatment=c(A=1,B=1),phase=NA) > > The NA value can be replaced by the specific contrast matrix that I would > like to use. For instance, the previous statement is equivalent to the > following one: > > > levelcomb <- list(treatment=c(0,1,1),phase=contr.sum(levels(phase))) > > And then, let's see an example of interaction contrast: I want to see if the > difference between the "pre-test" phase and the following two, itself differs > between the "control" group and the two treatments. This would be > > > levelcomb <- list(treatment=c(2,-1,-1),phase=c(2,-1,-1)) > or > > levelcomb <- > list(treatment=c(control=2,A=-1,B=-1),phase=c(pretest=2,posttest=- > 1,followup=-1)) > > etc. > > Now, to perform the test I just use this "levelcomb" list object with the > model and ANOVA previously performed, in the following way: > > > factorltest.mlm(mod.ok,levelcomb,av.ok) > > Or if I want to use idata and idesign directly: > > > factorltest.mlm(mod.ok,levelcomb,idata=idata,idesign=~phase*hour) > > Of course, this function only performs one test at a time. To perform > multiple tests as suggested by John, the user should run it in a loop (and > then correct the p-values as he or she see fit). > > If you find this useful, please let me know if you also find some error, > opportunity of improvement, or any commentary. > > Thanks, > Helios > > [Now follows the full commented code:] > > # factorltest.mlm(modmlm,levelcomb,aovmlm,...) > # > # Hypothesis testing for linear contrasts of factor levels, for multivariate > # linear models with response transformation, as in a repeated-measures > ANOVA. > # > # Arguments: > # modmlm: multivariate linear model (mlm object) # levelcomb: list with the > combinations of levels for each factor, in form of # a literal expression or > numeric matrix (see Details). > # aovmlm: result from Anova to modmlm (i.e. Anova.mlm object) # ...: > additional arguments passed to Anova() --from the car package-- to # perform > an ANOVA with an intra-subject design (see Details). > # > # Details: > # An intra-subject design is required for the response transformation of the > # linear model. The intra-subject design may be defined by the arguments # > idata, idesgin (and optionally icontrasts) defined for the function > Anova() > # in the car package, or alternatively by an Anova.mlm object resulting from > # that function (applied to modmlm). If both modes of defining the intra- > subjects # design are used, the design implied by aovmlm prevales. > # > # The name of each element in levelcomb must coincide with a factor (i.e. a > main # effect) analysed in the ANOVA (as defined in aovmlm or indirectly by > modmlm and # the other arguments). The value of each element can be: (1) the > name of one level # of that factor, (2) a vector with two level names of that > factor --to perform # a paired contrast between those levels--, (3) a named > vector of coefficients # whose names coincide with some levels of the factor, > for more specific linear # combinations --unspecified coefficients are > assumed to be zero--, (4) an unnamed # vector of coefficients whose length is > equal to the number of levels in the # factor --the values of the vector are > assigned to the factor levels in the # default order--, (5) a matrix with > named or unnamed rows, in which each column # represents a combination of > levels as in (3) or (4). > # > # This function depends on linearHypothesis() in the car package. The > combinations # specified in levelcomb are tested against zero or a matrix of > zeros. > > > > factorltest.mlm <- function(modmlm,levelcomb,aovmlm,...){ > ## AUXILIARY FUNCTIONS > # checkfactors() is used to check consistency of level combinations in > the > # factors defined in levelcomb, and convert them to numeric matrices > checkfactors <- function(fframe,fnames,levelcomb){ > for (fname in fnames){ > f <- fframe[[fname]] > fc <- levelcomb[[fname]] > # Literal expressions must represent one level or a > pair of > levels > if ((is.character(fc)) & (!any(is.na(fc)))){ > fcnames <- fc > if (length(fc)==1){ > # If there is one level, > evaluate the values at that level > fc <- 1 > names(fc) <- fcnames > }else if (length(fcnames)==2){ > # If there is a pair of levels, > evaluate the contrast > fc <- c(1,-1) > names(fc) <- fcnames > }else{ > stop("Incorrect number of > literal levels for factor \"",fname,"\" (must be 1 or 2).") > } > } > # Check the numeric coefficients of the factor levels > if ((is.numeric(fc)) & (!any(is.na(fc)))){ > # Make fc into a matrix (in case it was a > vector) > fc <- as.matrix(fc) > # Unnamed coefficient matrices must have the > same rows > as levels in the factor > if (is.null(rownames(fc))){ > if (nrow(fc) != nlevels(f)){ > stop("Mismatch in the > number of unnamed factor levels for \"",fname,"\".") > }else{ > rownames(fc) <- > levels(f) > } > }else{ > # Named coefficients are > assigned to a matrix of factor levels, filled in with zeros > flevels <- > match(rownames(fc),levels(f)) > fctmp <- fc > fc <- > matrix(0,nrow=nlevels(f),ncol=ncol(fctmp)) > if (any(is.na(flevels))){ > stop("Mismatch in the > names of factor levels for \"",fname,"\".") > } > fc[flevels,] <- fctmp > rownames(fc) <- levels(f) > } > # Convert NA into default contrast matrix > }else if (is.na(fc)){ > fc <- contrasts(f) > }else{ > stop("Invalid value assigned to factor > \"",fname,"\".") > } > levelcomb[[fname]] <- fc > } > return(levelcomb) > } > > # makeT() creates a transformation matrix, used to define the Linear > Hypothesis > # matrix (for between-subjects factors) or the response transformation > matrix > # (for within-subjects factors), depending on the combinations of > factor levels. > # The transformation matrix is created progressively, by "translating" > the combinations > # of each factor into matrices that are sequentially multiplied. > makeT <- function(fframe,fnames,levelcomb){ > # First matrix, defined from the identic rows of the > between/within-subjects > # model data frame, when unspecified factors are removed. > # A dummy column with ones is added to the model data frame, to > avoid > # problems when all columns be eventually removed. > # (The name of the dummy column is coerced to be different from > any other one.) > dummyname <- paste("z",max(fnames),sep=".") > fframe[[dummyname]] <- 1 > # All factors are interacting, the transformation matrix > (Tm) in the first step is diagonal. > if (length(fnames)==ncol(fframe)-1){ > m <- nrow(fframe) > Tm <- diag(m) > }else{ > # Remove columns of unspecified factors > fframe <- fframe[,c(fnames,dummyname)] > # Vector with a different value for each combination of > factors > fframe_1 <- apply(fframe,1,paste,collapse="") > n <- nrow(fframe) > # Subset of unique elements > fframe <- unique(fframe) > fframe_1.m <- unique(fframe_1) > m <- nrow(fframe) > # Matrix with coefficients for averaging identical > combinations of factors > # Rows in Mo are the original (repeated) combinations > To <- matrix(rep(fframe_1,each=m),nrow=m) > # Columns in Mb are the unique combinations > Tu <- matrix(rep(fframe_1.m,n),nrow=m) > Tm <- (To==Tu) * m/n > } > # Number of contrasts to calculate for the current factor > (initialized as 1) > nc <- 1 > # Labels for the final transformation matrix (only defined for > multiple contrasts) > Tlabels <- NULL > # Progressive transformation of Tm, factor by factor > for (fname in fnames){ > # f is the current factor vector > f <- fframe[[fname]] > # n is the factor vector length > n <- m*nc > # Remove the corresponding column from the model data > frame > fframe[[fname]] <- NULL > # nc is the number of contrasts for the current factor > (updated) > nc <- ncol(levelcomb[[fname]]) > if (nc > 1){ > # If there are multiple contrasts, the rows of > the > model data frame are > # replicated (by the number of > contrasts), and a new column is added to > # assign a contrast to each group of copied > rows. > fframe[[ncol(fframe)+1]] <- 1 > fframe <- fframe[rep(1:n,nc),] > fframe[[ncol(fframe)]] <- > rep(1:nc,each=n) > # Moreover, we create labels to identify what > contrast > is represented > # by each row of the final > transformation matrix > newTlabels <- > paste(fname,as.character(1:nc),sep="") > if (is.null(Tlabels)){ > Tlabels <- newTlabels > }else{ > # (In case there is more than > one factor with multiple contrasts) > Tlabels <- > paste(rep(Tlabels,nc), rep(newTlabels,each=length(Tlabels)), sep=":") > } > } > # The same routine as for the first Tm: vector with > combined values > # of the (transformed) model data frame... > fframe_1 <- apply(fframe,1,paste,collapse="") > # ... subset to unique rows > fframe <- unique(fframe) > fframe_1.m <- unique(fframe_1) > m <- nrow(fframe)/nc > # And create transformation matrix, depending on the > combination of factor > # levels defined in levelcomb > kf <- t(levelcomb[[fname]]) > kf <- matrix(rep(kf[,f],each=m),ncol=n) > To <- > matrix(rep(matrix(fframe_1,ncol=n,byrow=TRUE),each=m),ncol=n) > Tu <- matrix(rep(fframe_1.m,n),ncol=n) > Tm <- ((To==Tu) * kf) %*% Tm > } > # When the loop is finished, assign labels to final Tm, and > return it > rownames(Tm) <- Tlabels > return(Tm) > } > > ## END OF AUXILIARY FUNCTIONS > ## MAIN PROCEDURE > > # 1. Check class of modmlm and intra-subject design from input > arguments > if (missing(aovmlm)){ > aovmlm <- Anova(modmlm,...) > } > > # 2. Define model data frames: > # Between-subjects model data frame, with contrasts copied from linear > model > bframe <- expand.grid(modmlm$xlevels) > bfactors <- names(bframe) > for (bf in bfactors){ > contrasts(bframe[[bf]]) <- modmlm$contrasts[[bf]] > } > # Within-subjects model data frame, with contrasts copied from intra- > subject design > wframe <- aovmlm$idata > wfactors <- names(wframe) > for (wf in wfactors){ > if (is.null(attr(wframe[[wf]], "contrasts"))){ > contrasts(wframe[[wf]]) <- if > (is.ordered(wframe[[wf]])) aovmlm$icontrasts[2] else aovmlm$icontrasts[1] > } > } > > # 3. Check that interacting factors in levelcomb are included in both > the > # between-subjects and within-subject designs > ifactors <- names(levelcomb) > iterm <- paste(ifactors,collapse=":") > itermsort <- paste(sort(ifactors),collapse=":") > anovaterms <- > lapply(strsplit(aovmlm$terms,":"),function(x){paste(sort(x),collapse=":")}) > if (all(is.na(match(anovaterms,itermsort)))){ > warning("The term \"", iterm, "\" is not in the model > design.") > } > > # 4. Search factors of levelcomb in bframe and wframe, and convert the > level > # combinations into numeric matrix format > bfactor.interact <- match(ifactors,bfactors) > bfactors <- > bfactors[bfactor.interact[!is.na(bfactor.interact)]] > levelcomb <- checkfactors(bframe,bfactors,levelcomb) > wfactor.interact <- match(ifactors,wfactors) > wfactors <- > wfactors[wfactor.interact[!is.na(wfactor.interact)]] > levelcomb <- checkfactors(wframe,wfactors,levelcomb) > > # 5. Preliminary definition of the Linear Hypothesis matrix (L) > rhf <- paste(as.character(formula(modmlm))[c(1,3)],collapse=" > ") > L <- model.matrix(as.formula(rhf), data=bframe) > > # 6. Transformed Linear Hypothesis (L) and response transformation (P) > matrices > if (length(bfactors)){ > L <- makeT(bframe,bfactors,levelcomb) %*% L > }else{ > L <- colSums(L)/nrow(L) > } > if (length(wfactors)){ > P <- t(makeT(wframe,wfactors,levelcomb)) > }else{ > P <- matrix(rep(1/nrow(wframe),nrow(wframe))) > } > > # 7. Result, consisting in: > # levelcomb: numeric matrix values of levelcomb > # lsm: table of least square means for the tested > interactions > # test: test value, from LinearHypothesis > result <- list(levelcomb=levelcomb,lsm=NULL,test=NULL) > result$lsm <- L %*% modmlm$coefficients %*% P > result$test <- linearHypothesis(modmlm,L,P=P) > return(result) > } > > INSTITUTO DE BIOMECÁNICA DE VALENCIA > Universidad Politécnica de Valencia • Edificio 9C Camino de Vera s/n • 46022 > VALENCIA (ESPAÑA) Tel. +34 96 387 91 60 • Fax +34 96 387 91 69 www.ibv.org > > Antes de imprimir este e-mail piense bien si es necesario hacerlo. > En cumplimiento de la Ley Orgánica 15/1999 reguladora de la Protección de > Datos de Carácter Personal, le informamos de que el presente mensaje contiene > información confidencial, siendo para uso exclusivo del destinatario arriba > indicado. En caso de no ser usted el destinatario del mismo le informamos que > su recepción no le autoriza a su divulgación o reproducción por cualquier > medio, debiendo destruirlo de inmediato, rogándole lo notifique al remitente. > > ______________________________________________ > R-help@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. ______________________________________________ R-help@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.
