## Fan Jia, Terrence Jorgensen, and Sunthud Pornprasertmanit
## 16 February 2013
## CRMDA Saturday Seminar: ANOVA and Regression in R
## Answers for exercises

library(car)
library(psych)
library(QuantPsyc)
library(phia)
library(multcomp)
library(rockchalk)

##################################################################
## Exercise 1

# Item 1

library(car)
?Moore
head(Moore)
Moore$partner.status <- factor(Moore$partner.status)
Moore$fcategory <- factor(Moore$fcategory)

# Item 2

table(Moore$partner.status, Moore$fcategory)

# Item 3

boxplot(conformity ~ partner.status, data = Moore)
boxplot(conformity ~ fcategory, data = Moore)

# Item 4

bartlett.test(conformity ~ partner.status + fcategory, data = Moore)
leveneTest(conformity ~ fcategory*partner.status, data=Moore)
leveneTest(conformity ~ fcategory, data=Moore) 
leveneTest(conformity ~ partner.status, data=Moore) 

# Item 5

t.test(conformity ~ partner.status, data=Moore)
m5 <- aov(conformity ~ partner.status, data=Moore)
summary(m5)

# Item 6

m6 <- aov(conformity ~ fcategory, data=Moore)
summary(m6)

# Item 7

TukeyHSD(m6, "fcategory")
library(multcomp)
pairwise <- glht(m6, linfct = mcp(fcategory = "Tukey"))
summary(pairwise, test = adjusted(type = "bonferroni"))

# Item 8

linear <- c(-1, 0, 1)
ctr <- matrix(linear, 1)
linearmodel <- glht(m6, linfct = mcp(fcategory = ctr))
summary(linearmodel)

# Item 9

quadratic <- c(1, -2, 1)
ctr2 <- rbind(linear, quadratic)
polynomial <- glht(m6, linfct = mcp(fcategory = ctr2))
summary(polynomial, test = adjusted(type = "bonferroni"))

# Item 10

m10 <- aov(conformity ~ fcategory*partner.status, data=Moore)
summary(m10)
Anova(m10, type=2)

# Item 11

library(phia)
interactionMeans(m10)
plot(interactionMeans(m10))

# Item 12

testInteractions(m10, fixed = "partner.status", across = "fcategory", adjustment = "bonferroni")
testInteractions(m10, fixed = "fcategory", across = "partner.status", adjustment = "bonferroni")

# Item 13

testInteractions(m10, fixed = "partner.status", pairwise = "fcategory", adjustment = "bonferroni")


##################################################################
## Exercise 2

# Item 1a
library(car)
help(package="car")

data(Highway1)
# ?Highway1
head(Highway1)
mod1 <- lm(rate ~ acpt, data=Highway1)
summary(mod1)

# Item 1b
mod2 <- lm(rate ~ acpt + sigs1, data=Highway1)
summary(mod2)

# Item 1c
anova(mod1, mod2, test="F")

# Item 1d
mod3 <- lm(rate ~ acpt * sigs1, data=Highway1)
summary(mod3)
anova(mod2, mod3, test="F")

library(rockchalk)
plotSlopes(mod3, plotx = "acpt", modx = "sigs1")

# Item 2a
data(Salaries)
#?Salaries
head(Salaries)
levels(Salaries$discipline)
contrasts(Salaries$discipline)

# Item 2b
Salaries$discipline.rev <- relevel(Salaries$discipline, "B")
contrasts(Salaries$discipline.rev)
mod4 <- lm(salary ~ discipline.rev, data=Salaries)
summary(mod4)

# Item 2c
mod5 <- lm(salary ~ yrs.service * discipline, data=Salaries)
summary(mod5)

plotSlopes(mod5, plotx = "yrs.service", modx = "discipline")

# Item 2d
qqnorm(mod5$residuals)
qqline(mod5$residuals)
library(QuantPsyc)
plotNormX(mod5$residuals)

plot(mod5$residuals ~ mod5$fitted.values)
plot(mod5$residuals ~ mod5$model$yrs.service)