## function to test extract values from the ANOVA## results are from either lm() or from aov()## formula gives the DV ~ IV## wts is the weight matrix for the comparisonsdunnsTest <- function(results,wts,which) { ## Let's do error checking for the function ## error <- 0 ## First, check to make sure that a formula is passed if it is a multifactorial experiment if(length(results$contrasts) > 1 && missing(which)) { print("ERROR: you need to define the specific IV for a multifactorial experiment", quote = FALSE) error <- 1 } else { ## assign the IV and the DV from either the formula or directly from the model for a single factor design if (!missing(which)) { IV <- which if (!IV %in% names(ANOVA_test$xlevels)){ print("ERROR: Independent Variable is undefined in the model", quote = FALSE) error <- 1 } } else { IV <- as.character(names(results$xlevels[1])) } DV <- as.character(results$terms[[2]]) # check to make sure that the comparisons (weights) are valid grp_level <- length(results$xlevels[[IV]]) # get the number of levels for the factor num_comp <- nrow(wts) # get the number of comparisons being made if (grp_level == 2) {print("WARNING: Dunn's test is unnecessary when k = 2; report with caution", quote = FALSE)} if (num_comp > grp_level-1) {error=1; print(paste("ERROR: number of comparisons must be <=",grp_level-1), quote = FALSE)} if (ncol(wts) != grp_level) {error=1; print(paste("ERROR: number of weights must equal", grp_level ), quote = FALSE)} for (i in 1:num_comp) {if (sum(wts[i,]) != 0){error=1; print(paste("ERROR: sum of weights for comparison",i,"must equal 0"), quote = FALSE)}} } # if no errors, then continue on with the function if (error == 0) { grp_means <- tapply(results$model[[DV]],results$model[[IV]],mean) # Get the group means from the model data grp_n <- tapply(results$model[[DV]],results$model[[IV]],length) # Get the group n from the model data # get the MSe and dfe value from the results dfe <- results$df.residual MSe <- sum(results$residuals^2)/dfe ### compute the difference values and the Dunn's t ### # function to compute the adjusted p based on multiple comparisons with Sidak correction dunns.adj <- function(x,df,c) {1-(1-((abs(1-pt(abs(x),df)))*2))^c } dunn.value <- NULL # Difference of Means dunn.t <- NULL # Computed t-value for Dunn's test using adjusted MSe term dunn.p <- NULL # Adjusted p-value for (row in 1:num_comp) { dunn.value[row] <- sum(wts[row,]*grp_means) # calculate the comparison value ## need to compute the harmonic mean for n — in case of unequal n comp_n <- NULL; for (i in 1:ncol(wts)) {if(wts[row,i] != 0) {comp_n <- c(comp_n,grp_n[i])}} # get the n in each group dunn.n <- 1/mean(1/comp_n) # calculates harmonic mean of the group ns in analysis dunn.t[row] <- (dunn.value[row])/sqrt(sum(wts[row,]^2)*MSe/dunn.n) # computed Dunn's t-score dunn.p[row] <- dunns.adj(dunn.t[row],dfe,num_comp) # compute adjusted p } group_comp <- matrix(NA,nrow=num_comp,ncol=3) # set up results for output comp_names <- matrix(NA,nrow=num_comp,ncol=1) # get the comparison names for (row in 1:num_comp) { negcomp <- NULL # set up variable for negative comparisons poscomp <- NULL # set up variable for positive comparisons for (i in 1:ncol(wts)) { if (wts[row,i] < 0) {negcomp <- cbind (negcomp,names(grp_means[i]))} else if (wts[row,i] > 0) {poscomp <- cbind(poscomp, names(grp_means[i]))}} n_len <- length(negcomp[1,]) p_len <- length(poscomp[1,]) if(n_len > 1 && p_len >1 ){ comp_names[row] <- paste(paste("(",paste(poscomp[1,],collapse="+"),")",sep=""),"-",paste("(",paste(negcomp[1,],collapse="+"),")",sep=""),sep="") } else if (n_len == 1 && p_len > 1) { comp_names[row] <- paste(paste("(",paste(poscomp[1,],collapse="+"),")",sep=""),"-",negcomp[1,],sep="") } else if (n_len > 1 && p_len == 1){ comp_names[row] <- paste(poscomp[1,],"-",paste("(",paste(negcomp[1,],collapse="+"),")",sep=""),sep="") } else { comp_names[row] <- paste(poscomp[1,],"-",negcomp[1,],sep="") } group_comp[row,] <- c(round(dunn.value[row],4),round(dunn.t[row],4),round(dunn.p[row],4)) } rownames(group_comp) = comp_names colnames(group_comp) = c("Difference", "Dunns t-Value", "Adjusted P-Value") return(as.table(group_comp)) } }