## R-code for HW #3
## AHW3
## Feb 20, 2008
rm(list = ls())

y1 <- c(16.03, 16.01, 16.04, 15.96, 16.05, 15.98,
	   16.05, 16.02, 16.02, 15.99)
y2 <- c( 16.02, 16.03,  15.97, 16.04,  15.96, 16.02,
	    16.01, 16.01, 15.99, 16.00)
n1 <- length(y1)
n2 <- length(y2)

y <- c(y1,y2)
x <- c(rep("1",n1), rep("2", n2))
		
## randomization test
t.stat.obs <- t.test(y[x=="1"],y[x=="2"],var.equal=T)$stat

t.stat.sim<-real()
for(nsim in 1:5000) {
xsim<-sample(x)
t.stat.sim[nsim]<- t.test(y[xsim=="1"],y[xsim=="2"],var.equal=T)$stat
}

p.value <-mean( abs(t.stat.sim) >= abs(t.stat.obs) )
p.value
	
					

## diagnostic plots
pdf("diag.pdf",family="Times",height=3,width=6)
par(mar=c(3,3,1,1), mfrow=c(2,2))

qqnorm(y1,main="Normal Q-Q plot for y1"); qqline(y1)
qqnorm(y2,main="Normal Q-Q plot for y2"); qqline(y2)
boxplot(y1, col="yellow",main="y1")
boxplot(y2, col="yellow", main="y2")

dev.off()



## t-test
alpha <- 0.05
sp.sq <- ((n1-1)/(n1+n2-2))*var(y1) + ((n2-1)/(n1+n2-2))*var(y2)
sp <- sqrt(sp.sq)
t.obs <- (mean(y1)-mean(y2))/(sp*sqrt((1/n1 + 1/n2)))
t.crit <- qt(1-alpha/2, n1+n2-2)
p.value <- 2*(1-pt(t.obs, n1+n2-2))
p.value

## or
t.test(y[x=="1"],y[x=="2"],var.equal=T)

## 95% CI
(mean(y1) - mean(y2)) - sp*sqrt(1/n1 + 1/n2)*qt(1-0.05/2, n1+n2-2)
(mean(y1) - mean(y2)) + sp*sqrt(1/n1 + 1/n2)*qt(1-0.05/2, n1+n2-2)

## 
x <- seq(-8,8,0.01)
plot(x,dt(x,n1+n2-2),type="l", lwd=2, main="H0 distribution")		
abline(v=c(-t.crit, t.crit), lwd=2, col="red") 
abline(v=c(-abs(t.obs), abs(t.obs)), lwd=2, col="blue")


## power
delta <- 0.5
alpha <- 0.05
sigma <- sp

# NULL
pdf("power.pdf",family="Times",height=3,width=6)
par(mar=c(3,3,1,1))
x <- seq(-75,75,0.01)
plot(x,dt(x,n1+n2-2),type="l", lwd=2, main="Null and Alternative Distributions with Specific delta=0.5")	
abline(v=c(-t.crit, t.crit), lwd=2, col="red") 

## Alternative
t.gamma <- delta/(sigma*sqrt(1/n1+1/n2))
lines(x,dt(x,n1+n2-2, ncp=t.gamma),type="l", lwd=2)
dev.off()

power <- pt(-t.crit, n1+n2-2, ncp=t.gamma) + 1-pt(t.crit, n1+n2-2, ncp=t.gamma)
power