### R code to calculate the Tango's CI using a non-iterative method; ### Code to create the results in Table 2 of the manuscript; ### "A non-iterative implementation of Tango's score confidence interval for a paired difference of proportions"; ### for Statistics in Medicine; ### by Zhao Yang, Xuezheng Sun, James W. Hardin; ### Version Date: 12Jul2012 TangoCI <- function(N, b, c, alpha = 0.05) { options(digits = 12) g <- qnorm(1 - alpha/2)^2 m0 <- N^2*(N/g + 1)^2 m1 <- -N*(b - c)*(N/g + 1)*(4*N/g + 1) m2 <- 2*N*(b - c)^2*(3*N/g + 2)/g - N^2*((b + c)/g + 1) m3 <- N*(b - c)*(2*(b + c)/g + 1) - (b - c)^3*(4*N/g + 1)/g m4 <- (b - c)^2*((b - c)^2/g - (b + c) )/g u1 <- m1/m0; u2 <- m2/m0; u3 <- m3/m0; u4 <- m4/m0; if (b != c){ nu1 <- -u2 nu2 <- u1*u3 - 4*u4 nu3 <- -(u3^2 + u1^2*u4 - 4*u2*u4) d1 <- nu2 - nu1^2/3 d2 <- 2*nu1^3/27 - nu1*nu2/3 + nu3 CritQ <- d2^2/4 + d1^3/27 y1A <- h1 <- h2 <- h3 <- h4 <- root1 <- root2 <- root3 <- root4 <- c() if (CritQ > 0) { ## keep one real root; BigA <- -d2/2 + sqrt(CritQ) BigB <- -d2/2 - sqrt(CritQ) x1 <- sign(BigA)*abs(BigA)^(1/3) + sign(BigB)*abs(BigB)^(1/3) y1A <- x1 - nu1/3 } if (CritQ == 0) { ## keep two real roots; BigA <- -d2/2 + sqrt(CritQ) BigB <- -d2/2 - sqrt(CritQ) Omega <- complex(real = -1/2, imaginary = sqrt(3)/2) Omega2 <- complex(real = -1/2, imaginary = -sqrt(3)/2) x1 <- sign(BigA)*abs(BigA)^(1/3) + sign(BigB)*abs(BigB)^(1/3) x2 <- Omega*sign(BigA)*abs(BigA)^(1/3) + Omega2*sign(BigB)*abs(BigB)^(1/3) y1A[1] <- x1 - nu1/3 y1A[2] <- x2 - nu1/3 } if (CritQ <0) { ## keep three real roots; BigA <- -d2/2 + sqrt(as.complex(CritQ)) BigB <- -d2/2 - sqrt(as.complex(CritQ)) Omega <- complex(real = -1/2, imaginary = sqrt(3)/2) Omega2 <- complex(real = -1/2, imaginary = -sqrt(3)/2) x1 <- BigA^(1/3) + BigB^(1/3) x2 <- Omega*BigA^(1/3) + Omega2*BigB^(1/3) x3 <- Omega2*BigA^(1/3) + Omega*BigB^(1/3) y1A[1] <- x1 - nu1/3 y1A[2] <- x2 - nu1/3 y1A[3] <- x3 - nu1/3 } y1 <- Re(y1A) #keep the real part; ny <- length(y1) for (i in 1:ny){ h1[i] <- sqrt(u1^2/4 - u2 + y1[i]) h2[i] <- (u1 * y1[i]/2 - u3)/(2*h1[i]) h3[i] <- (u1/2 + h1[i])^2 - 4*(y1[i]/2 + h2[i]) h4[i] <- (h1[i] - u1/2)^2 - 4*(y1[i]/2 - h2[i]) if (h3[i] >= 0){ root1[i] <- ( - (u1/2 + h1[i]) + sqrt( h3[i] ) )/2 root2[i] <- ( - (u1/2 + h1[i]) - sqrt( h3[i] ) )/2 } if (h4[i] >= 0){ root3[i] <- ( (h1[i] - u1/2) + sqrt( h4[i] ) )/2 root4[i] <- ( (h1[i] - u1/2) - sqrt( h4[i] ) )/2 } } lower <- max(-1,min(root1, root2, root3, root4, na.rm = TRUE)) upper <- min(max(root1, root2, root3, root4, na.rm = TRUE), 1) if (b == N & c == 0) root <- c(lower, 1) else if (b == 0 & c == N) root <- c(-1, upper) else root <- c(lower, upper) } if (b == c){ root1 <- -sqrt(-u2) root2 <- sqrt(-u2) root <- c(root1, root2) } return(root) } TangoCI(44,0,1) TangoCI(14,3,1) TangoCI(32,9,3) TangoCI(50,12,2) TangoCI(50,14,0) TangoCI(100,97,1) TangoCI(30,29,1) TangoCI(100,98,0) TangoCI(30,30,0) TangoCI(54,0,0) TangoCI(350,254,2) TangoCI(350,297,3) TangoCI(605,290,73) TangoCI(350,101,220) TangoCI(10,9,0) TangoCI(10,1,0) TangoCI(5,4,0) TangoCI(5,1,0) TangoCI(3,2,0) TangoCI(3,1,0) TangoCI(3,1,1) TangoCI(2,1,0) TangoCI(10,0,0) TangoCI(9,0,0) TangoCI(7,0,0) TangoCI(5,0,0) TangoCI(3,0,0) TangoCI(2,0,0)