--- title: "fruit flies" author: "Chris Parrish" date: "January 10, 2016" output: pdf_document --- fruit flies reference: - Cannon, et al., Stat2, chapter 07, example 7.1, 7.4-7.9. 7.13-7.14 - Cannon, et al., Student's R Manual, chapter 07 Import the data. {r} data <- read.csv("FruitFlies.csv", header=TRUE) head(data, 4) dim(data)  Scatterplot matrix. {r} pairs(~ Longevity + Thorax + Sleep + Treatment, data=data, col="darkred")  Levene's test. {r} library(car) leveneTest(Longevity ~ Treatment, data=data)  Multiple comparisons: Fisher's LSD {r} with(data, pairwise.t.test(Longevity, Treatment, p.adj='none'))  Multiple comparisons: Bonferroni {r} with(data, pairwise.t.test(Longevity, Treatment, p.adj='bonferroni'))  Multiple comparisons: Tukey's HSD {r} fruitfly.aov <- aov(Longevity ~ Treatment, data=data) options(show.signif.stars=FALSE) summary(fruitfly.aov) TukeyHSD(fruitfly.aov) oldpar <- par(mar=c(4, 10, 4, 2)) # make more room for names plot(TukeyHSD(fruitfly.aov), las=1, col="firebrick") par(oldpar)  Multiple comparisons: Fisher's LSD using mc.R {r} # i j n.i n.j i-j.diff lcl ucl signif? # 1 1pg 1vir 25 25 8.04 -0.2526709 16.332671 no # 2 1pg 8pg 25 25 1.44 -6.8526709 9.732671 no # 3 1pg 8vir 25 25 26.08 17.7873291 34.372671 YES # 4 1pg none 25 25 1.24 -7.0526709 9.532671 no # 5 1vir 8pg 25 25 -6.60 -14.8926709 1.692671 no # 6 1vir 8vir 25 25 18.04 9.7473291 26.332671 YES # 7 1vir none 25 25 -6.80 -15.0926709 1.492671 no # 8 8pg 8vir 25 25 24.64 16.3473291 32.932671 YES # 9 8pg none 25 25 -0.20 -8.4926709 8.092671 no # 10 8vir none 25 25 -24.84 -33.1326709 -16.547329 YES  Multiple comparisons: Bonferroni using mc.R {r} # i j n.i n.j i-j.diff lcl ucl signif? # 1 1pg 1vir 25 25 8.04 -3.938157 20.018157 no # 2 1pg 8pg 25 25 1.44 -10.538157 13.418157 no # 3 1pg 8vir 25 25 26.08 14.101843 38.058157 YES # 4 1pg none 25 25 1.24 -10.738157 13.218157 no # 5 1vir 8pg 25 25 -6.60 -18.578157 5.378157 no # 6 1vir 8vir 25 25 18.04 6.061843 30.018157 YES # 7 1vir none 25 25 -6.80 -18.778157 5.178157 no # 8 8pg 8vir 25 25 24.64 12.661843 36.618157 YES # 9 8pg none 25 25 -0.20 -12.178157 11.778157 no # 10 8vir none 25 25 -24.84 -36.818157 -12.861843 YES  Multiple comparisons: Tukey's HSD using mc.R {r} # i j n.i n.j i-j.diff lcl ucl signif? # 1 1pg 1vir 25 25 8.04 -3.560486 19.640486 no # 2 1pg 8pg 25 25 1.44 -10.160486 13.040486 no # 3 1pg 8vir 25 25 26.08 14.479514 37.680486 YES # 4 1pg none 25 25 1.24 -10.360486 12.840486 no # 5 1vir 8pg 25 25 -6.60 -18.200486 5.000486 no # 6 1vir 8vir 25 25 18.04 6.439514 29.640486 YES # 7 1vir none 25 25 -6.80 -18.400486 4.800486 no # 8 8pg 8vir 25 25 24.64 13.039514 36.240486 YES # 9 8pg none 25 25 -0.20 -11.800486 11.400486 no # 10 8vir none 25 25 -24.84 -36.440486 -13.239514 YES  Compare interval lengths for $8v - none$ {r} Fisher.CI.length <- abs(-33.1326709 + 16.547329) Bonferroni.CI.length <- abs(-36.818157 + 12.861843) Tukey.CI.length <- abs(-36.440486 + 13.239514) CI.lengths <- round(cbind(Fisher=Fisher.CI.length, Bonferroni=Bonferroni.CI.length, Tukey=Tukey.CI.length), 3) rownames(CI.lengths) <- "CI.length" CI.lengths  Group statistics. See "Cannon notes/chap 05/example 5.1-5.6, 5.10-5.11 fruit flies/fruitflies.Rmd" {r} n <- with(data, tapply(Longevity, Treatment, length)) mean <- with(data, round(tapply(Longevity, Treatment, mean), 3)) sd <- with(data, round(tapply(Longevity, Treatment, sd), 3)) idx <- c(5, 1, 3, 2, 4) # idx orders the rows in the table fruitfly.statistics <- cbind(n, mean, sd)[idx, ] fruitfly.statistics  Comparison: $\mu_{8v} - \mu_{none} = 0$ ? See "Cannon notes/chap 05/example 5.1-5.6, 5.10-5.11 fruit flies/fruitflies.Rmd" HT: $H_0 : \mu_{8v} - \mu_{none} = 0$ $H_a : \mu_{8v} - \mu_{none} \not= 0$ Test statistic: $$t = \frac{\bar{y}_{8v} - \bar{y}_{none}}{\sqrt{MSE\Big(\frac{1}{n_{8v}} + \frac{1}{n_{none}})}}$$ {r} y.bar.8v <- fruitfly.statistics[5, 2] y.bar.none <- fruitfly.statistics[1, 2] point.estimate <- y.bar.8v - y.bar.none alpha <- 0.05 df <- 120 mse <- 219.3 # from ANOVA n.8v <- n.none <- 25 se <- sqrt(mse * (1 / n.8v + 1 / n.none)) t <- point.estimate / se t p.value <- 2 * pt(t, df=df) p.value  Contrast: $\frac{1}{2}(\mu_{1v} + \mu_{8v}) - \frac{1}{2}(\mu_{1p} + \mu_{8p}) = 0$ ? See "Cannon notes/chap 05/example 5.1-5.6, 5.10-5.11 fruit flies/fruitflies.Rmd" HT: $H_0 : \frac{1}{2}(\mu_{1v} + \mu_{8v}) - \frac{1}{2}(\mu_{1p} + \mu_{8p}) = 0$ $H_a : \frac{1}{2}(\mu_{1v} + \mu_{8v}) - \frac{1}{2}(\mu_{1p} + \mu_{8p}) \not= 0$ Test statistic: $$t = \frac{\frac{1}{2}\bar{y}_{1v} + \frac{1}{2}\bar{y}_{8v} - \frac{1}{2}\bar{y}_{1p} - \frac{1}{2}\bar{y}_{8p}}{\sqrt{MSE \cdot \frac{1}{4} \cdot \Big(\frac{1}{n_{1v}} + \frac{1}{n_{8v}} + \frac{1}{n_{1p}} + \frac{1}{n_{8p}}\Big)}}$$ {r} y.bar.1v <- fruitfly.statistics[4, 2] y.bar.8v <- fruitfly.statistics[5, 2] y.bar.1p <- fruitfly.statistics[2, 2] y.bar.8p <- fruitfly.statistics[3, 2] point.estimate <- (y.bar.1v + y.bar.8v - y.bar.1p - y.bar.8p) / 2 alpha <- 0.05 df <- 120 mse <- 219.3 # from ANOVA n.1v <- n.8v <- n.1p <- n.8p <- 25 se <- sqrt(mse * (1 / 4) * (1 / n.1v + 1 / n.8v + 1 / n.1p + 1 / n.8p)) t <- point.estimate / se t p.value <- 2 * pt(t, df=df) p.value  ANOVA and regression with indicators: pooled two-sample t-test. {r fig.width=6, fig.height=3.5} group.8v <- data[data$Treatment=="8 virgin", "Longevity"] group.none <- data[data$Treatment=="none", "Longevity"] stripchart(list(group.8v, group.none), pch=20, col="darkred", at=c(1.25, 1.7), method="stack", las=1, ylim=1:2, group.names=c("8v", "none"), xlab="Longevity") t.test(group.8v, group.none, var.equal=TRUE) # pooled t-test  ANOVA and regression with indicators: ANOVA with two groups. {r} two.groups <- with(data, Treatment=="8 virgin" | Treatment=="none") fruitfly.aov <- aov(Longevity ~ Treatment, data=data[two.groups, ]) summary(fruitfly.aov)  ANOVA and regression with indicators: regression with an indicator. {r} fruitfly.lm <- lm(Longevity ~ Treatment, data=data[two.groups, ]) summary(fruitfly.lm) confint(fruitfly.lm)  Treatment effects. {r message=FALSE} library(alr4) plot(allEffects(fruitfly.aov))  One-way ANOVA for means as regression: ANOVA {r} fruitfly.aov5 <-aov(Longevity ~ Treatment, data=data) # 5 levels of Treatment summary(fruitfly.aov5)  One-way ANOVA for means as regression: regression {r} fruitfly.lm5 <-lm(Longevity ~ Treatment, data=data) # 5 levels of Treatment summary(fruitfly.lm5)  Reorder the factors of Treatment to make "none" the reference. {r} data$Treatment <- factor(data$Treatment, levels=c("none", "1 pregnant", "1 virgin", "8 pregnant", "8 virgin")) fruitfly.lm5 <-lm(Longevity ~ Treatment, data=data) # 5 levels of Treatment summary(fruitfly.lm5) aov(fruitfly.lm5)