# Orange Example -- source code file orange.r

# Load and print orange_orig data frame.
orange_orig = read.table("orange.txt", header=T, skip=2)
cat("orange_orig data frame\n")
print(orange_orig)

# Send plots to a single pdf file.
pdf("orange.pdf")

# Create plots of yield vs. factor1, using factor2
# as plotting symbol.  Then reverse factors.
# Connect the means of equal plotting symbols with
# line segments to produce interaction plots.
attach(orange_orig)
plot(pesticide, yield, pch=as.character(variety),
     main="Orange Data Frame, plotting symbol=variety",
     xlab="Pesticide",
     ylab="Yield")
plot(variety, yield, pch=as.character(pesticide),
     main="Orange Data Frame, plotting symbol=pesticide",
     xlab="Variety",
     ylab="Yield")
detach(orange_orig)

# Create orange data frame where variety and 
# pesticide are discrete factors.
orange = data.frame(variety=as.factor(orange_orig$variety),
                     pesticide=as.factor(orange_orig$pesticide),
                     yield=orange_orig$yield)
cat("Verify that variety and pesticide are discrete factors:\n")
print(orange$variety)
print(orange$pesticide)

# Attach data frame for plots.
attach(orange)

# Compute additive model.
additive = lm(yield ~ variety + pesticide)
print(summary(additive))
p = fitted(additive)
r = residuals(additive)


# Create residual plots.
plot(p, r, main="Residual Plot for Orange Data Frame",
	         xlab="Predicted Values",
           ylab="Residuals")
abline(h=0, lty="dashed")
plot(p, r, pch=as.character(variety),
	         main="Residual Plot for Orange Data Frame",
	         sub="Use variety as plotting symbol",
	         xlab="Predicted Values",
           ylab="Residuals")
abline(h=0, lty="dashed")

dev.off( )

# Compute model with interaction.
interact = lm(yield ~ variety + pesticide + variety*pesticide)
print(summary(interact))

summary(aov(yield ~ variety + pesticide + variety*pesticide))

# Show ANOVA sum of squares tables for both models
cat("ANOVA table for additive model:\n")
print(aov(yield ~ variety + pesticide))
cat("ANOVA table for model with interaction:\n")
print(aov(yield ~ variety + pesticide + variety*pesticide))