#####################################
### Computes confidence Interval  ###
### for the difference of         ###
### two proportions based on      ###
### inverting the score test      ###
#####################################
#based on relationship/identity in Andres et al. (2011) Stat. Meth. Med. Research
#if function value <0 than score statistic is not significant at 1-conflev level

score.sig <- function(delta,p,n,comb,z2) {
  C= z2/(sum(comb*p) - delta)
  b=1-2*p
  R=sqrt(n^2+2*n*comb*b*C+comb^2*C^2)
  sum(n)+(sum(comb)-2*delta)*C-sum(R)
}

score.int <-function(y,n,conflev=0.95, alternative="two-sided", comb=c(1,-1)) {
  delta.max=sum(comb[comb>=0])
  delta.min=sum(comb[comb<0])
  p=y/n
  point.est=sum(comb*p)
  eps=10^(-5)
  if (alternative=="two-sided") {
      z=qnorm(1-(1-conflev)/2)
      if(point.est==delta.min) delta.lb=delta.min   #to handle y=0 or y=n
      else delta.lb=uniroot(score.sig,interval=c(delta.min,point.est-eps),p=p,n=n,comb=comb,z2=z^2)$root
      if(point.est==delta.max) delta.ub=delta.max   #to handle y=0 or y=n
      else delta.ub=uniroot(score.sig,interval=c(point.est+eps,delta.max),p=p,n=n,comb=comb,z2=z^2)$root
      return(c(delta.lb,delta.ub))
      }
  else {
      z=qnorm(1-(1-conflev))
      if (alternative=="lower") {
        if(point.est==delta.min) delta.lb=delta.min   #to handle y=0 or y=n
        else delta.lb=uniroot(score.sig,interval=c(delta.min,point.est-eps),p=p,n=n,comb=comb,z2=z^2)$root
      return(delta.lb)
      }
      else {
        if(point.est==delta.max) delta.ub=delta.max   #to handle y=0 or y=n
        else delta.ub=uniroot(score.sig,interval=c(point.est+eps,delta.max),p=p,n=n,comb=comb,z2=z^2)$root
        return(delta.ub)
        }
      }
}

#Example
y=c(2,48)
n=c(50,50)
score.int(y,n)