###################################################
###################################################
##R tutorial 01

my.url<-'http://www.math.mcgill.ca//dstephens//204//Data//BWdata.txt'

BWdat <- read.table(url(my.url),header=T)

str(BWdat)

?str

names(BWdat)

summary(BWdat$BW)

tapply(BWdat$BW,BWdat$Group,mean)
tapply(BWdat$BW,BWdat$Group,var)

tapply(BWdat$BW,BWdat$Group,summary)

boxplot(BW~Group,data=BWdat)

boxplot(BW~Group,data=BWdat,col="red")
boxplot(BW~Group,data=BWdat,col=c("red","green"))

t.test(BW~Group,data=BWdat,var.equal=T)

##########################################################################
#Small simulation exercise

set.seed(3424)
nreps<-10000

n.A<-10
n.B<-7

mu.A<-3000
mu.B<-3000
var.A<-var.B<-200^2

tstored<-numeric(length=nreps)

#Carry out 10000 repeat experiments
#under the null assumption of equal mus.

for(irep in 1:nreps){

	x.A<-rnorm(n.A,mu.A,sqrt(var.A))
	x.B<-rnorm(n.B,mu.B,sqrt(var.B))

	mean.A<-mean(x.A)
	mean.B<-mean(x.B)
	ssq.A<-var(x.A)
	ssq.B<-var(x.B)
	
	ssq<-((n.A-1)*ssq.A+(n.B-1)*ssq.B)/(n.A+n.B-2)

	tstat<-(mean.A-mean.B)/sqrt(ssq*(1/n.A+1/n.B))

	tstored[irep]<-tstat

}

hist(tstored)

#Theory says that the distribution of the statistic
$should be Student-t(15) - is it ?

del<-0.25
br<-seq(-6,6,by=del)
hist(tstored,breaks=br,xlab="t",main="Histogram of test statistics")

xv<-seq(-6,6,by=0.01)
yv<-dt(xv,df=n.A+n.B-2)

lines(xv,yv*nreps*del,col="red",lwd=2)

#A good match bewteen theory and observed results !

#Alternative plotting method: get hist to plot on the 'density' scale.

del<-0.25
br<-seq(-6,6,by=del)
hist(tstored,breaks=br,xlab="t",main="Histogram of test statistics",freq=FALSE)

xv<-seq(-6,6,by=0.01)
yv<-dt(xv,df=n.A+n.B-2)
lines(xv,yv,col="red",lwd=2)