[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