##############################
#
#  assignment06a.R
#
##############################




####################
# Problem 06.2

tm <- 1:9
ex <- c(51,42,34,28,23,19,15,13,10)
ob <- c(49,38,30,23,18,14,11,8,7)
X2 <- sum((ob-ex)^2/ex)
p  <- 1-pchisq(X2,df=length(ex)-1)
p


# Alternatively

source("http://courses.kvasaheim.com/stat40x3/scripts/chisq.1d.test.R")
chisq.1d.test(ex,ob)

# So, we conclude that there is no statistically significant 
# effect due to the caffeine.




# Graph

png("062plot.png", width=600, height=600)
par(mar=c(5,4,2,2)+0.1)
plot(ex~tm, pch=16, col=2, las=1, ylim=c(0,65), ylab="Caffeine Concentration (total)",
     xlab="Time Elapsed (h)", xlim=c(0,10), type="b")
points(ob~tm, pch=16, col=4, type="b")
legend("topright", legend=c("Without Caffeine", "With Caffeine"), col=c(2,4), 
     pch=16, lty=1, bty="n")
abline(h=62.5/2)
text(10,62.5/2+1.1, "Half-Concentration Line", adj=c(1,0.5) )
dev.off()



####################
# Problem 06.3

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


##########
# One way to tidy up the data


rejex <- sri$REJECTED
total <- sri$TOTAL
rajwn <- sri$WONPROVINCE
prvnc <- sri$PROVINCE

drop <- which(total<6)
rejex <- rejex[-drop]
total <- total[-drop]
rajwn <- rajwn[-drop]
prvnc <- prvnc[-drop]

drop <- which(rejex<6)
rejex <- rejex[-drop]
total <- total[-drop]
rajwn <- rajwn[-drop]
prvnc <- prvnc[-drop]


pinvalid <- rejex/total


##########
# Graph

w <- "antique white"
l <- "tomato"

boxplot(pinvalid~rajwn, ylim=c(0,0.05),
   las=1, ylab="Proportion of Votes Invalidated",
   xaxt="n", xlab="Rajapaksa Performance", col=c(l,w) )
axis(1, label=c("Lost","Won"), at=c(1,2))

# OR

par(mar=c(5,8,2,2)+0.1)
boxplot(pinvalid~rajwn, horizontal=TRUE, las=1, ylim=c(0,0.05),
     xlab="Proportion of vote declared Invalid", yaxt="n", col=c(l,w) )
axis(2, label=c("Rajapaksa Lost","Rajapaksa Won"), at=c(1,2), las=1)

# OR

png("063boxplota.png", width=600, height=160)
par(mar=c(0,12,2,2)+0.1)
boxplot(pinvalid~rajwn, horizontal=TRUE, las=1, ylim=c(0,0.05),
     xaxt="n", names=c("Rajapaksa Lost","Rajapaksa Won"), col=c(l,w),
	at=c(2,1) )
dev.off()


##########
# Statistics

mean(pinvalid[rajwn==0])
mean(pinvalid[rajwn==1])

median(pinvalid[rajwn==0])
median(pinvalid[rajwn==1])

var(pinvalid[rajwn==0])
var(pinvalid[rajwn==1])




shapiro.test(pinvalid[rajwn==0])  # Not Normal (barely)
shapiro.test(pinvalid[rajwn==1])  # Not Normal


var.test(pinvalid~rajwn)          # Not equal variance
bartlett.test(pinvalid,rajwn)     # Not equal variance

wilcox.test(pinvalid~rajwn)       # Not equal medians (or means)



# If you sought a transformation, you probably did not find one. 
# Besides, we reject so strongly that there is little hope that
# the analysis of variance will give a different conclusion.

# Therefore, we conclude that the proportion of the vote declared 
# invalid (REJECTED/TOTAL) is much higher in provinces supporting 
# Fonseka than in provinces supporting Rajapaksa.

# Also:

png("063boxplotb.png", width=600, height=640)
par(mar=c(5,12,0,2)+0.1)
w <- "antique white"
l <- "tomato"

boxplot(pinvalid~prvnc, horizontal=TRUE, las=1, col=c(w,l,w,w,l,w,w,w,w),
	xlab="Proportion of vote invalidated", ylim=c(0,0.05),
	at=c(6,8,4,3,9,2,1,7,5))

dev.off()

kruskal.test(pinvalid~prvnc)
library(agricolae)
kruskal(pinvalid,prvnc)