##############################
#
# Script: Solutions 3 (assignment03a.R)
#
##############################



####################
# Problem 03.1

png("boxplot031.png",height=800,width=600)
par(mar=c(5,4,2,2)+0.1)
boxplot(Petal.Width~Species, data=iris, las=1, ylab="Petal Width (cm)", xlab="Species", ylim=c(0,3))
dev.off()

# Test of equal variances
bartlett.test(Petal.Width~Species, data=iris)   # p << 0.0001

# Should not use AOV, but will anyway
model1 <- aov(Petal.Width~Species, data=iris)
summary(model1)                                 # Reject equal means (p << 0.0001)

# Should use Kruskal-Wallis
model1b <- kruskal.test(Petal.Width~Species, data=iris)
print(model1b)                                  # Reject equal means (p << 0.0001)


# Why did I not test the Normality assumption of AOV?
# If the data failed the assumptions of AOV, why did I tell you to do it anyway??




####################
# Problem 03.2

d2 <- read.csv("http://courses.kvasaheim.com/stat40x3/data/gdpcap.csv", header=TRUE)

png("boxplot032.png",height=800,width=600)
par(mar=c(5,7,2,2)+0.1)
boxplot(hig~region, data=d2, las=1, horizontal=TRUE, xlab="Honesty in Government index")
dev.off()

# Equality of variances test
bartlett.test(hig~region, data=d2)              # p << 0.0001

model2 <- kruskal.test(hig~region, data=d2)
print(model2)                                   # Reject equal means (p << 0.0001)

# Why did I not test the Normality assumption of AOV?




####################
# Problem 03.3

png("boxplot033.png",height=800,width=600)
par(mar=c(5,6,2,2)+0.1)
boxplot(weight~feed, data=chickwts, horizontal=TRUE, las=1, xlab="Chick weight (g)")
dev.off()

# Equal variances test
bartlett.test(weight~feed, data=chickwts)                   # p=0.66
fligner.test(weight~feed, data=chickwts)                    # p=0.58

# Normality tests
shapiro.test(chickwts$weight[chickwts$feed=="sunflower"])   # p=0.36
shapiro.test(chickwts$weight[chickwts$feed=="soybean"])     # p=0.51
shapiro.test(chickwts$weight[chickwts$feed=="meatmeal"])    # p=0.96
shapiro.test(chickwts$weight[chickwts$feed=="linseed"])     # p=0.90
shapiro.test(chickwts$weight[chickwts$feed=="horsebean"])   # p=0.53
shapiro.test(chickwts$weight[chickwts$feed=="casein"])      # p=0.26


# Assumptions not violated, so use AOV
model3 <- aov(weight~feed, data=chickwts)
summary(model3)                                 # Reject equal means (p << 0.0001)

# Why do we 'want' to use AOV instead of Kruskal-Wallis? (Two reasons.)