[R-sig-ME] Help with mixed models function
Cody Likavec
codylikavec at gmail.com
Tue Aug 13 22:01:43 CEST 2013
Hello,
I'm trying to write a function to automate a mixed-model ANOVA (one-between,
one-within). I've gotten most of it to work. I can throw in my data and get
output for the within and between subjects effects. However, when I try to
do the pairwise comparisons, I'm running into a wall on how to actually run
them, since they're embedded in the function. I get stuck once I get to
"pair2". Because the "iv" is going to dynamically change, I can't seem to
figure out what to use to give me output that will dynamically change based
on what I enter into the function. Below is what I have for the full
function.
onewithin.onebetween <- function(dvs, iv, id, data){
dataframe <- cbind(data[id], data[dvs], data[iv])
means <- sapply(split(dataframe[dvs], data[iv]), colMeans,
na.rm=TRUE)
sds <- sapply(split(dataframe[dvs], data[iv]), sd, na.rm=TRUE)
shapiros <- sapply(data[dvs], shapiro.test);
myexpr.df <- make.rm(constant=c(id, iv), repeated=c(dvs),
data=dataframe);
formula <- as.formula(paste("repdat~contrasts +",iv, "+
contrasts*",iv));
random <- as.formula(paste("~1|",id))
owobanova <- lme(formula, data=myexpr.df, random=random);
anova <- anova(owobanova);
#pairwise1 <- summary(glht(owobanova , linfct = mcp(contrasts =
"Tukey")));
#pair2 <- paste(iv,'= "Tukey" ');
#pairwise2 <- summary(glht(owobanova , linfct = mcp(pair2)));
#pair3 <- paste(contrasts,"*",iv,'= "Tukey" ');
#pairwise3 < summary(glht(owobanova , linfct = mcp(pair3)));
output <- list();
output$means <- as.list(as.data.frame(means));
output$sds <- as.list(as.data.frame(sds));
output$shapiros <- as.list(as.data.frame(shapiros));
output$anova2 <- anova;
#output$pairwise <- as.list(unclass(pairwise1$test$coefficients));
#output$pairwise2 <- as.list(unclass(pairwise1$test$pvalues));
#output$pairwise3 <- as.list(unclass(pairwise2$test$coefficients));
#output$pairwise4 <- as.list(unclass(pairwise2$test$pvalues));
#output$pairwise5 <- as.list(unclass(pairwise3$test$coefficients));
#output$pairwise6 <- as.list(unclass(pairwise3$test$pvalues));
return(output)
}
More information about the R-sig-mixed-models
mailing list