#########################################################
#########################################################
##########                                     ##########
##########	R Scripts for MATH 204 Winter 2013 ##########
##########                                     ##########
#########################################################
#########################################################

#Understanding the ANOVA F-statistic

n<-18;ni<-rep(6,3)
k<-3
X<-c(-0.88,0.24,-0.46,0.78,-0.47,-0.38,
-0.75,0.11,0.64,1.98,-1.03,1.84,
1.38,1.20,0.42,0.05,-1.29,-0.04)

tmt<-factor(rep(1:3,each=6),labels=c("TMT1","TMT2","TMT3"))

#Compute treament means
tapply(X,tmt,mean)

#Compute treatment variances
tapply(X,tmt, var)

#Compute overall mean
mean(X)

#Compute pooled variance
sum((ni-1)*tapply(X,tmt,var))/(n-k)

#Compute SST
SST<-sum(ni*(tapply(X,tmt,mean)-mean(X))^2)

#Compute SSE
SSE<-sum((ni-1)*tapply(X,tmt,var))

#Compute SS
SS<-sum((n-1)*var(X))

#Compute F
F<-(SST/(k-1))/(SSE/(n-k))

#Compute Critical Value
alpha<-0.05
CR<-qf(1-alpha,k-1,n-k)

#Compute the p-value
pval<-1-pf(F,k-1,n-k)

#Plot the null distribution
xv<-seq(0,5,by=0.01)
yv<-df(xv,k-1,n-k)
plot(xv,yv,type="l",xlab="x",ylab="f(x)")
abline(v=CR,col="red")
text(F,0,"x")

#Means plot
library(gplots)
plotmeans(X~tmt)

#########################################
#Second Data Set
X<-c(-0.88,0.24,-0.46,0.78,-0.47,-0.38,
9.25,10.11,10.64,11.98,8.97,11.84,
21.38,21.20,20.42,20.05,18.71,19.96)

#Compute treament means
tapply(X,tmt,mean)

#Compute treatment variances
tapply(X,tmt, var)

#Compute overall mean
mean(X)

#Compute pooled variance
sum((ni-1)*tapply(X,tmt,var))/(n-k)

#Compute SST
SST<-sum(ni*(tapply(X,tmt,mean)-mean(X))^2)

#Compute SSE
SSE<-sum((ni-1)*tapply(X,tmt,var))

#Compute SS
SS<-sum((n-1)*var(X))

#Compute F
F<-(SST/(k-1))/(SSE/(n-k))

#Compute Critical Value
alpha<-0.05
CR<-qf(1-alpha,k-1,n-k)

#Compute the p-value
pval<-1-pf(F,k-1,n-k)


#Plot the null distribution and test statistoc
xv<-seq(0,5,by=0.01)
yv<-df(xv,k-1,n-k)
plot(xv,yv,type="l",xlab="x",ylab="f(x)")
abline(v=CR,col="red")
text(F,0,"x")

#Means plot
require(gplots)
plotmeans(X~tmt)


##################################################################
##################################################################
#Diet vs Protein Level

boxplot(protein~Diet,data=Milk,xlab="Diet",ylab="Protein Level (%)")

#A summary function
my.summary<-function(x){

	my.sum<-data.frame(n=length(x),Mean=mean(x),
                         StDev=sd(x),StdErr=sd(x)/sqrt(length(x)),
                         Lower=mean(x)-1.96*sd(x)/sqrt(length(x)),
                         Upper=mean(x)+1.96*sd(x)/sqrt(length(x)),
                         Min=min(x),Max=max(x))
	return(round(my.sum,5))

}

tapply(Milk$protein,Milk$Diet,my.summary)
my.summary(Milk$protein)

#Levene's Test for Equality of variances

library(car)
leveneTest(protein~Diet,center='mean',data=Milk)
leveneTest(protein~Diet,center='median',data=Milk)

#ANOVA

#Total
summary(aov(protein~1,data=Milk))
anova(lm(protein~1,data=Milk))

#Fit and Error (Residual)
summary(aov(protein~Diet,data=Milk))
anova(lm(protein~Diet,data=Milk))
summary(lm(protein~Diet,data=Milk))

#Means plot
require(gplots)
plotmeans(protein~Diet,data=Milk)

##################################################################
##################################################################
#Pull strength

boxplot(Pull~Dose,data=Pull,xlab="Dose level",ylab="Pull strength")

tapply(Pull$Pull,Pull$Dose,my.summary)
my.summary(Pull$Pull)

#Levene's Test for Equality of variances

library(car)
leveneTest(Pull~Dose,center='mean',data=Pull)
leveneTest(Pull~Dose,center='median',data=Pull)

#ANOVA

#Total
summary(aov(Pull~1,data=Pull))
anova(lm(Pull~1,data=Pull))

