# R-program for Supplementary Exercise 13.31 of IPS7e
adv <- read.csv("r:/Chapter 13/ex13_031.csv")

library(doBy) # package to facilitate processing by groups
adv.means <- summaryBy( eprice ~ promo+percent, data=adv, FUN=c(length,mean,sd))
adv.means

# ANOVA
adv.full <- lm(eprice ~ as.factor(promo)*as.factor(percent), data=adv)
summary(adv.full)
anova(adv.full)

# interaction plots
library(lattice)
xyplot(eprice.mean ~ percent, groups=promo, type="b", data=adv.means)
xyplot(eprice.mean ~ promo, groups=percent, type="b", data=adv.means)

# residual analysis
adv.diag <- adv # all diagnostics to be stored in adv.diag (for convenience, not at all necessary)
adv.diag$fit <- fitted(adv.full) 
adv.diag$stdres <- rstandard(adv.full)
plot(stdres~fit, data=adv.diag)
hist(adv.diag$stdres, nclass=12) # no overlaid normal distribution
qqnorm(adv.diag$stdres)
qqline(adv.diag$stdres)
shapiro.test(adv.diag$stdres)

# calculation of tstar values from t(144)
qt(0.025, 144, lower.tail=F)
qt(0.025/6, 144, lower.tail=F) # for Bonferroni method

# general method for pairwise comparisons: manual calculation
# first refit model without interaction
adv.red <- <- lm(eprice ~ as.factor(promo)+as.factor(percent), data=adv)
summary(adv.red)
anova(adv.red)
coef <- coef(adv.red)
vcov <- vcov(adv.red)
# comparisons for promo
wpromo<- matrix(c(0,1,0,0,0,0,0, 0,0,1,0,0,0,0, 0,0,0,1,0,0,0,
                  0,1,-1,0,0,0,0, 0,1,0,-1,0,0,0, 0,0,1,-1,0,0,0), ncol=6); wpromo
est<- t(wpromo)%*%coef
var <- diag(t(wpromo)%*%vcov%*%wpromo)
se <- matrix(sqrt(var), ncol=1)
t<- est/se
p<- 2*(1-pt(abs(t),df.residual(adv.red)))
# multiple testing
t(p) # unadjusted P-values
p.adjust(p, method="bonferroni")
# comparisons for perc
wperc<- matrix(c(0,0,0,0,1,0,0, 0,0,0,0,0,1,0, 0,0,0,0,0,0,1,
                  0,0,0,0,1,-1,0, 0,0,0,0,1,0,-1, 0,0,0,0,0,1,-1), ncol=6); wperc
est<- t(wperc)%*%coef
var <- diag(t(wperc)%*%vcov%*%wperc)
se <- matrix(sqrt(var), ncol=1)
t<- est/se
p<- 2*(1-pt(abs(t),df.residual(adv.red)))
# multiple testing
t(p) # unadjusted P-values
p.adjust(p, method="bonferroni")