score.z <- function(y,n,delta.null, MN=TRUE) { 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)) a=(3.14159265358979+acos(c/b^3))/3 p2=2*b*cos(a)-L2/(3*L3) p1=p2+delta.null if (!MN) z=(y[1]/n[1]-y[2]/n[2]-delta.null)/sqrt(p1*(1-p1)/n[1] + p2*(1-p2)/n[2]) #Mee else z=(y[1]/n[1]-y[2]/n[2]-delta.null)/sqrt((p1*(1-p1)/n[1] + p2*(1-p2)/n[2])*N/(N-1)) #Mietinnen and Nurminen return(z^2) } #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^(-6) ) {#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^(-6) ) {#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) )) }