## April 5, 2007 550.413 Section


# Prepping the data
popcorn = read.table("c:/413/popcorn.txt", header=T)
attach(popcorn)
names(popcorn)
#Make sure all factors are treated as categorical variables.
PopcornBrand=as.factor(PopcornBrand)
Power = as.factor(Power)
OvenBrand = as.factor(OvenBrand)
PoppingTime = as.factor(PoppingTime)
# Create new variable which indicates which treatment an observation is in.
Treatment = as.factor(paste(PopcornBrand, Power, OvenBrand, PoppingTime, sep=""))
popcorn=data.frame(popcorn, Treatment)
popcorn
# Note: We won't use Treatment this week.  We will use it when we look at testing again.

# Problem 1
# part a
PopModel = aov(PctEdible~PopcornBrand)
tapply(PctEdible, PopcornBrand, mean) # find the mean of PctEdible for each brand.
PopFit=PopModel$fitted.values # save fitted values of model as PopFit
PopResid = round(PopModel$residuals, 2)  # save residuals of model as PopResid
data.frame(PctEdible, PopFit, PopResid) 
# put response variable, fitted values and residuals in a data frame and display

#part b
TukeyHSD(PopModel, "PopcornBrand")

# part c
plot(PopcornBrand, PopResid)
plot(as.numeric(PopcornBrand), PopResid)

#part d
newEdible = (PctEdible)^0.5  #create transformed response variable
PopModel2 = aov(newEdible~PopcornBrand)
PopFit2 = PopModel2$fitted.values
PopResid2 = round(PopModel2$residuals, 2)
data.frame(newEdible, PopFit2, PopResid2)

#part e
plot(as.numeric(PopcornBrand), PopResid2) #residual plot
qqnorm(PopResid2) # normal probability plot

#Problem2
#part b
# two ways to fit this model:
PopModel3 = aov(PctEdible ~ PopcornBrand*OvenBrand)
anova(PopModel3)
PopModel4 = aov(PctEdible ~ PopcornBrand + OvenBrand + PopcornBrand:OvenBrand)
anova(PopModel4)

PopFit3 = PopModel3$fitted.values
PopResid3 = round(PopModel3$residuals, 2)
data.frame(PctEdible, PopFit3, PopResid3)

#part c
plot(PopFit3, PopResid3)

#part d
qqnorm(PopResid3)
hist(PopResid3)

#part e
interaction.plot(PopcornBrand, OvenBrand, PctEdible)

detach(popcorn)