# -----------------------------------------------
cat("\n# -------------------  butterfly -------------------------------\n")
butterfly = read.table("http://www-stat.stanford.edu/~rag/stat141/exs/monarch.txt",header=T)
attach(butterfly)

writeLines("\n----  Summary ----\n")
print(summary(butterfly))
writeLines("\n")
print(str(butterfly))

boxplot(butterfly[,2:7])

X11()
boxplot(wingspan, main="wingspan")

b = butterfly
b$sx[b$sex == "Female"] = -1 
b$sx[b$sex == "Male"] = 1 
b = b[,2:8]
X11()
plot(b)

X11()
boxplot(wetwt~sex, main="wetwt")

cat("\n# -- method 1 --:\n")
female = butterfly[sex=="Female",]   # whats the comma for?
print(t.test(female$wetwt,conf.level=.99))
cat("\n# -- method 2 --:\n")
female = wetwt[sex=="Female"]
print(t.test(female,conf.level=.99))

female = wetwt[sex=="Female" & wingspan < 10]
print(t.test(female,conf.level=.99))

male = wetwt[sex=="Male"]
print(t.test(male,conf.level=.99))

male = dryabd[sex=="Male"]
female = dryabd[sex=="Female"]
t.test(male,female,alternative="less")

X11()
boxplot(dryabd~sex)

# -----------------------------------------------
cat("\n# --------------  heights ----------------------------------\n")

CIsimulate = function(true.mean=70, true.sd=5, sample.size=30, 
	num.samples=500, conf.level=.95){
	
	# generate random variables
	tmp = rnorm(sample.size*num.samples, mean=true.mean,sd = true.sd)
	# make a matrix from the random variables
	samples = matrix(tmp, ncol=sample.size)  
	
	alpha = 1-conf.level
	crit.val = qnorm(c(alpha/2, 1-alpha/2))
	
	# sample mean +/- t*SE
	CI.lower = apply(samples, 1, mean) + crit.val[1]*apply(samples,1,sd)/sqrt(sample.size)
	CI.upper = apply(samples, 1, mean) + crit.val[2]*apply(samples,1,sd)/sqrt(sample.size)
	
	CIs = cbind(CI.lower, CI.upper)
	captured = (true.mean > CIs[,1])*(true.mean < CIs[,2])
	num.correct = sum(captured)
	prop.correct = num.correct/num.samples
	cat("\nDesired confidence level: ", conf.level, sep="")
	cat("\nTotal number of simulations: ", num.samples, sep="")
	cat("\nNumber of intervals capturing the true mean: ", num.correct, sep="")
	cat("\nProportion of intervals that correctly captured: ", prop.correct, "\n", sep="")
	}
	
# -----------------------------------------------
cat("\n# ---------------------------------- centipedes --------------------\n")
lengths = c(14.5, 17.2, 14.2, 14.6, 15.4, 13.7, 12.3)
print(t.test(lengths, mu=16))

# -----------------------------------------------
cat("\n# ----------------------- knee ------------------------\n")
knee = read.table("http://www-stat.stanford.edu/~rag/stat141/ass/knee.dat", header=TRUE)
print (summary(knee))
print(str(knee))
X11()
boxplot(knee$days~knee$fitness, main="Knee: fitness vs days", xlab = "fitness", ylab="days")

cat("\nmean:\n")
print(tapply(knee$days, knee$fitness, mean))
cat("\nvar:\n")
print(tapply(knee$days, knee$fitness, var))

writeLines("\n\n --- Incorrect ---\n")
AOV = aov(knee$days ~ knee$fitness)
writeLines("\n\n --- AOV ----\n")
print(AOV)
writeLines("\n --- Summary AOV ----\n")
print (summary(AOV))

writeLines("\n\n --- Correct ---\n")
knee$fitness = as.factor(knee$fitness)   # fitness must be a factor, not numeric
print(str(knee))
AOV = aov(knee$days ~ knee$fitness)
writeLines("\n\n --- AOV ----\n")
print(AOV)
writeLines("\n --- Summary AOV ----\n")
print (summary(AOV))

writeLines("\n\n --- TukeyHSD ---\n")
print(TukeyHSD(AOV))

# -----------------------------------------------
cat("\n# ------- butterfly ANOVA ----------\n")
detach(butterfly)
butterfly = read.table("http://www-stat.stanford.edu/~rag/stat141/exs/monarch.txt",header=T)
butterfly = butterfly[1:10,]
attach(butterfly)


X11()
#jpeg(file="anova.jpeg",height=900, width=800)
par(mfrow=c(2,2))

# --- A
boxplot(wetwt,drywt, names=c("wetwt","drywt"), main= "A")

# --- B
both = c(wetwt,drywt)
grandmean = mean(c(drywt,wetwt))
plot(both, xlab="", ylab="weight", main= "B")
abline(grandmean,0)
for (i in seq(1,20)) lines( c(i,i), c(grandmean, both[i]))
text(15,grandmean+10, "Grand Mean")

# --- C
plot(wetwt,xlim = c(1,21),ylim=c(200,650),xlab="",ylab="weight", main= "C")
abline(mean(wetwt),0)
for (i in 1:length(wetwt)) lines(c(i,i),c(mean(wetwt),wetwt[i]))
points(drywt~seq(11,20))
abline(mean(drywt),0)
for (i in 1:length(drywt)) lines(c(i+10,i+10),c(mean(drywt),drywt[i]))

abline(grandmean,0,lty=2)
text(14,grandmean+10,"Grand Mean")
text(14,mean(wetwt)+10,"mean wetwt")
text(3,mean(drywt)+10,"mean drywt")

# --- D
plot(both, xlab="", ylab="weight", main= "D")
points(rep(mean(wetwt),10), pch=16)
points(rep(mean(drywt),10)~seq(11,20), pch=16)
abline(grandmean,0)
for (i in seq(1,10)) lines( c(i,i), c(grandmean, mean(wetwt)))
for (i in seq(11,20)) lines( c(i,i), c(grandmean, mean(drywt)))
text(15,grandmean+10, "Grand Mean")
#dev.off()

wt = c(wetwt,drywt)
label = as.factor(c(rep("wet",10),rep("dry",10)))
AOV = aov(wt~label)
print (summary(AOV))

cat("\n# -- identical groups --:\n")
wt = c(wetwt,wetwt)
AOV = aov(wt~label)
print (summary(AOV))