##############################
#
# Script: Solutions 4 (assignment04a.R)
#
##############################



####################
# Problem 04.1

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

png("boxplot041.png",width=600,height=800)
par(mar=c(5,4,2,2)+0.1)
boxplot(height~fertilizer, data=wheat, las=1, xlab="Fertilizer Type", ylab="Plant Height (in)")
dev.off()

# Normality tests
with(wheat, shapiro.test(height[fertilizer=="N18P51K20"]) )   # p = 0.57
with(wheat, shapiro.test(height[fertilizer=="N13P00K44"]) )   # p = 0.42
with(wheat, shapiro.test(height[fertilizer=="N18P18K18"]) )   # p = 0.70

# Equal variance tests
bartlett.test(height~fertilizer, data=wheat)    # p = 0.0242
fligner.test(height~fertilizer, data=wheat)     # p = 0.1579

# So, which results to use?

# Do both means tests and see if there is an appreciable difference

model041a <- aov(height~fertilizer, data=wheat)
summary(model041a)       # p << 0.0001

model041b <- kruskal.test(height~fertilizer, data=wheat)
print(model041b)         # p << 0.0001

# As the conclusions are the same, we can feel more comfortable about our conclusions





####################
# Problem 04.2

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

png("boxplot042.png",width=600,height=800)
par(mar=c(5,4,2,2)+0.1)
boxplot(mfri~biome, data=bio, las=1, xlab="Biome", ylab="Mean Fire Return Interval (yr)")
dev.off()


# Normality tests
with(bio, shapiro.test(mfri[biome=="STR"]) )   # p = 0.18
with(bio, shapiro.test(mfri[biome=="TBF"]) )   # p = 0.63
with(bio, shapiro.test(mfri[biome=="TCF"]) )   # p = 0.73
with(bio, shapiro.test(mfri[biome=="TS"]) )    # p = 0.55
with(bio, shapiro.test(mfri[biome=="XER"]) )   # p = 0.85

# Equal variance tests
bartlett.test(mfri~biome, data=bio)            # p << 0.0001
fligner.test(mfri~biome, data=bio)             # p = 0.002

model042 <- kruskal.test(mfri~biome, data=bio)
print(model042)                                # p << 0.0001