score.sig <- function(delta,p,n,comb,z2) {
#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
  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)

## for just getting the score statistic
score.z <- function(delta, y, n) {#gives MLEs for p1,p2 under restriction p1-p2=delta and returns score statistic z (Agresti 2002, p77)
        N=sum(n)
        C=sum(y)
        L3=N
        L2=(n[1]+2*n[2])*delta-N-C
        L1=(n[2]*delta-N-2*y[2])*delta+C
        L0=y[2]*delta*(1-delta)
        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))
        a=(pi+acos(c/b^3))/3
        p2=2*b*cos(a)-L2/(3*L3)
        p1=p2+delta
        z=(y[1]/n[1]-y[2]/n[2]-delta)/sqrt(p1*(1-p1)/n[1] + p2*(1-p2)/n[2]) #Mee
        #z=(y[1]/n[1]-y[2]/n[2]-delta)/sqrt((p1*(1-p1)/n[1] + p2*(1-p2)/n[2])*N/(N-1)) #Mietinnen and Nurminen
        return(z)
}