## 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