######################################################################
# Wheaties data #
# #
######################################################################
######################################################################
#Read in data
wheaties <- read.csv(file = "wheaties.csv")
wheaties
######################################################################
# Plots
#boxplot(formula = Response ~ Design, data = wheaties, main = "Box and dot plot", ylab = "Sales",
# xlab = "Design", pars = list(outpch=NA)) #Not helpful here due to small sample size for each treatment level
stripchart(x = wheaties$Response ~ wheaties$Design, lwd = 2, col = "red",
method = "jitter", vertical = TRUE, pch = 1, main =
"Dot plot", ylab = "Number sold", xlab = "Design type")
aggregate(formula = Response ~ Design, data = wheaties, FUN = mean)
aggregate(formula = Response ~ Design, data = wheaties, FUN = sd)
trt.means <- aggregate(formula = Response ~ Design, data = wheaties, FUN = mean)
trt.means
set.seed(1819) #Keep same jittering
stripchart(x = wheaties$Response ~ wheaties$Design, lwd = 2, col = "red",
method = "jitter", vertical = TRUE, pch = 1, main =
"Dot plot for SSE and SST", ylab = "Number sold", xlab = "Design type")
move.amount <- 0.25
segments(x0 = trt.means[,1] - move.amount, y0 = trt.means[,2],
x1 = trt.means[,1] + move.amount, y1 = trt.means[,2], lwd = 2, lty = "solid", col = "blue")
abline(h = mean(wheaties$Response), lwd = 2, lty = "dashed", col = "darkgreen")
legend(x = 1, y = 30, legend = c("Design means", "Overall mean"), col = c("blue", "darkgreen"),
lwd = c(2,2), lty = c("solid", "dashed"), bty = "n")
######################################################################
# Analysis
t <- 4
N <- 10
qf(p = 0.95, df1 = t - 1, N - t)
pf(q = 86/7.67, df1 = t - 1, df2 = N - t)
#ANOVA model
mod.fit <- aov(formula = Response ~ factor(Design), data = wheaties)
mod.fit
summary(object = mod.fit)
names(x = mod.fit)
#Another way - this will make more sense after our chapter on regression
mod.fit.lm <- lm(formula = Response ~ factor(Design), data = wheaties)
names(x = mod.fit.lm)
mod.fit.lm$coefficients
summary(object = mod.fit.lm)
anova(object = mod.fit.lm)
#Multiple comparisons
pairwise.t.test(x = wheaties$Response, g = wheaties$Design, p.adjust.method = "none",
alternative = "two.sided")
pairwise.t.test(x = wheaties$Response, g = wheaties$Design, p.adjust.method = "bonferroni",
alternative = "two.sided")
TukeyHSD(x = mod.fit, conf.level = 0.95)
temp <- TukeyHSD(x = mod.fit, conf.level = 0.95)
names(temp)
#Compare to Bonferroni p-values
LSD <- pairwise.t.test(x = wheaties$Response, g = wheaties$Design, p.adjust.method = "none",
alternative = "two.sided")
names(LSD)
LSD$p.value
LSD$p.value*6
################################
#Other ways to do the multiple comparisons in R
# Could not get to work in R 4.0.2
library(agricolae) #http://cran.r-project.org/web/packages/agricolae/index.html
#LSD - save results into an object and then print. Oddly, nothing is given
# in the console window when LSD.test() is just run without first saving results
# into an object. This is due to the last line of the function being "invisible(output)".
save.LSD <- LSD.test(y = mod.fit, trt = "factor(Design)", alpha = 0.05,
group = FALSE, p.adj = "none", main = "Wheaties box design")
names(save.LSD)
save.LSD
# Verification of LCL value in $means first row: 15-2.446912*sqrt(7.6666667)/sqrt(2)
#Bonferroni
save.Bon <- LSD.test(y = mod.fit, trt = "factor(Design)", alpha = 0.05, group = FALSE, p.adj = "bonferroni", main = "Wheaties box design")
save.Bon
save.HSD <- HSD.test(y = mod.fit, trt = "factor(Design)", alpha = 0.05,
group = FALSE, main = "Wheaties box design")
save.HSD
library(package = multcomp)
#HSD
HSD <- glht(model = mod.fit, linfct = mcp("factor(Design)" = "Tukey"))
summary(object = HSD)
confint(object = HSD, level = 0.95)
plot(HSD)
#The numeric levels of Design make things harder to work with in the multcomp package for some multiple comparisons
# Below is how I changed these levels
wheaties$Design2 <- factor(x = wheaties$Design, levels = 1:4, labels = c("D1", "D2", "D3", "D4"))
wheaties
#ANOVA model
mod.fit2 <- aov(formula = Response ~ Design2, data = wheaties)
summary(object = mod.fit2)
#List out the comparisons - be careful about syntax
comparisons <- rbind("D1 - D2 = 0", "D1 - D3 = 0", "D1 - D4 = 0",
"D2 - D3 = 0", "D2 - D4 = 0", "D3 - D4 = 0")
comparisons
compare.means <- glht(model = mod.fit2, linfct = mcp(Design2 = comparisons))
compare.means
summary(object = compare.means, test = adjusted(type = "none"))
summary(object = compare.means, test = adjusted(type = "bonferroni"))
plot(compare.means)
confint(object = compare.means, level = 0.95, calpha = qt(p = 0.975, df = mod.fit$df.residual))
#This is another way to do the multiple comparisons using the original model fit
#Construct contrasts
mod.fit$coefficients
K <- rbind("1 - 2" = c(0, 1, 0, 0),
"1 - 3" = c(0, 0, 1, 0),
"1 - 4" = c(0, 0, 0, 1),
"2 - 3" = c(0, 1, -1, 1),
"2 - 4" = c(0, 1, 0, -1),
"3 - 4" = c(0, 0, 1, -1))
K
#A little different way
#matrix(data = c(0, 1, 0, 0,
# 0, 0, 1, 0,
# 0, 0, 0, 1,
# 0, 1, -1, 0,
# 0, 1, 0, -1,
# 0, 0, 1, -1),
# nrow = 6, ncol = 4, byrow = TRUE)
compare.means <- glht(model = mod.fit, linfct = K)
compare.means
summary(object = compare.means, test = adjusted(type = "none"))
summary(object = compare.means, test = adjusted(type = "bonferroni"))
plot(compare.means)
confint(object = compare.means, level = 0.95)
#HSD
HSD <- glht(model = mod.fit, linfct = mcp("factor(Design)" = "Tukey"))
summary(object = HSD)
confint(object = HSD, level = 0.95)
plot(HSD)
#