if !require(fmsb) { install.packages("fmsb") library(fmsb) } # Figure 8-1 TAB1 <- matrix(c(4,4,386,1250), 2) T8.1 <- pvalueplot(TAB1, plot.OR=TRUE, plot.log=TRUE, xrange=c(0.1, 100)) res <- oddsratio(TAB1) segments(1, 0, 1, 1, lty=1) segments(0.1, res$p.value, 100, res$p.value, lty=2) text(res$estimate, 1, "point estimate") text(1, res$p.value, sprintf("p-value=%3.2f", res$p.value)) # Figure 8-2 TAB2 <- matrix(c(1090, 1000, 14910, 15000), 2) T8.2 <- pvalueplot(TAB2, plot.OR=TRUE, plot.log=TRUE, xrange=c(0.1, 100)) lines(T8.1$OR, T8.1$p.value) segments(1, 0, 1, 1, lwd=2) # Figure 8-3 CI95 <- oddsratio(TAB1, conf.level=0.95)$conf.int CI90 <- oddsratio(TAB1, conf.level=0.90)$conf.int CI80 <- oddsratio(TAB1, conf.level=0.80)$conf.int pvalueplot(TAB1, plot.OR=TRUE, plot.log=TRUE, xrange=c(0.1, 100)) segments(1, 0, 1, 1, lty=1) arrows(CI95[1], 0.05, CI95[2], 0.05, code=3, lty=3) text(res$estimate, 0.05, "95% CI") arrows(CI90[1], 0.10, CI90[2], 0.10, code=3, lty=3) text(res$estimate, 0.10, "90% CI") arrows(CI80[1], 0.20, CI80[2], 0.20, code=3, lty=3) text(res$estimate, 0.20, "80% CI") # Table 8-3 TAB3 <- matrix(c(468, 229, 480, 205), 2) res3 <- oddsratio(TAB3) # https://www.nejm.org/doi/full/10.1056/NEJM199810083391504 # https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(05)74403-2/fulltext ORMH <- function(TBL, conf.level=0.95) { TT <- rowSums(TBL) GG <- TBL[,1]*TBL[,4]/TT HH <- TBL[,2]*TBL[,3]/TT OR <- sum(GG)/sum(HH) PP <- (TBL[,1]+TBL[,4])/TT QQ <- (TBL[,2]+TBL[,3])/TT VARlnOR <- sum(GG*PP)/(2*sum(GG)^2) + sum(GG*QQ+HH*PP)/(2*sum(GG)*sum(HH)) + sum(HH*QQ)/(2*sum(HH)^2) SElnOR <- sqrt(VARlnOR) ORL <- exp(log(OR)-qnorm(1-(1-conf.level)/2)*SElnOR) ORU <- exp(log(OR)+qnorm(1-(1-conf.level)/2)*SElnOR) return(list(estimate=OR, conf.int=c(ORL, ORU), conf.level=conf.level)) } # TBL has to be matrix # each line is one study # columns are exposed cases, unexposed cases, exposed controls, and unexposed controls # Table 10-6 ORMH(matrix(c(3, 9, 104, 1059, 1, 3, 5, 86), 2, 4, byrow=TRUE), conf.level=0.9) TenStudies <- matrix( c(215, 229, 311-215, 306-229, 38, 33, 59-38, 51-33, 161, 174, 293-161, 293-174, 76, 88, 164-76, 163-88, 103, 105, 129-103, 133-105, 65, 67, 120-65, 125-67, 81, 75, 113-81, 110-75, 48, 63, 160-48, 159-63, 22, 21, 60-22, 62-21, 56, 51, 137-56, 140-51 ), 10, 4, byrow=TRUE) ElevenStudies <- rbind(TenStudies, c(468, 480, 229, 205)) ORMH(TenStudies) ORMH(ElevenStudies) pvpORMH <- function (XTAB, xrange = c(0.6, 1.2)) { cp <- c(1:9/1000, 1:9/100, 10:90/100, 0.9 + 1:9/100, 0.99 + 1:9/1000) cl <- 1-cp lu <- uu <- numeric(length(cl)) for (i in 1:length(cl)) { res <- ORMH(XTAB, conf.level=cl[i]) lu[i] <- res$conf.int[1] uu[i] <- res$conf.int[2] } cpx <- c(cp, 1, rev(cp)) cOR <- c(lu, res$estimate, rev(uu)) rval <- data.frame(OR = cOR, p.value = cpx) plot(rval, type = "l", xlim = xrange) return(rval) } R11 <- pvpORMH(ElevenStudies) pvpORMH(TenStudies) lines(R11$OR, R11$p.value, lty=2) segments(1, 0, 1, 1, lwd=2) # Figure 8-4 is made by the code above ################################################### ## Alternative way ################################################### # library(meta) # res <- metabin(TenStudies[,1],TenStudies[,1]+TenStudies[,3],TenStudies[,2],TenStudies[,2]+TenStudies[,4], sm="OR") # print(res) # forest(res) # drapery(res, type="pval") # Draw p-value plots ################################################### # Table 8-4 res4 <- riskratio(12, 5, 59, 50, conf.level=0.9) TAB4 <- matrix(c(12, 47, 5, 45), 2) pvalueplot(TAB4, plot.OR=FALSE, plot.log=TRUE, xrange=c(0.1, 10)) segments(0.1, res4$p.value, 10, res4$p.value, lty=2) text(0.5, res4$p.value, sprintf("p-value=%3.2f", res4$p.value)) segments(1, 0, 1, 1, lwd=2) arrows(res4$conf.int[1], 0.1, res4$conf.int[2], 0.1, code=3, lty=3) text(res4$estimate, 0.1, paste("90%", sprintf("CI: %3.2f-%3.2f", res4$conf.int[1], res4$conf.int[2]), sep="")) # https://jamanetwork.com/journals/jama/fullarticle/193754