# Rosner, chapter 12 # ================================= # 12.2 one-way anova -- fixed effects model group <- 1:6 name <-c("NS", "PS", "NI", "LS", "MS", "HS") FEF.mean <- c(3.78, 3.30, 3.32, 3.23, 2.73, 2.59) FEF.sd <- c(0.79, 0.77, 0.86, 0.78, 0.81, 0.82) ns <- c(200, 200, 50, 200, 200, 200) (smoking <- cbind(group, FEF.mean, FEF.sd, ns)) plot(group, FEF.mean, pch=19, col="dark red", main="FEF smoking data") # ================================= # 12.3 hypothesis testing in one-way anova -- fixed effects model # ================================= # Example 12.3 pulmonary disease ys <- FEF.mean ss <- FEF.sd (SS.between <- sum(ns * ys^2) - sum(ns * ys)^2 / sum(ns)) (SS.within <- sum((ns - 1) * ss^2)) n <- sum(ns) k <- 6 MS.between <- SS.between / (k - 1) MS.within <- SS.within / (n - k) # ================================= # Example 12.4 pulmonary disease # f test for one-way anova f.test.for.one.way.anova <- function(ys, ss, ns){ call <- match.call() SS.between <- sum(ns * ys^2) - sum(ns * ys)^2 / sum(ns) SS.within <- sum((ns - 1) * ss^2) n <- sum(ns) k <- length(ys) MS.between <- SS.between / (k - 1) MS.within <- SS.within / (n - k) f <- MS.between / MS.within p <- pf(f, df1=k-1, df2=n-k) return(list(SS.between=SS.between, SS.within=SS.within, MS.between=MS.between, MS.within=MS.within, n=n, k=k, f=f, p=p, call=call)) } # example ys <- FEF.mean ss <- FEF.sd ns f.test.for.one.way.anova(ys, ss, ns) # => \$SS.between # [1] 184.3762 # # \$SS.within # [1] 663.8665 # # \$MS.between # [1] 36.87524 # # \$MS.within # [1] 0.6358875 # # \$n # [1] 1050 # # \$k # [1] 6 # # \$f # [1] 57.9902 # # \$p # [1] 1 # # \$call # f.test.for.one.way.anova(ys = ys, ss = ss, ns = ns) # ================================= # Example 12.5 and 12.6 pulmonary disease # scatter plot with error bars # http://rss.acs.unt.edu/Rdoc/library/sfsmisc/html/errbar.html library(sfsmisc) k <- 6 x <- 1:k y <- FEF.mean ss <- FEF.sd SS.between <- sum(ns * ys^2) - sum(ns * ys)^2 / sum(ns) MS.between <- SS.between / (k - 1) s.squared <- MS.within se <- sqrt(s.squared / ns) yplus <- y + se yminus <- y - se errbar(x, y, yplus, yminus, pch=19, col="steelblue1", xlab="group", ylab="FEF (L/s)", main="FEF means for six smoking groups") # ================================= # t test for comparison of pairs of groups in one-way anova (LSD procedure) # Least Significant Difference procedure t.test.for.comparison.of.pairs.of.groups.in.one.way.anova <- function(y1.bar, y2.bar, n1, n2, ys, ss, ns){ call <- match.call() SS.between <- sum(ns * ys^2) - sum(ns * ys)^2 / sum(ns) SS.within <- sum((ns - 1) * ss^2) n <- sum(ns) k <- length(ys) MS.between <- SS.between / (k - 1) MS.within <- SS.within / (n - k) s.squared <- MS.within t <- (y1.bar - y2.bar) / sqrt(s.squared * (1 / n1 + 1 / n2)) if (t < 0){ p <- 2 * pt(t, df=n-k) } else { p <- 2 * (1 - pt(t, df=n-k)) } return(list(t=t, p=p, call=call)) } # example smoking y1.bar <- FEF.mean[1] y2.bar <- FEF.mean[2] n1 <- ns[1] n2 <- ns[2] ss <- FEF.sd # ns == ns t.test.for.comparison.of.pairs.of.groups.in.one.way.anova(y1.bar, y2.bar, n1, n2, ys, ss, ns) # => \$t # [1] 6.019371 # # \$p # [1] 2.421314e-09 # # \$call # t.test.for.comparison.of.pairs.of.groups.in.one.way.anova(y1.bar = y1.bar, # y2.bar = y2.bar, n1 = n1, n2 = n2, ys=ys, ss = ss, ns = ns) # ================================= # do all pairwise comparison t tests for (i in 1:5){ for (j in (i+1):6){ print(t.test.for.comparison.of.pairs.of.groups.in.one.way.anova( y1.bar=FEF.mean[i], y2.bar=FEF.mean[j], n1=ns[i], n2=ns[j], ys, ss, ns)); print("==============") } } # ================================= # Example 12.10 pulmonary disease # linear contrasts # t test for linear contrasts in one-way anova t.test.for.linear.contrasts.in.one.way.anova <- function(ys, ss, ns, cs){ call <- match.call() SS.between <- sum(ns * ys^2) - sum(ns * ys)^2 / sum(ns) SS.within <- sum((ns - 1) * ss^2) n <- sum(ns) k <- length(ys) MS.between <- SS.between / (k - 1) MS.within <- SS.within / (n - k) s.squared <- MS.within linear.contrast <- sum(cs * ys) t <- linear.contrast / sqrt(s.squared * (sum(cs^2 / ns))) if (t < 0){ p <- 2 * pt(t, df=n-k) } else { p <- 2 * (1 - pt(t, df=n-k)) } return(list(t=t, p=p, call=call)) } # example smoking ys <- FEF.mean ss <- FEF.sd # ns == ns cs <- c(1, 0, 0, -0.1, -0.7, -0.2) t.test.for.linear.contrasts.in.one.way.anova(ys, ss, ns, cs) # => \$t # [1] 14.69121 # # \$p # [1] 0 # # \$call # t.test.for.linear.contrasts.in.one.way.anova(ys = ys, ss = ss, # ns = ns, cs = cs) # ================================= # Example 12.11 pulmonary disease cs.original <- c(5.5, 25, 40) (cs.new <- cs.original - sum(cs.original) / 3) sum(cs.new) (cs <- c(c(0, 0, 0), cs.new)) ys <- FEF.mean ss <- FEF.sd # ns # cs t.test.for.linear.contrasts.in.one.way.anova(ys, ss, ns, cs) # => \$t # [1] -8.198896 # # \$p # [1] 7.064023e-16 # # \$call # t.test.for.linear.contrasts.in.one.way.anova(ys = ys, ss = ss, # ns = ns, cs = cs) # ================================= # Example 12.13 # multiple comparisons -- Bonferroni approach # comparison of pairs of groups in one-way anova -- Bonferroni multiple comparisons procedure Bonferroni.multiple.comparisons.procedure.in.one.way.anova <- function(y1.bar, y2.bar, n1, n2, ys, ss, ns, alpha){ call <- match.call() SS.between <- sum(ns * ys^2) - sum(ns * ys)^2 / sum(ns) SS.within <- sum((ns - 1) * ss^2) n <- sum(ns) k <- length(ys) MS.between <- SS.between / (k - 1) MS.within <- SS.within / (n - k) s.squared <- MS.within t <- (y1.bar - y2.bar) / sqrt(s.squared * (1 / n1 + 1 / n2)) alpha.star <- alpha / choose(k, 2) t.cv.upper <- qt(1 - alpha.star/2, df=n-k) t.cv.lower <- qt(alpha.star/2, df=n-k) if ((t.cv.lower <= t) & (t <= t.cv.upper)){ accept.H0 <- TRUE } else { accept.H0 <- FALSE } return(list(t=t, t.cv.lower=t.cv.lower, t.cv.upper=t.cv.upper, accept.H0=accept.H0, alpha.star=alpha.star, call=call)) } # example # experiment-wise type I error == .05 smoking y1.bar <- FEF.mean[1] y2.bar <- FEF.mean[2] n1 <- ns[1] n2 <- ns[2] # ys <- FEF.mean # ss <- FEF.sd # ns == ns alpha <- .05 Bonferroni.multiple.comparisons.procedure.in.one.way.anova(y1.bar, y2.bar, n1, n2, ys, ss, ns, alpha) # => \$t # [1] 6.019371 # # \$t.cv.lower # [1] -2.941972 # # \$t.cv.upper # [1] 2.941972 # # \$accept.H0 # [1] FALSE # # \$alpha.star # [1] 0.003333333 # # \$call # Bonferroni.multiple.comparisons.procedure.in.one.way.anova(y1.bar = y1.bar, # y2.bar = y2.bar, n1 = n1, n2 = n2, ys = ys, ss = ss, ns = ns, # alpha = alpha) # conclusion : the critical region using Bonferroni is (-2.941972, 2.941972) ys.sorted <- sort(ys) e <- .01 stripchart(ys.sorted, pch=19, col="dark red", xlab="FEF level", main="Bonferoni multiple comparisons procedure") segments(ys.sorted[1] - e, 1.02, ys.sorted[2] + e, 1.02, lwd=3, col="pink") segments(ys.sorted[3] - e, 1.02, ys.sorted[5] + e, 1.02, lwd=3, col="pink") # ================================= # Example 12.15 # Sheffe's multiple-comparisons procedure for linear contrasts Sheffe.multiple.comparisons.procedure.for.linear.contrasts <- function(ys, ss, ns, cs){ call <- match.call() SS.between <- sum(ns * ys^2) - sum(ns * ys)^2 / sum(ns) SS.within <- sum((ns - 1) * ss^2) n <- sum(ns) k <- length(ys) MS.between <- SS.between / (k - 1) MS.within <- SS.within / (n - k) s.squared <- MS.within linear.contrast <- sum(cs * ys) t <- linear.contrast / sqrt(s.squared * (sum(cs^2 / ns))) c.upper <- sqrt((k - 1) * qf(1 - alpha, df1=k-1, df2=n-k)) c.lower <- -c.upper if ((c.lower <= t) & (t <= c.upper)){ accept.H0 <- TRUE } else { accept.H0 <- FALSE } return(list(t=t, c.lower=c.lower, c.upper=c.upper, accept.H0=accept.H0, call=call)) } # example smoking ys <- FEF.mean ss <- FEF.sd # ns == ns # cs, from example 12.11 cs.original <- c(5.5, 25, 40) (cs.new <- cs.original - sum(cs.original) / 3) sum(cs.new) (cs <- c(c(0, 0, 0), cs.new)) Sheffe.multiple.comparisons.procedure.for.linear.contrasts(ys, ss, ns, cs) # => \$t # [1] -8.198896 # # \$c.lower # [1] -3.333672 # # \$c.upper # [1] 3.333672 # # \$accept.H0 # [1] FALSE # # \$call # Sheffe.multiple.comparisons.procedure.for.linear.contrasts(ys = ys, # ss = ss, ns = ns, cs = cs) # ================================= # the false-discovery rate # http://compdiag.molgen.mpg.de/docs/storeypp4.pdf # http://strimmerlab.org/notes/fdr.html # http://strimmerlab.org/software/fdrtool/ # http://www.gipsa-lab.grenoble-inp.fr/~sophie.achard/brainwaver/compute.FDR.html # ================================= # Example 12.16