###################################################################### # Computer notes Statistics for Laboratory Scientists (2) # Lecture 11 Johns Hopkins University ###################################################################### # These notes provide R code related to the course lectures, to # further illustrate the topics in the lectures and to assist the # students to learn R. # # Lines beginning with the symbol '#' are comments in R. All other # lines contain code. # # In R for Windows, you may wish to open this file from the menu bar # (File:Display file); you can then easily copy commands into the # command window. (Use the mouse to highlight one or more lines; then # right-click and select "Paste to console".) ###################################################################### ############################## # The anova() function in R doesn't # do quite what we want for nested ANOVA; # it compares each MS to the residual # (within-subgroup) MS. # # The following function does what we want. ############################## nested.anova <- function(aov.output) { anova.output <- anova(aov.output) nr <- nrow(anova.output) # number of rows in ANOVA table ms <- anova.output[,3] # mean squares df <- anova.output[,1] # degrees of freedom F <- ms[-nr] / ms[-1] # new F ratios pval <- 1-pf(F, df[-nr], df[-1]) # new p-values anova.output[-nr,4] <- F anova.output[-nr,5] <- pval anova.output } ############################## # Mosquito example: # data at http://www.biostat.jhsph.edu/~kbroman/teaching/labstat/fourth/mosq.csv ############################## # read data; convert "individual" column to factor mosq <- read.csv("mosq.csv") mosq$individual <- factor(mosq$individual) # the nested anova; we use the symbol "/" in the aov() function # to indicate the nesting mosq.aov <- aov(length ~ cage / individual, data=mosq) nested.anova(mosq.aov) # what happens when we use the average of the two measurements # for each mosquito? aves <- tapply(mosq$length, list(mosq$individual, mosq$cage), mean) aves <- as.numeric(aves) # turn it into a vector (rather than a matrix) cage <- factor(rep(c("A","B","C"), rep(4,3))) anova(aov(aves ~ cage)) # what happens when we ignore the cages and just do # anova with individual mosquitos as the groups? ind <- factor( paste(as.character(mosq$cage), as.character(mosq$ind), sep=":") ) anova(aov(mosq$length ~ ind)) # what happens when we ignore the individual mosquitoes and just # do ANOVA with the cages as the groups? anova(aov(length ~ cage, data=mosq)) ############################## # 2nd example: flies # data at http://www.biostat.jhsph.edu/~kbroman/teaching/labstat/fourth/flies.csv ############################## flies <- read.csv("flies.csv") nested.anova(aov(rsp ~ strain / jar, data=flies)) ################## # End of comp11.R ##################