# R program for GO Problem 13.8, VHM 802
rockets <- read.csv("r:/ch13pr8.csv")
head(rockets)
# recode as factors in preparation for modelling
rockets$days <- as.factor(rockets$days)
rockets$temp <- as.factor(rockets$temp)
rockets$slant <- as.factor(rockets$slant)

# ANOVA for model with volley error term
summary (aov( azimuth~ temp*slant+days+Error(temp:slant:days), data=rockets))
# diagnostics
rockets.within <- lm( azimuth ~ temp*slant+days+temp:slant:days, data=rockets)
summary(rockets.within)
within.diag <- rockets
within.diag$fit <- fitted(rockets.within)
within.diag$stdres <- rstandard(rockets.within)
plot(stdres~fit, data=within.diag)
qqnorm(within.diag$stdres)
qqline(within.diag$stdres)
shapiro.test(within.diag$stdres)

# analysis of mean azimuth error in a group
library(doBy)
rockets.mean <- summaryBy( azimuth~ temp+slant+days, data=rockets, FUN=mean)
rockets.mean.lm <- lm( azimuth.mean ~ temp*slant+days, data=rockets.mean)
summary(rockets.mean.lm)
anova(rockets.mean.lm)
# diagnostics and further analysis by usual methods

# analysis of mean azimuth error without one extreme observation
r <- rockets.mean # to facilitate notation
rockets.mean.excl1 <- lm( azimuth.mean ~ temp*slant+days, 
   data=r[!((r$days==2)&(r$slant==2)&(r$temp==3)),])
summary(rockets.mean.excl1)
anova(rockets.mean.excl1) 
# note: ANOVA table based on sequential instead of partial sum of squares

# analysis of mean azimuth error without two extreme observations
rockets.mean.excl2 <- lm( azimuth.mean ~ temp*slant+days, 
   data=r[!((r$days==2)&(r$slant==2)&(r$temp==3))|((r$days==2)&(r$slant==1)&(r$temp==3)),])
summary(rockets.mean.excl2)
anova(rockets.mean.excl2)
# note: ANOVA table based on sequential instead of partial sum of squares

# analysis of median azimuth error in a group
rockets.median <- summaryBy( azimuth~ temp+slant+days, data=rockets, FUN=median)
rockets.median.lm <- lm( azimuth.median ~ temp*slant+days, data=rockets.median)
summary(rockets.median.lm)
anova(rockets.median.lm)

# random effects model for the full data, ignoring the bad-looking residuals
library(nlme)
rockets$volley <- 100*as.numeric(rockets$temp)+10*as.numeric(rockets$slant)+as.numeric(rockets$days) # unique id for volleys
head(rockets) # to check calculation
rockets.mixed <- lme( azimuth ~ temp*slant+days, random = ~1|volley, data=rockets)
summary(rockets.mixed)
anova(rockets.mixed)

# full data model allowing for negative correlations
rockets.gls <- gls( azimuth ~ temp*slant+days, correlation = corCompSymm(form=~1|volley), data=rockets)
summary(rockets.gls)
anova(rockets.gls)