#####################################
### Computes confidence Interval  ###
### for the difference of         ###
### two proportions based on      ###
### inverting the score test      ###
### of Mee or MN                 ###
#####################################
score.z <- function(y,n,delta.null,MN=TRUE,alternative = "two.sided") {
  N=sum(n)
  C=sum(y)
  L3=N
  L2=(n[1]+2*n[2])*delta.null-N-C
  L1=(n[2]*delta.null-N-2*y[2])*delta.null+C
  L0=y[2]*delta.null*(1-delta.null)
  c=L2^3/(3*L3)^3 - L1*L2/(6*L3^2) + L0/(2*L3)
  b=ifelse(c>=0,1,-1)*sqrt(L2^2/(3*L3)^2 - L1/(3*L3))
  d <- min(c/b^3,1)
  d <- max(d, -1)
  a=(3.14159265358979+acos(d))/3
  p2=2*b*cos(a)-L2/(3*L3)
  p2 = max(p2, 0)
  p2 = min(p2,1)
  p1=p2+delta.null
  p1 = max(p1, 0)
  p1 = min(p1,1)
  var0 <- p1*(1 - p1)/n[1] + p2*(1 - p2)/n[2]
  if (!MN) z2=(y[1]/n[1]-y[2]/n[2]-delta.null)^2/var0 #Mee
  else z2=(y[1]/n[1]-y[2]/n[2]-delta.null)^2/(var0*N/(N-1)) #Mietinnen and Nurminen
  if (var0 == 0) z2 = 0
  return(round(z2,7))
}

#computes lower or upper score bounds for diff of prop p1 - p2
score.int <- function(y,n,conflev=0.95, type="two.sided", MN=TRUE) {
  if (type=="two.sided") c=qnorm(1-(1-conflev)/2)^2 else c=qnorm(conflev)^2
  delta1=(y[1]+1)/(n[1]+2) - (y[2]+1)/(n[2]+2) #starting point
  if (any(type=="lower",type=="two.sided")) {
    delta2=-1
    while( abs(delta1-delta2)>10^(-7) ) {#Bisection for LB
      delta=(delta1+delta2)/2
      z=score.z(y,n,delta, MN)
      if (z>c) delta2=delta else delta1=delta
    }
  }
  delta.l=delta1
  delta1=(y[1]+1)/(n[1]+2) - (y[2]+1)/(n[2]+2) #starting point
  if (any(type=="upper",type=="two.sided")) {
    delta2=1;
    while( abs(delta1-delta2)>10^(-7) ) {#Bisection for UB
      delta=(delta1+delta2)/2
      z=score.z(y,n,delta,MN)
      if (z>c) delta2=delta else delta1=delta
    }
  }
  delta.u=delta1
    return(switch(type,
                  "lower"=delta.l,
                  "upper"=delta.u,
                  "two.sided"=cbind(delta.l,delta.u)
  ))
}

## Example:
y=c(12,4)
n=c(20,20)
score.int(y,n)