#Fit and Error (Residual)
summary(aov(Pull~Dose,data=Pull))
anova(lm(Pull~Dose,data=Pull))
summary(lm(Pull~Dose,data=Pull))

#Means plot
require(gplots)
plotmeans(Pull~Dose,data=Pull)

##################################################################
##################################################################
#Mannitol gut permeability

boxplot(mann~diag,data=AIDS,xlab="Diagnosis",ylab="Permeability")

tapply(AIDS$mann,AIDS$diag,my.summary)
my.summary(AIDS$mann)

#Levene's Test for Equality of variances

library(car)
leveneTest(mann~diag,center='mean',data=AIDS)
leveneTest(mann~diag,center='median',data=AIDS)

#ANOVA

#Total
summary(aov(mann~1,data=AIDS))
anova(lm(mann~1,data=AIDS))

#Fit and Error (Residual)
summary(aov(mann~diag,data=AIDS))
anova(lm(mann~diag,data=AIDS))
summary(lm(mann~diag,data=AIDS))

#Means plot
require(gplots)
plotmeans(mann~diag,data=AIDS)

##################################################################
##################################################################
#Milk bacteria level

boxplot(BacLevel~Batch,data=MilkBac,xlab="Batch",ylab="bacteria Level (per million)")

tapply(MilkBac$BacLevel,MilkBac$Batch,my.summary)
my.summary(MilkBac$BacLevel)

#Levene's Test for Equality of variances

library(car)
leveneTest(BacLevel~Batch,center='mean',data=MilkBac)
leveneTest(BacLevel~Batch,center='median',data=MilkBac)

#ANOVA

#Total
summary(aov(BacLevel~1,data=MilkBac))
anova(lm(BacLevel~1,data=MilkBac))

#Fit and Error (Residual)
summary(aov(BacLevel~Batch,data=MilkBac))
anova(lm(BacLevel~Batch,data=MilkBac))
summary(lm(BacLevel~Batch,data=MilkBac))

#Means plot
require(gplots)
plotmeans(BacLevel~Batch,data=MilkBac,ylab="bacteria Level (per million)")

##################################################################
##################################################################
#PTSD Treatment

boxplot(PTSD~GROUP,data=PTSD,xlab="GROUP",ylab="PTSD Score")

tapply(PTSD$PTSD,PTSD$GROUP,my.summary)
my.summary(PTSD$PTSD)

#Levene's Test for Equality of variances

library(car)
leveneTest(PTSD~GROUP,center='mean',data=PTSD)
leveneTest(PTSD~GROUP,center='median',data=PTSD)

#ANOVA

#Total
summary(aov(PTSD~1,data=PTSD))
anova(lm(PTSD~1,data=PTSD))

#Fit and Error (Residual)
summary(aov(PTSD~GROUP,data=PTSD))
anova(lm(PTSD~GROUP,data=PTSD))
summary(lm(PTSD~GROUP,data=PTSD))

#Means plot
require(gplots)
plotmeans(PTSD~GROUP,data=PTSD,xlab="GROUP",ylab="PTSD Score")


##################################################################
##################################################################
#Memory Task

boxplot(Memory~Task,data=Memory,xlab="Task",ylab="Memory Score")

tapply(Memory$Memory,Memory$Task,my.summary)
my.summary(Memory$Memory)

#Levene's Test for Equality of variances

library(car)
leveneTest(Memory~Task,center='mean',data=Memory)
leveneTest(Memory~Task,center='median',data=Memory)

#ANOVA

#Total
summary(aov(Memory~1,data=Memory))
anova(lm(Memory~1,data=Memory))

#Fit and Error (Residual)
summary(aov(Memory~Task,data=Memory))
anova(lm(Memory~Task,data=Memory))
summary(lm(Memory~Task,data=Memory))

#Means plot
require(gplots)
plotmeans(Memory~Task,data=Memory,xlab="Task",ylab="Memory Score")

##################################################################
##################################################################
#Post-hoc tests

summary(aov(Pull~Dose,data=Pull))

pairwise.t.test(Pull$Pull, Pull$Dose, p.adj = "none")
pairwise.t.test(Pull$Pull, Pull$Dose, p.adj = "bonferroni")

TukeyHSD(aov(Pull~Dose,data=Pull))

plot(TukeyHSD(aov(Pull~Dose,data=Pull)))



summary(aov(PTSD~GROUP,data=PTSD))

pairwise.t.test(PTSD$PTSD, PTSD$GROUP, p.adj = "none")
pairwise.t.test(PTSD$PTSD, PTSD$GROUP, p.adj = "bonferroni")

TukeyHSD(aov(PTSD~GROUP,data=PTSD))

plot(TukeyHSD(aov(PTSD~GROUP,data=PTSD)))