## Stat 642 - Feb 15, 2005 
## R - code for Homework 1

#problem 1
D <- data.frame(start = c(22,25,27,31,32,35,37,40,42,46),
                end=c(37,35,42,38,47,49,48,55,57,53),
                status=c(0,1,0,0,0,1,1,0,1,0))
D
time <- apply(pmax(outer(D$end, 2:5*10+10, pmin) -
                    outer(D$start, 2:5*10, pmax),0), 2, sum)
events <- tapply(D$status, cut(D$end, 2:6*10), sum)
events[is.na(events)] <- 0

lambda <- events/time
lambda*100
sqrt(events)/time*100
cumsum(lambda*10)  ## Cumulative hazard
1-exp(-cumsum(lambda*10))  ## P_j

w <- c(30,35,20,15)
wrate <- c(1,2,3,7)
sum(lambda*w)              ## Directly Standardized cum incidence
sum(wrate*w)/100

sum(lambda*w)/(sum(wrate*w)/100)  ## CMF
sum(events)/(sum(wrate*time)/100)  ## SMR
events2 <- c(0,2,1,1)
time2 <- c(30,90,24,15)
events2/time2*100
sum(events2/time2*w)/(sum(wrate*w)/100)  ## CMF 2
sum(events2)/(sum(wrate*time2)/100)      ## SMR 2


## problem 2
matrix(c(0,0,3,4), 2)  # generate example matrix
fisher.test(matrix(c(0,0,3,4), 2))
fisher.test(matrix(c(3,0,0,4), 2))
names(fisher.test(matrix(c(3,0,0,4), 2)))
fisher.test(matrix(c(3,0,0,4), 2))$p.value

fe.test <- matrix(NA, 4,5)
for(x in 0:3) for(y in 0:4)
  fe.test[x+1,y+1] <- fisher.test(matrix(c(x,y,3-x,4-y), 2))$p.value

## alternatively
fe.test <- outer(0:3,0:4,function(x,y)
                {X <- array(c(x,y,3-x,4-y), c(length(x), 2,2))
                  apply(X, 1, function(U) fisher.test(U)$p.value)
                 })
fe.test
dimnames(fe.test) <- list(0:3,0:4)
fe.test

fe.test < 0.05  # which tables are rejected

pc.test <- outer(0:3,0:4,function(x,y)
                {z <- x
                 for(i in seq(along=x)) z[i] <-
                   chisq.test(matrix(c(x[i],y[i],3 - x[i],4-y[i]),2),
                              correct=F)$statistic
                 z})
pc.test[is.na(pc.test)] <- 0

## alternatively - vectorized version
pchisq.stat <- function(x,y,n1,n0) {
  m1 <- x + y ; N <- n1 + n0
  ifelse(m1==0 | m1==N, 0, (x*N - n1*m1)^2*N/n1/(N-n1)/m1/(N-m1))}
  
pc.test <- outer(0:3,0:4, pchisq.stat, n1=3,n0=4)
dimnames(pc.test) <- list(0:3,0:4)
pc.test

outer(dbinom(0:3,3,.1), dbinom(0:4,4,.1))  ## table probabilities
outer(dbinom(0:3,3,.1), dbinom(0:4,4,.1))*(fe.test<.05)  
sum(outer(dbinom(0:3,3,.1), dbinom(0:4,4,.1))*(fe.test<.05))  ## rejection probability
for(p in c(.1,.3,.5)) cat(p, sum(outer(dbinom(0:3,3,p),
     dbinom(0:4,4,p))*(fe.test<.05)),"\n")
for(p in c(.1,.3,.5)) cat(p, sum(outer(dbinom(0:3,3,p),
     dbinom(0:4,4,p))*(pc.test>3.84)),"\n")


#problem 4
DH.table <- rbind(c(7,49,516,445,299,41),c(61,91,615,408,162,20))
DH.table
DH.table[,-1]  ## chop off first column
DH.OR <- DH.table[1,-1]/DH.table[2,-1]*DH.table[2,1]/DH.table[1,1]
DH.OR <- apply(DH.table[,-1],2,function(x,y) x[1]/x[2]*y[2]/y[1],DH.table[,1])
DH.OR
DH.se <- sqrt(apply(1/DH.table[,-1],2,sum) + sum(1/DH.table[,1]))*DH.OR
DH.se
DH.se.log <- apply(DH.table[,-1],2,function(x,y) sqrt(sum(1/c(x,y))),DH.table[,1])
DH.OR*exp(outer(DH.se.log,c(-1,1))*1.96)  ## 95% CI's
apply(DH.table[,-1],2,function(x,y)
      chisq.test(cbind(x,y),correct=F)$statistic, y=DH.table[,1])

##alternatively use pchisq.stat defined above
pchisq.stat(DH.table[1,-1], DH.table[2,-1], DH.table[1,-1]+DH.table[1,1], DH.table[2,-1]+DH.table[2,1])

log(DH.OR)^2/DH.se.log^2