######################################################################
# Computer notes                                 Biostatistics 140.776
# Homework 2                                  Johns Hopkins University 
######################################################################
# Lines beginning with the symbol '#' are comments in R.  All other
# lines contain code.
#
# In R for Windows, you may wish to open this file from the menu bar
# (File:Display file); you can then easily copy commands into the
# command window.  (Use the mouse to highlight one or more lines; then
# right-click and select "Paste to console".)
######################################################################

###############################################
# 1
###############################################

x=sample(c("A","C","G","T"),10000,rep=T)

ktuple=function(x,k=3){
  if (length(x)==1) x=unlist(strsplit(x,""))
  if (length(x)<k) stop("k is larger than the length of the input string")
  n=length(x)
  z=matrix(0,nrow=k,ncol=n-k+1)
  for (j in 1:k) z[j,]=x[(k:n)+j-k]
  return(table(apply(z,2,paste,collapse="")))}

ktuple(x)
ktuple(x,2)

###############################################
# 2
###############################################

mylm=function(x,y,hatmat=F,addint=T){
  if(addint) x=cbind(1,x)
  if(det(t(x)%*%x)==0) stop("The design matrix is rank deficient!")
  xtxinv=solve(t(x)%*%x)
  beta=as.vector(xtxinv%*%t(x)%*%y)
  fitted=x%*%beta
  residuals=as.vector(y-fitted)
  rss=sum(residuals^2)
  sigma=sqrt(rss/(nrow(x)-ncol(x)))
  varbeta=sigma^2*xtxinv
  hat=x%*%xtxinv%*%t(x)
  if(hatmat){
  z=list(beta=beta,sigma=sigma,varbeta=varbeta,fitted=fitted,
         residuals=residuals,hat=hat)}
  else{
  z=list(beta=beta,sigma=sigma,varbeta=varbeta,fitted=fitted,
         residuals=residuals)}
  return(z)}

install.packages(VR")
library(MASS)
# note: the data frame hills is already loaded!
attach(hills)
summary(lm(time~dist+climb))
mylm(cbind(dist,climb),time)
detach(hills)

###############################################
# 3
###############################################

ct=rnorm(50,150,15)
tr=rnorm(50,140,15)

t.test(ct,tr)
names(t.test(ct,tr))
t.test(ct,tr)$p.value

nit=10000;pv=rep(0,nit)
for (i in 1:nit){
  ct=rnorm(50,150,15)
  tr=rnorm(50,140,15)
  pv[i]=t.test(ct,tr)$p.value}
mean(pv<0.05)

power.t.test(n=50,delta=10,sd=15)
# or:
power.t.test(n=50,delta=2/3)

nit=10000;pv=rep(0,nit)
for (i in 1:nit){
  ct=rnorm(50,150,15)
  tr=rnorm(50,140,20)
  pv[i]=t.test(ct,tr)$p.value}
mean(pv<0.05)

###############################################
# 4
###############################################

n=1000;m=5000;x=matrix(10+rnorm(n*m),c(n,m));z=x
system.time(for (i in 1:m){z[,i]=x[,i]-mean(x[,i])})
system.time(t(apply(x,1,"-",apply(x,2,mean))))
system.time(sweep(x,2,apply(x,2,mean)))
system.time(sweep(x,2,rowsum(x,rep(1,n))/n))
# try also
system.time(apply(x,2,function(x) x-mean(x)))
system.time(sweep(x,2,colMeans(x)))

###############################################
# 5
###############################################

data(airquality)
attach(airquality)

par(las=1)
plot(Ozone,type="h",xaxt="n",xlab="",ylab="Ozone concentration [ppb]",ylim=c(0,175))
axis(1,subset(1:153,Day==1),paste(month.abb[5:9],"1",sep=" "))
axis(4,c(120,140,160),cex.axis=0.7)
abline(h=c(120,140,160),col=c("yellow3","orange","red"),lty=2)
wh=subset(1:153,Ozone>=120)
text(wh,Ozone[wh]+5,paste(month.abb[Month][wh],Day[wh],sep=" "))

plot(Solar.R,Ozone,xlab="Solar radiation",ylab="Ozone concentration [ppb]")

coplot(Ozone~Solar.R|Wind)
coplot(Ozone~Solar.R|Temp)
coplot(Ozone~Solar.R|Wind*Temp)

library(lattice)
wind=equal.count(Wind,number=6)
temp=equal.count(Temp,number=6)

trellis.device(col=F)
xyplot(Ozone~Solar.R|wind)
xyplot(Ozone~Solar.R|temp)
xyplot(Ozone~Solar.R|wind*temp)

##################
# End of hw2.R
##################