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 from and article by Yang, Sun, and Hardin for a non-iterative R function for finding Tango's score confidence interval for a difference of proportions with matched pairs (Statistics in Medicine 1998), described in Agresti and Min (Statistics in Medicine, 2005); for details, see https://works.bepress.com/zyang/1/

-------------------------------------------------------------------------------------------------
### R code to calculate the Tango's CI using a non-iterative method;
### "A non-iterative implementation of Tango's score confidence interval for a paired difference of proportions";
###  Statistics in Medicine;
### by Zhao Yang, Xuezheng Sun, James W. Hardin;
### Version Date: 12Jul2012

### b is the cell count in the (Success,Failure) cell, and c the cell count in the (Failure,Success) cell

TangoCI <- function(N, b, c, alpha = 0.05) {
  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)
-------------------------------------------------------------------------------------------------

This is older R code (by Yongyi Min) for Tango's score confidence interval for a difference of proportions with matched pairs, (Statistics in Medicine 1998),
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
-----------------------------------------------------------------------------