# This function is already included in fmsb package. # Since version 0.4.3, the method of calculation has largely been changed. pvalueplot <- function(XTAB, plot.OR=FALSE, plot.log=FALSE, xrange=c(0.01,5)) { # XTAB must be 2x2 cross table. # ref. Rothman KJ (2002) Epidemiology: An introduction. Oxford Univ. Press # ref. Rothman KJ (2012) Epidemiology: An introduction. 2nd Ed. Oxford Univ. Press # According to 2nd ed. p.156, p-value function must match with nested confidence intervals. # Limitation: p-values less than 0.0005 are not calculated. x.a <- XTAB[1,1] x.b <- XTAB[1,2] x.c <- XTAB[2,1] x.d <- XTAB[2,2] x.N1 <- sum(XTAB[,1]) x.N0 <- sum(XTAB[,2]) cp <- c(1:9/1000, 1:9/100, 10:90/100, 0.9+1:9/100, 0.99+1:9/1000) cpx <- c(cp, 1, rev(cp)) cpy <- c(cp/2, 0.5, 0.5+cp/2) cRR <- exp(log(x.a*x.N0/x.b/x.N1)+qnorm(cpy)*sqrt(1/x.a-1/x.N1+1/x.b-1/x.N0)) cOR <- exp(log(x.a*x.d/x.b/x.c)+qnorm(cpy)*sqrt(1/x.a+1/x.b+1/x.c+1/x.d)) if (plot.OR) { rval <- data.frame(OR=cOR, p.value=cpx) } else { rval <- data.frame(RR=cRR, p.value=cpx) } OpLog <- ifelse(plot.log, "x", "") plot(rval, type="l", xlim=xrange, log=OpLog) return(rval) } XTAB <- matrix(c(12, 5, 47, 45), 2, byrow=TRUE) pvalueplot(XTAB, plot.log=TRUE, xrange=c(0.1, 10))