################################################# #Cohen, B. (2000). One-way independent ANOVA #In Explaining psychological statistics(2nd ed., pp.324--361) #New York, NY: John Wiley & Sons ################################################# #-----------------------------------------# # An example with three equal-sized group (pp.330-334) #-----------------------------------------# ## (a) Table 12.1 m <- c(9, 7, 5.5) #mean s <- c(3.5,3, 2.5) #standard deviation n <- 10 #each sample size (orthogonal case) sig.level <- .05 #significance level ##(b) numerater of F ratio k <- length(m) #number of groups MSb <- n * ( sum((m - mean(m))^2) / (k -1) ) #Formula 12.5 A ##(c) F value MSw <- sum(s^2) / k f.value <- MSb/MSw ##(d) significance testing dfb <- k - 1 dfw <- n * k - k qf(p=sig.level, df1=dfb, df2=dfw, lower.tail=F) #critical F p.value <- pf(q=f.value, df1=dfb, df2=dfw, lower.tail=F) #p.value ##(e) output output <- data.frame(matrix("",ncol=5, nrow=3)) rownames(output) <- c("Between", "Within", "Total") colnames(output) <- c("SS", "df", "MS", "F", "p") digits <- 2 output$SS <- c(MSb*dfb, MSw*dfw, MSb*dfb+MSw*dfw) output$df <- c(dfb, dfw, dfb+dfw) output$MS <- c(MSb, MSw,NA) output$F <- c(f.value, NA,NA) output$p <- c(p.value, NA,NA) class(output) <- c("anova", "data.frame") round(output, digits) #-----------------------------------------# #pp.330-334の関数を定義 (2次データから一元配置母数モデルの分散分析) #-----------------------------------------# oneway.anova.second <- function(m, s, n, sig.level=.05, digits=2){ ##(a) numerater of F ratio k <- length(m) #number of groups MSb <- n * ( sum((m - mean(m))^2) / (k -1) ) #Formula 12.5 A ##(b) F value MSw <- sum(s^2) / k f.value <- MSb/MSw ##(c) significance testing dfb <- k - 1 dfw <- n * k - k qf(p=sig.level, df1=dfb, df2=dfw, lower.tail=F) #critical F p.value <- pf(q=f.value, df1=dfb, df2=dfw, lower.tail=F) #p.value ##(d) output output <- data.frame(matrix("",ncol=5, nrow=3)) rownames(output) <- c("Between", "Within", "Total") colnames(output) <- c("SS", "df", "MS", "F", "p") output$SS <- c(MSb*dfb, MSw*dfw, MSb*dfb+MSw*dfw) output$df <- c(dfb, dfw, dfb+dfw) output$MS <- c(MSb, MSw,NA) output$F <- c(f.value, NA,NA) output$p <- c(p.value, NA,NA) class(output) <- c("anova", "data.frame") return(round(output, digits)) }#end function ## test oneway.anova.second(m=m, s=s, n=n) #-----------------------------------------# #Exercises A (pp.334--335) #-----------------------------------------# ##(1) s <- c(10, 15, 12, 11, 10) sum(s^2) / length(s) #MSw ##(3) dfb <- 4 dfw <- 80 (k <- 4 + 1) (n <- (dfw + k)/k ) ##(5) oneway.anova.second(m=c(28.7,34.3,37.2), s=c(11.2,12.0,13.5), n=80) ##(7) oneway.anova.second(m=c(7, 8, 5, 4), s=c(3.25,3.95,3.16,2.07), n=8) qf(p=.05, df1=3, df2=28, lower.tail=F) #critical F #-----------------------------------------# # An ANOVA example with unequal sample sizes (pp.336-340) #-----------------------------------------# ## (a) Table 12.2 m <- c(28, 30.6, 35.83) #mean s <- c(5.477, 4.393, 4.875) #standard deviation n <- c(4, 5, 6) #sample size sig.level <- .05 #significance level ##(b) Find region of rejection k <- length(m) #number of groups dfb <- k - 1 dfw <- sum(n) - k qf(p=sig.level, df1=dfb, df2=dfw, lower.tail=F) #critical F ##(c) Caluculate test statistic MSw <- sum((n-1) * s^2) /dfw #Formula 12.6 Xg <- sum(n*m)/sum(n) #grand mean MSb <- sum(n*(m-Xg)^2)/dfb #Formula 12.7 f.value <- MSb/MSw p.value <- pf(q=f.value, df1=dfb, df2=dfw, lower.tail=F) #p.value ##(d) output output <- data.frame(matrix("",ncol=5, nrow=3)) rownames(output) <- c("Between", "Within", "Total") colnames(output) <- c("SS", "df", "MS", "F", "p") digits <- 2 output$SS <- c(MSb*dfb, MSw*dfw, MSb*dfb+MSw*dfw) output$df <- c(dfb, dfw, dfb+dfw) output$MS <- c(MSb, MSw,NA) output$F <- c(f.value, NA,NA) output$p <- c(p.value,NA,NA) class(output) <- c("anova", "data.frame") round(output, digits) #-----------------------------------------# #pp.336-340の関数を定義 #-----------------------------------------# oneway.anova.second2 <- function(m, s, n, sig.level=.05, digits=2){ ##(a) Find region of rejection k <- length(m) #number of groups dfb <- k - 1 dfw <- sum(n) - k qf(p=sig.level, df1=dfb, df2=dfw, lower.tail=F) #critical F ##(b) Caluculate test statistic MSw <- sum((n-1) * s^2) /dfw #Formula 12.6 Xg <- sum(n*m)/sum(n) #grand mean MSb <- sum(n*(m-Xg)^2)/dfb #Formula 12.7 f.value <- MSb/MSw p.value <- pf(q=f.value, df1=dfb, df2=dfw, lower.tail=F) #p.value ##(c) output output <- data.frame(matrix("",ncol=5, nrow=3)) rownames(output) <- c("Between", "Within", "Total") colnames(output) <- c("SS", "df", "MS", "F", "p") output$SS <- c(MSb*dfb, MSw*dfw, MSb*dfb+MSw*dfw) output$df <- c(dfb, dfw, dfb+dfw) output$MS <- c(MSb, MSw,NA) output$F <- c(f.value, NA,NA) output$p <- c(p.value,NA,NA) class(output) <- c("anova", "data.frame") return(round(output, digits)) }#End function oneway.anova.second2(m=m, s=s, n=n) #-----------------------------------------# # Raw-score formulas (pp.340-342) #-----------------------------------------# ##(a) Table 12.2 x <- c(29,35,26,22,30,37,29,32,25,35,38,33,41,28,40) group <- factor(c(rep("A", 4), rep("B", 5), rep("C", 6))) m <- tapply(x, group, mean) s <- tapply(x, group, sd) n <- tapply(x, group, length) ##(b) SSw and MSw T <- tapply(x, group, sum) SSw <- sum(x^2)- sum(T^2/n) #Formula 12.8 k <- nlevels(group) dfw <- length(x) - k MSw <- SSw/dfw ##(c) SSb and MSb SSb <- sum(T^2/n) - sum(T)^2/sum(n) #Formula 12.9 dfb <- k - 1 MSb <- SSb/dfb ##(d) SStotal SSt <- SSb + SSw #Formula 12.10 sum(x^2) - sum(T)^2/sum(n) #Formula 12.11 #-----------------------------------------# #Exercises B (pp.348--350) #-----------------------------------------# ##(1) oneway.anova.second2(m=c(23,30,34,29,26), s=c(6.5,7.2,7,5.8,6), n=c(12,15,14,12,15)) ##(3) x <- c(86,23,47,51,63,75,42,35,56,70,46,49,28,68,52,63,82,36) group <- factor(c(rep("A", 5), rep("C", 6), rep("N", 7))) m <- tapply(x, group, mean) s <- tapply(x, group, sd) n <- tapply(x, group, length) T <- tapply(x, group, sum) SSw <- sum(x^2)- sum(T^2/n) #Formula 12.8 k <- nlevels(group) dfw <- length(x) - k MSw <- SSw/dfw SSb <- sum(T^2/n) - sum(T)^2/sum(n) #Formula 12.9 dfb <- k - 1 MSb <- SSb/dfb MSb/MSw #f ratio oneway.anova.second2(m=m, s=s, n=n) #f ratio ##(5) x <- c(3,7,1,0,9,2, 3,4,5,6,4,6, 2,0,4,6,4,1) group <- rep(factor(c("Green","Red","Blue"), levels=c("Green","Red","Blue")), each=6) m <- tapply(x, group, mean) s <- tapply(x, group, sd) n <- tapply(x, group, length) oneway.anova.second2(m=m, s=s, n=n, sig.level=.01) #f ratio qf(p=.01, df1=2, df2=15, lower.tail=F) #critical F ##(8) x <- c(8,7,7,9,4, 2,3,2,4,4, 6,8,9,5,8, 7,9,8,5,7) group <- rep(factor(c("A","B","C","D")), each=5) m <- tapply(x, group, mean) s <- tapply(x, group, sd) n <- tapply(x, group, length) oneway.anova.second2(m=m, s=s, n=n, sig.level=.05) #f ratio oneway.anova.second2(m=m, s=s, n=n, sig.level=.01) #f ratio ##(9) m <- c(14.5,12.2,17.8,15.1,13.4,12.0) n <- c(8,7,8,6,10,7) sum.sq <- 9820 k <- length(m) T <- m*n SSw <- 9820 - sum(T^2/n) #Formula 12.8 dfw <- sum(n) - k MSw <- SSw/dfw SSb <- sum(T^2/n) - sum(T)^2/sum(n) #Formula 12.9 dfb <- k - 1 MSb <- SSb/dfb MSb/MSw #f ratio qf(p=.05, df1=dfb, df2=dfw, lower.tail=F) #critical F #-----------------------------------------# # Effect size and proportion of variance accounted for(pp.351-354) #-----------------------------------------# SSb <- 161.97 SSt <- 448 k <- 3 MSw <- 23.84 eta.sq <- SSb/SSt #Formula 12.12 omega.sq <- (SSb - (k-1) * MSw)/(SSt + MSw) #Formula 12.14