################################################### ### chunk number 1: ################################################### system("rm -rf figs/*.*") library(meta) library(rmeta) library(tools) # chunk=1 ################################################### ### chunk number 2: ################################################### s <- sessionInfo() s$locale <- Sys.getlocale('LC_TIME') cat(toLatex(s), sep='\n') # chunk=2 ################################################### ### chunk number 3: ################################################### rosig <- read.csv("rosiglitazone.csv", header=T, sep=",") rosig[,1:7] # chunk=3 ################################################### ### chunk number 4: ################################################### meta.a <- metabin(e.mi.t, n.t, e.mi.c, n.c, studlab=study, data=rosig, method="Peto", sm="OR", incr=0.0, allstudies=TRUE) summary(meta.a) # chunk=4 ################################################### ### chunk number 5: ################################################### meta.a # chunk=5 ################################################### ### chunk number 6: ################################################### meta.b <- metabin(e.mi.t, n.t, e.mi.c, n.c, studlab=study, data=rosig, method="MH", sm="OR", incr=0.0, allstudies=TRUE) meta.b # chunk=6 ################################################### ### chunk number 7: ################################################### indx <- c(1:42) ind.1 <- rosig$e.mi.t+rosig$e.mi.c # if 0 no event in both arms ind.2 <- rosig$e.mi.t*rosig$e.mi.c # if 0 at least one zero cell n.c <- rosig$n.c[indx<41&ind.1>0] y.c <- rosig$e.mi.c[indx<41&ind.1>0] pi.c <- sum(y.c)/sum(n.c) n.t <- rosig$n.t[indx<41&ind.1>0] y.t <- rosig$e.mi.t[indx<41&ind.1>0] pi.t <- sum(y.t)/sum(n.t) tapply(n.c, y.c==0, summary) tapply(n.t, y.t==0, summary) wilcox.test(n.c~(y.c==0)) wilcox.test(n.t~(y.t==0)) # chunk=7 ################################################### ### chunk number 8: ################################################### z.c <- z.t <- r <- NULL for (i in 1:5000) { y.0 <- rbinom(length(n.c), n.c, pi.c) y.1 <- rbinom(length(n.t), n.t, pi.t) z.c <- c(z.c, sum(y.0==0)) z.t <- c(z.t, sum(y.1==0)) r <- c(r, sum(y.0==0)/sum(y.1==0)) } summary(z.c) summary(z.t) summary(r) # chunk=8 ################################################### ### chunk number 9: ################################################### # function to generate zero-truncated binomial variates rztb <- function(n, size, prob) { # n, prob are scalars, size is vector y <- rep(0, n) if (length(size)==1) { for (i in 1:n) { y[i] <- rbinom(1, size, prob) while(y[i]==0) { y[i] <- rbinom(1, size, prob) } } } else { for (i in 1:n) { y[i] <- rbinom(1, size[i], prob) while(y[i]==0) { y[i] <- rbinom(1, size[i], prob) } } } y } # chunk=9 ################################################### ### chunk number 10: ################################################### # use rosig data to generate some parameter values n.c <- rosig$n.c[indx<41] n.t <- rosig$n.t[indx<41] y.c <- rosig$e.mi.c[indx<41] y.t <- rosig$e.mi.t[indx<41] p.c <- sum(y.c>0)/length(y.c) p.t <- sum(y.t>0)/length(y.t) pi.c <- sum(y.c[y.c>0])/sum(n.c[y.c>0]) r <- 1.5 # true odds ratio pi.t <- pi.c*r N <- length(n.c) nsim <- 1000 out <- NULL OR.Peto <- OR.MH <- OR.Peto.zc <- OR.MH.zc <- NULL for (i in 1:nsim) { z.c <- rbinom(N, 1, p.c) z.t <- rbinom(N, 1, p.t) y.c <- y.t <- rep(0, N) y.c[z.c==1] <- rztb(sum(z.c), n.c[z.c==1], pi.c) y.t[z.t==1] <- rztb(sum(z.t), n.t[z.t==1], pi.t) dat <- data.frame(study=c(1:N), n.t, y.t, n.c, y.c) # Peto method out <- metabin(y.t, n.t, y.c, n.c, studlab=study, data=dat, method="Peto", sm="OR", incr=0.0, allstudies=TRUE) OR.Peto <- c(OR.Peto, exp(summary(out)$fixed$TE)) # MH method out <- metabin(y.t, n.t, y.c, n.c, studlab=study, data=dat, method="MH", sm="OR", incr=0.0, allstudies=TRUE) OR.MH <- c(OR.MH, exp(summary(out)$fixed$TE)) # zero-correction ind.1 <- y.t+y.c # if 0 no event in both arms ind.2 <- y.t*y.c # if 0 at least one zero cell # add 0.5 correction to zero cells dat$n.t[ind.1>0&ind.2==0] <- dat$n.t[ind.1>0&ind.2==0]+1 dat$n.c[ind.1>0&ind.2==0] <- dat$n.c[ind.1>0&ind.2==0]+1 dat$y.t[ind.1>0&ind.2==0] <- dat$y.t[ind.1>0&ind.2==0]+0.5 dat$y.c[ind.1>0&ind.2==0] <- dat$y.c[ind.1>0&ind.2==0]+0.5 # Peto method out <- metabin(y.t, n.t, y.c, n.c, studlab=study, data=dat, method="Peto", sm="OR", incr=0.0, allstudies=TRUE) OR.Peto.zc <- c(OR.Peto.zc, exp(summary(out)$fixed$TE)) # MH method out <- metabin(y.t, n.t, y.c, n.c, studlab=study, data=dat, method="MH", sm="OR", incr=0.0, allstudies=TRUE) OR.MH.zc <- c(OR.MH.zc, exp(summary(out)$fixed$TE)) } summary(OR.Peto) summary(OR.MH) summary(OR.Peto.zc) summary(OR.MH.zc) # chunk=10 ################################################### ### chunk number 11: ################################################### method <- c(rep("Peto",nsim),rep("MH", nsim),rep("Peto-0.5",nsim),rep("MH-0.5",nsim)) OR <- c(OR.Peto, OR.MH, OR.Peto.zc, OR.MH.zc) boxplot(OR~method, border="grey",col="lightgrey", ylab="Risk Ratio", boxwex=0.5) points(jitter(rep(1:4, each=nsim),0.5), unlist(split(OR, method)), cex=0.5, pch='.') abline(h=1.5, col="red") # chunk=11 ################################################### ### chunk number 12: ################################################### n.t <- n.c <- rep(100, 20) y.t <- c(rep(1,16),rep(0,4)) y.c <- c(rep(0,16),rep(4,4)) dat <- data.frame(study=c(1:N), n.t, y.t, n.c, y.c) out <- metabin(y.t, n.t, y.c, n.c, studlab=study, data=dat, method="Peto", sm="OR", incr=0.0, allstudies=TRUE) summary(out) out <- metabin(y.t, n.t, y.c, n.c, studlab=study, data=dat, method="MH", sm="OR", incr=0.0, allstudies=TRUE) summary(out) out <- metabin(y.t, n.t, y.c, n.c, studlab=study, data=dat, method="Peto", sm="OR", incr=0.5, allstudies=TRUE) summary(out) # chunk=12 ################################################### ### chunk number 13: ################################################### # work around for dump bewteen R 2.5.0 and JAGS dump.jags <- function(list, file='dumpdata.R') { # dump in R-2.5.0 has different parsings, see ?.deparseOpts for details dump(list, file='temp', control=c("all","S_compatible")) # may add "S_compatible" fixedfile <- file # 1) change ` to "; 2) add as.integer(c()) on .Dim=n:(n+1); 3) change : to , system(paste("sed -e 's/\`/\"/g' -e 's/[0-9]*\:[0-9]*/as.integer(c(&))/g' -e 's/\:/\, /g' temp >", fixedfile)) file.remove("temp") } ################################################### ### chunk number 14: ################################################### indx <- c(1:42) dat <- rosig ind.1 <- dat$e.mi.t+dat$e.mi.c # if 0 no event in both arms ind.2 <- dat$e.mi.t*dat$e.mi.c # if 0 at least one zero cell k <- 36 r.c <- dat$e.mi.c[indx<41 & ind.1>0] r.t <- dat$e.mi.t[indx<41 & ind.1>0] n.c <- dat$n.c[indx<41 & ind.1>0] n.t <- dat$n.t[indx<41 & ind.1>0] datta <- c("k","r.t","n.t","r.c","n.c") delta <- rep(1.43, k) theta <- 1.43 tau <- 1 pars <- c("delta","theta","tau") # dump.jags(datta, "meta.dat") dump.jags(pars, "meta.in") system("jags OR-a.cmd") ################################################### ### chunk number 15: ################################################### require(coda) out <- read.jags("rosig.out") summary(out) summary(exp(out[,1])) ################################################### ### chunk number 16: ################################################### indx <- c(1:42) dat <- rosig ind.1 <- dat$e.mi.t+dat$e.mi.c # if 0 no event in both arms ind.2 <- dat$e.mi.t*dat$e.mi.c # if 0 at least one zero cell # confirm the denominators for combined small studies, 10285 and 6106 y <- dat$n.t[ind.1>0&indx<41] sum(y) y <- dat$n.c[ind.1>0&indx<41] sum(y) # add 0.5 correction to zero cells dat <- rosig dat$n.t[ind.1>0&ind.2==0] <- dat$n.t[ind.1>0&ind.2==0]+1 dat$n.c[ind.1>0&ind.2==0] <- dat$n.c[ind.1>0&ind.2==0]+1 dat$e.mi.t[ind.1>0&ind.2==0] <- dat$e.mi.t[ind.1>0&ind.2==0]+0.5 dat$e.mi.c[ind.1>0&ind.2==0] <- dat$e.mi.c[ind.1>0&ind.2==0]+0.5 # chunk=12 ################################################### ### chunk number 17: ################################################### meta.1 <- metabin(e.mi.t, n.t, e.mi.c, n.c, studlab=study, subset=(ind.1[indx<41]>0), data=dat[indx<41,],method="Peto", sm="OR",incr=0,allstudies=TRUE) summary(meta.1) # chunk=13 ################################################### ### chunk number 18: ################################################### meta.2 <- metabin(e.mi.t, n.t, e.mi.c, n.c, studlab=study,subset=(ind.1>0), data=dat,method="Peto",sm="OR",incr=0.0, allstudies=TRUE) meta.2 # chunk=14 ################################################### ### chunk number 19: ################################################### meta.3 <- metabin(e.mi.t, n.t, e.mi.c, n.c, studlab=study,subset=(ind.1>0), data=dat,method="MH",sm="OR",incr=0.0, allstudies=TRUE) summary(meta.3) # chunk=15 ################################################### ### chunk number 20: ################################################### # fixed effects model a <- meta.MH(n.t,n.c,e.mi.t,e.mi.c, data=dat,names=study,subset=(ind.1>0)) summary(a) # random effect model b <- meta.DSL(n.t,n.c,e.mi.t,e.mi.c, data=dat,names=study,subset=(ind.1>0)) summary(b) # chunk=16 ################################################### ### chunk number 21: ################################################### plot(a) #, colors=meta.colors(summary="black", lines="#FFFFF0", # box="#FFFF50",zero="grey90", text="black",background="grey", # axes="grey90")) # chunk=17 ################################################### ### chunk number 22: ################################################### meta.c <- metabin(e.death.t, n.t, e.death.c, n.c, studlab=study, data=rosig[indx<41,], method="Peto", sm="OR", incr=0.0, allstudies=TRUE) summary(meta.c) # chunk=17 ################################################### ### chunk number 23: ################################################### meta.d <- metabin(e.death.t, n.t, e.death.c, n.c, studlab=study, data=rosig[,], method="Peto", sm="OR", incr=0.0, allstudies=TRUE) summary(meta.d) # chunk=18 ################################################### ### chunk number 24: ################################################### dat <- rosig ind.1 <- dat$e.death.t+dat$e.death.c ind.2 <- dat$e.death.t*dat$e.death.c sum(dat$n.t[ind.1>0&indx<41]) sum(dat$n.c[ind.1>0&indx<41]) dat$n.t[ind.1>0&ind.2==0] <- dat$n.t[ind.1>0&ind.2==0]+1 dat$n.c[ind.1>0&ind.2==0] <- dat$n.c[ind.1>0&ind.2==0]+1 dat$e.death.t[ind.1>0&ind.2==0] <- dat$e.death.t[ind.1>0&ind.2==0]+0.5 dat$e.death.c[ind.1>0&ind.2==0] <- dat$e.death.c[ind.1>0&ind.2==0]+0.5 # chunk=19 ################################################### ### chunk number 25: ################################################### meta.4 <- metabin(e.death.t, n.t, e.death.c, n.c, studlab=study, subset=(ind.1[indx<41]>0), data=dat[indx<41,],method="Peto", sm="OR",allstudies=TRUE) summary(meta.4) # chunk=20 ################################################### ### chunk number 26: ################################################### meta.5 <- metabin(e.death.t, n.t, e.death.c, n.c, studlab=study,subset=(ind.1>0), data=dat,method="Peto",sm="OR", incr=0.0, allstudies=TRUE) meta.5 # chunk=21