#
#Add storm as a factor to the data frame "tacoma"
#
storm<-tacoma[,4]
storm<-factor(storm)
tacoma.frame<-data.frame(tacoma.frame,storm=storm)
#
#Check for approximate normality of residuals
#
hist(tacoma.aov$residuals)
qqnorm(tacoma.aov$residuals)
#
#Did not look too good--try a transformation of the dependent variable
#
tacoma.rt3<-aov(so4^(1/3)~year*station,data=tacoma)
hist(tacoma.rt3$residuals)
qqnorm(tacoma.rt3$residuals)
#
#Looked better. Now look at the residuals for each value of the factors
#
plot(year, tacoma.rt3$residuals)
plot(station, tacoma.rt3$residuals)
#
#Look at the residuals for each of the 10 storms
#
plot(storm,tacoma.rt3$residuals)
#
#Look at residuals against observed and against fitted values
#
plot(so4^(1/3),tacoma.rt3$residuals)
plot(tacoma.rt3$residuals,tacoma.rt3$fitted.values)
#
#We can do interaction plots as well
#
interaction.plot(year,station,so4^(1/3))
interaction.plot(station,year,so4^(1/3))
#
#See what the consequence is of using an additive model (without 
#interaction between year and station)
#
plot(tacoma.rt3$residuals,aov(so4^(1/3)~year+station)$residuals)