next up previous
Next: About this document ...

R code for the adjustment of the Wald confidence interval for a difference of proportions, with matched pairs.
This is the interval called Wald+2 in Agresti and Min, Statistics in Medicine, 2005, which adds
0.5 to each cell before constructing the Wald CI. The CI is truncated when it overshoots the boundary.

--------------------------------------------------------------------------
diffpropci <- function(b,c,n,conflev)
{
  z  <- qnorm(1-(1-conflev)/2)
  diff <- (c-b)/(n+2)
  sd <- sqrt((b+c+1)-(c-b)^2/(n+2))/(n+2)
  ll <- diff - z*sd
  ul <- diff + z*sd
  if(ll < -1) ll = -1
  if(ul > 1) ul = 1 
  c(ll, ul)
}
# Adjusted Wald interval for difference of proportions with matched pairs
# "conflev"=confidence coefficient, n=sample size, b,c = off-diag counts
---------------------------------------------------------------------------


R code for Tango's score confidence interval for a difference of proportions (Statistics in Medicine 1998),
with matched pairs, described in Agresti and Min (Statistics in Medicine, 2005)

-----------------------------------------------------------------------------
scoreci <- function(b,c,n,conflev)
{
   pa = 2*n
   z = qnorm(1-(1-conflev)/2)

   if(c == n) {ul = 1}
   else{ 
     proot = (c-b)/n
     dp = 1-proot
     niter = 1
     while(niter <= 50){
       dp = 0.5*dp
       up2 = proot+dp
       pb = - b - c + (2*n-c+b)*up2
       pc = -b*up2*(1-up2)
       q21 = (sqrt(pb^2-4*pa*pc)-pb)/(2*pa)
       score = (c-b-n*up2)/sqrt(n*(2*q21+up2*(1-up2)))
       if(abs(score)<z){ proot = up2 }
       niter=niter+1
       if((dp<0.0000001) || (abs(z-score)<.000001)){
       niter=51
       ul=up2
       }
      } 
   }

   if(b == n) {ll = -1}
   else{
     proot = (c-b)/n
     dp = 1+proot
     niter = 1
     while(niter <= 50){  
       dp = 0.5*dp
       low2 = proot-dp
       pb = - b - c + (2*n-c+b)*low2
       pc = -b*low2*(1-low2)
       q21 = (sqrt(pb^2-4*pa*pc)-pb)/(2*pa)
       score = (c-b-n*low2)/sqrt(n*(2*q21+low2*(1-low2)))
       if(abs(score) < z){proot = low2}
       niter = niter+1
       if((dp<0.0000001) || (abs(z-score)<.000001)){
       ll = low2
       niter = 51
       }
      }
  }
  c(ll,ul)
}
# Tango score interval for difference of proportions with matched
# pairs, confidence coefficient = conflev, n = sample size,
# b,c = off-diagonal counts
-----------------------------------------------------------------------------


R code for the adapted binomial score confidence interval for the subject-specific odds ratio,
with matched pairs, described in Agresti and Min, Statistics in Medicine, 2004.
This uses the Wilson score CI for a binomial parameter with the off-diagonal counts.

-----------------------------------------------------------------------------
oddsratioci <- function(b,c,conflev)
{
  z  <- qchisq(conflev,1)
  A <- b + c + z
  B <- 2*c + z
  C <- c^2/(b+c)
  l <- (B - sqrt(B^2-4*A*C))/(2*A)
  u <- (B + sqrt(B^2-4*A*C))/(2*A)
  ll <- l/(1-l)
  ul <- u/(1-u)
  c(ll,ul)
}
# adapts Wilson binomial score CI to form CI for subject-specific
# odds ratio with matched pairs data
# "conflev"=confidence coefficient, b,c = off-diag counts
-----------------------------------------------------------------------------





Alan Agresti 2003-11-13