##################################### ### 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)