vgrabovets
12/12/2016 - 3:54 PM

ANOVA

ANOVA

# one way ANOVA ####
plant.df = PlantGrowth
plant.df$group = factor(plant.df$group, labels = c("Control", "Treatment 1", "Treatment 2"))

require(ggplot2)

ggplot(plant.df, aes(x = group, y = weight)) +
     geom_boxplot(fill = "grey80", colour = "blue") +
     scale_x_discrete() + 
     xlab("Treatment Group") +
     ylab("Dried weight of plants")

plant.mod1 = lm(weight ~ group, data = plant.df)
plant.mod2 = aov(weight ~ group, data = plant.df)
summary(plant.mod1)
anova(plant.mod1)
TukeyHSD(plant.mod2, which = "group")
confint(plant.mod1)

plant.mod = data.frame(Fitted = fitted(plant.mod1), Residuals = resid(plant.mod1), Treatment = plant.df$group)

ggplot(plant.mod, aes(Fitted, Residuals, colour = Treatment)) + geom_point()
#__________________________________________________________________________________________________________________________________

# two way ANOVA ####
delivery.df = data.frame(
     Service = c(rep("Carrier 1", 15), rep("Carrier 2", 15),
                 rep("Carrier 3", 15)),
     Destination = c(rep(c("Office 1", "Office 2", "Office 3",
                           "Office 4", "Office 5"), 9)),
     Time = c(15.23, 14.32, 14.77, 15.12, 14.05,
              15.48, 14.13, 14.46, 15.62, 14.23, 15.19, 14.67, 14.48, 15.34, 14.22,
              16.66, 16.27, 16.35, 16.93, 15.05, 16.98, 16.43, 15.95, 16.73, 15.62,
              16.53, 16.26, 15.69, 16.97, 15.37, 17.12, 16.65, 15.73, 17.77, 15.52,
              16.15, 16.86, 15.18, 17.96, 15.26, 16.36, 16.44, 14.82, 17.62, 15.04)
)

ggplot(delivery.df, aes(Time, Destination, colour = Service)) + geom_point()

delivery.mod1 = aov(Time ~ Destination*Service, data = delivery.df)
summary(delivery.mod1)

delivery.res = delivery.df
delivery.res$M1.Fit = fitted(delivery.mod1)
delivery.res$M1.Resid = resid(delivery.mod1)

ggplot(delivery.res, aes(M1.Fit, M1.Resid, colour = Service)) + geom_point() +
     xlab("Fitted Values") + ylab("Residuals")

ggplot(delivery.res, aes(M1.Fit, M1.Resid, colour = Service)) +
     geom_point() + xlab("Fitted Values") + ylab("Residuals") +
     facet_wrap( ~ Destination)

ggplot(delivery.res, aes(M1.Fit, M1.Resid, colour = Destination)) +
     geom_point() + xlab("Fitted Values") + ylab("Residuals") +
     facet_wrap( ~ Service)

ggplot(delivery.res, aes(sample = M1.Resid)) + stat_qq()
TukeyHSD(delivery.mod1, which = "Service")

delivery.hsd = data.frame(TukeyHSD(delivery.mod1, which = "Service")$Service)
delivery.hsd$Comparison = row.names(delivery.hsd)

ggplot(delivery.hsd, aes(Comparison, y = diff, ymin = lwr, ymax = upr)) +
     geom_pointrange() + ylab("Difference in Mean Delivery Time by Service") +
     coord_flip()