rm(list=ls())

## TO run "meta.RD" for the risk difference between two proportions, one needs to install "binom" package first

library(binom)



## ==================================================================================================== ##
## ==================================================================================================== ##
## === meta.RD is the function for Obtaining the Combined Confidence Interval and P-value           === ##
## ===                              for the Risk Difference based on exact study-specific intervals === ##
## === The parameter of interest is the difference between two proportions                          === ##
## ==================================================================================================== ##
## ==================================================================================================== ##


## ==================================================================================================== ##
## ==================================================================================================== ##
## === meta.pRD is the function for Obtaining the Combined Confidence Interval and P-value          === ##
## ===                              for the Rate Difference based on exact study-specific intervals === ##
## === The parameter of interest is the difference between two poisson rates                        === ##
## ==================================================================================================== ##
## ==================================================================================================== ##




poisson.confint=function(x, prob)
                {alpha=1-prob
                 n=length(x)
                 lower=upper=rep(0, n)
                 for(i in 1:n)
                    {x0=x[i]
                     if(x0==0)
                        {f=function(lambda){ppois(0, lambda)-alpha}
                         fit=uniroot(f, c(0, 1000))
                         upper[i]=fit$root  
                         }
                     if(x0>0)
                        {f1=function(lambda){ppois(x0, lambda)-dpois(x0, lambda)*0.5-alpha/2}
                         f2=function(lambda){1-ppois(x0, lambda)+dpois(x0, lambda)*0.5-alpha/2}
                         fit1=uniroot(f1, c(0, 1000))
                         fit2=uniroot(f2, c(0, 1000))
                         upper[i]=fit1$root
                         lower[i]=fit2$root
                         }
                      }
                 return(list(lower=lower, upper=upper))
                 }



diff.pexact=function(x1,x2,e1,e2,delta.grd, n.grd=15, midp=T)
      {
        fit1=poisson.confint(x1,0.9995); 
        fit2=poisson.confint(x2,0.9995);

        rev=F
        if(fit2$upper/e2-fit2$lower/e2<fit1$upper/e1-fit1$lower/e1)
          {fit1=fit2
           e3=e2
           x3=x2

           e2=e1
           x2=x1                
           e1=e3
           x1=x3

           rev=T}  

        l1=fit1$lower 
        u1=fit1$upper
        N1=qpois(0.999999, u1)

        lambda1.grd=seq(l1, u1, length=n.grd)/e1; 
        pnull1.tot=matrix(0, N1+1, n.grd)
        for(b in 1:n.grd)
           {lambda1=lambda1.grd[b]
            pnull1.tot[,b]=dpois((0:N1),lambda1*e1)}

        fit2=poisson.confint(x2,0.9995); 
        u2=max(fit2$upper, e2*(u1/e1+max(delta.grd)))
        N2=qpois(0.999999, u2)

      
        dfnull=matrix(0, N1+1,N2+1)
        sdnull=matrix(0, N1+1,N2+1)
        for(i in 0:N1)
           {lambda1=(i+0.5)/(e1+1)
            lambda2=(0:N2+0.5)/(e2+1)
            dfnull[i+1, ]=(0:N2)/e2-i/e1
            sdnull[i+1, ]=sqrt(lambda1/(e1+1)+lambda2/(e2+1))
           }

        pv1=pv2=numeric(0)
        for(theta in delta.grd)
           {lambda1=(x1+0.5)/(e1+1)
            lambda2=(x2+0.5)/(e2+1)
            t=(x2/e2-x1/e1-theta)/sqrt(lambda1/(e1+1)+lambda2/(e2+1))
            tnull=(dfnull-theta)/sdnull


            pvalue1=pvalue2=rep(0,n.grd); error=1e-6
            
            if(max(lambda1.grd)+theta<=0)
                    { pv1=c(pv1, 0)
                      pv2=c(pv2, 1) }
            if(max(lambda1.grd)+theta>0)
            {
            for(b in 1:n.grd)
               {lambda1=lambda1.grd[b]
                lambda2=lambda1+theta
                if(lambda2>=0)
                  {pnull1=pnull1.tot[,b]
                   pnull2=dpois((0:N2), lambda2*e2)
                   n1.adj=N1+1-max(c(1,(1:(N1+1))[cumsum(sort(pnull1))<error]))
                   n2.adj=N2+1-max(c(1,(1:(N2+1))[cumsum(sort(pnull2))<error]))
                   id1=order(pnull1)
                   id2=order(pnull2)
                   id1=(id1[(N1+1):1])[1:n1.adj]
                   id2=(id2[(N2+1):1])[1:n2.adj]
                   pnull=pnull1[id1]%*%t(pnull2[id2])
                   if(midp==T)
                     {pvalue1[b]=sum(pnull[tnull[id1,id2]>t])+sum(pnull[tnull[id1,id2]==t])*0.5
                      pvalue2[b]=sum(pnull[tnull[id1,id2]<t])+sum(pnull[tnull[id1,id2]==t])*0.5
                     }else{
                          pvalue1[b]=sum(pnull[tnull[id1,id2]>t])+sum(pnull[tnull[id1,id2]==t])
                          pvalue2[b]=sum(pnull[tnull[id1,id2]<t])+sum(pnull[tnull[id1,id2]==t])
                          }
                    }
    
                
             }
             pv1=c(pv1, max(pvalue1)+(1-0.9995))
             pv2=c(pv2, max(pvalue2)+(1-0.9995))
           }
         }
        
        if(rev==T)
          {pv3=pv2
           pv2=rev(pv1)
           pv1=rev(pv3)
           }
          
        return(list(pv1=pv1, pv2=pv2))
      }


meta.pRD=function(data.mi, BB=1000, alpha=seq(0.1,0.95,length=20), B=20000, cov.prob=0.95, MH.imputation=F)
  {
    n=length(data.mi[,1])
    e1=data.mi[,3];             e2=data.mi[,4]
    lambda1=data.mi[,1]/data.mi[,3]; lambda2=data.mi[,2]/data.mi[,4]
 
    if(MH.imputation==T)
      {id=(1:n)[lambda1*lambda2==0]
       lambda1[id]=(data.mi[id,1]+0.5)/(data.mi[id,3]+1);  lambda2[id]=(data.mi[id,2]+0.5)/(data.mi[id,4]+1)
       e1[id]=data.mi[id,3]+1;                        e2[id]=data.mi[id,4]+1
       }

    deltalambda=lambda2-lambda1
    varlambda=lambda2/data.mi[,4]+lambda1/data.mi[,3]
    weight=(e1*e2/(e1+e2))/sum(e1*e2/(e1+e2))
    mu.MH=sum(deltalambda*weight);  sd.MH=sqrt(sum(weight^2*varlambda))
    ci.MH=c(mu.MH-qnorm((1+cov.prob)/2)*sd.MH, mu.MH+qnorm((1+cov.prob)/2)*sd.MH)
    p.MH=1-pchisq(mu.MH^2/sd.MH^2,1)
 
    d0=max(abs(ci.MH))
    delta.grd=seq(-d0*5, d0*5,length=BB-1); delta.grd=sort(c(0,delta.grd))

    

    pv1.pool=pv2.pool=numeric(0)
    for(kk in 1:n)
      { xx1=data.mi[kk,1]
        xx2=data.mi[kk,2] 
        ee1=data.mi[kk,3] 
        ee2=data.mi[kk,4]
        fit=diff.pexact(xx1, xx2, ee1, ee2, delta.grd)
        pv1.pool=rbind(pv1.pool, fit$pv1); pv2.pool=rbind(pv2.pool, fit$pv2)
        print(kk)
      }


    for(i in 1:n)
      { for(j in 1:BB)
          { pv1.pool[i,(BB-j+1)]=max(pv1.pool[i,1:(BB-j+1)]);pv2.pool[i,j]=max(pv2.pool[i,j:BB])
          }
      }
  

    K=length(alpha); sigma0=1/(data.mi[,3]+data.mi[,4])
    weight=matrix(0,K,4); weight[,1]=1; weight[,2]=1/(alpha*(1-alpha)); weight[,3]=1/(1-alpha)

    set.seed(100)
    t0=matrix(0,B,3); c=matrix(0,K,n)
    for(b in 1:B)
      { y=runif(n)
        for(k in 1:K){c[k,]=(y<alpha[k])-alpha[k]}
        t0[b,1]=sum(t(c*weight[,1])/sigma0); t0[b,2]=sum(t(c*weight[,2])/sigma0)
        t0[b,3]=sum(t(c*weight[,3])/sigma0)
      }

    t0=t0/n; alpha0=(1+cov.prob)/2; cut0=rep(0,3)
    for(b in 1:3){cut0[b]=quantile(t0[,b], 1-alpha0)}


    t1=t2=matrix(0,BB,3); c1=c2=c

    for(b in 1:BB)
      { for(k in 1:K)    
          { c1[k,]=(pv1.pool[,b]>=(1-alpha[k]))-alpha[k]
            c2[k,]=(pv2.pool[,b]>=(1-alpha[k]))-alpha[k]
          }   
        t1[b,1]=sum(t(c1*weight[,1])/sigma0);   t2[b,1]=sum(t(c2*weight[,1])/sigma0)
        t1[b,2]=sum(t(c1*weight[,2])/sigma0);   t2[b,2]=sum(t(c2*weight[,2])/sigma0)
        t1[b,3]=sum(t(c1*weight[,3])/sigma0);   t2[b,3]=sum(t(c2*weight[,3])/sigma0) 
      }

    t1=t1/n; t2=t2/n

    ci.cons=c(min(delta.grd[t1[,1]>=cut0[1]]),max(delta.grd[t2[,1]>=cut0[1]]))
    ci.iv=c(min(delta.grd[t1[,2]>=cut0[2]]),max(delta.grd[t2[,2]>=cut0[2]]))
    ci.fisher=c(min(delta.grd[t1[,3]>=cut0[3]]),max(delta.grd[t2[,3]>=cut0[3]]))
    ci.range=c(min(delta.grd), max(delta.grd))

    est.cons=delta.grd[abs(t2[,1]-t1[,1])==min(abs(t2[,1]-t1[,1]))][1]    
    est.iv=delta.grd[abs(t2[,2]-t1[,2])==min(abs(t2[,2]-t1[,2]))][1]    
    est.fisher=delta.grd[abs(t2[,3]-t1[,3])==min(abs(t2[,3]-t1[,3]))][1]    
    est.MH=mu.MH
    est.range=99999

    n0=max((1:BB)[delta.grd==0])
    c1=t1[n0,]; c2=t2[n0,]

    p.cons=min(1, 2*min(c(1-mean(t0[,1]>=c1[1]), 1-mean(t0[,1]>=c2[1]))))
    p.iv=min(1, 2*min(c(1-mean(t0[,2]>=c1[2]), 1-mean(t0[,2]>=c2[2]))))
    p.fisher=min(1, 2*min(c(1-mean(t0[,3]>=c1[3]), 1-mean(t0[,3]>=c2[3]))))

    pvalue=c(p.cons, p.iv, p.fisher, p.MH)
    names(pvalue)=c("constant", "inverse-variance", "fisher", "asymptotical-MH")
    ci=cbind(ci.cons, ci.iv, ci.fisher,ci.MH, ci.range)
    ci=rbind(c(est.cons, est.iv, est.fisher, est.MH,  est.range), ci)
    rownames(ci)=c("point est", "lower end of CI", "upper end of CI")
    colnames(ci)=c("constant", "inverse-variance", "fisher", "asymptotical-MH", "search range")
  
###########################################################################################
    id=(1:n)[lambda1*lambda2==0]
    lambda1[id]=(data.mi[id,1]+0.5)/(data.mi[id,3]+1);  lambda2[id]=(data.mi[id,2]+0.5)/(data.mi[id,4]+1)
    deltalambda=lambda2-lambda1
    varlambda=lambda2/data.mi[,4]+lambda1/data.mi[,3]
    sigmadelta=sqrt(varlambda)
    study.ci=matrix(0, n, 4)
    colnames(study.ci)=c("point est", "lower end of CI", "upper end of CI", "search limit")
    for(kk in 1:n) 
       {deltastudy.grd=seq(deltalambda[kk]-10*sigmadelta[kk], deltalambda[kk]+10*sigmadelta[kk],length=BB)
        xx1=data.mi[kk,1]
        xx2=data.mi[kk,2] 
        ee1=data.mi[kk,3] 
        ee2=data.mi[kk,4]
        fit=diff.pexact(xx1, xx2, ee1, ee2, deltastudy.grd)
        pv1=fit$pv1
        pv2=fit$pv2
        for(j in 1:BB)
           {pv1[BB-j+1]=max(pv1[1:(BB-j+1)])
            pv2[j]=max(pv2[j:BB])
            }
        study.ci[kk,2:3]=range(deltastudy.grd[pv1>(1-cov.prob)/2 & pv2>(1-cov.prob)/2])
        study.ci[kk,1]=deltastudy.grd[abs(pv2-pv1)==min(abs(pv2-pv1))][1]
        study.ci[kk,4]="OK"
        if(min(deltastudy.grd)==study.ci[,2] || max(deltastudy.grd)==study.ci[,3])
        study.ci[kk,4]="Reach Limit"
        }
################################################################################################

    plot(range(delta.grd), c(0, n+5), xlab="rate difference", ylab="study", type="n")
    for(i in 1:n)
       {lines( study.ci[i,2:3], rep(n+5-i+1, 2))}
    lines(ci.cons, rep(4, 2), lwd=2, col=2)
    lines(ci.iv, rep(3, 2), lwd=2, col=3)
    lines(ci.fisher, rep(2, 2), lwd=2, col=4)
    lines(ci.MH, rep(1, 2), lwd=2, col=5)
    lines(c(0, 0), c(0, n+6))

    text(min(delta.grd), 4, "const", col=2)
    text(min(delta.grd), 3, "iv", col=3)
    text(min(delta.grd), 2, "fisher", col=4)
    text(min(delta.grd), 1, "MH", col=5)


   

    return(list(ci.fixed=ci, study.ci=study.ci, accuracy=(max(delta.grd)-min(delta.grd))/BB))
 }





##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################


binom.confint=function(x, n, alpha, method="exact")
      {error.upper=1-(1-alpha)/2
       error.lower=(1-alpha)/2
       eq1=function(p){pbinom(x, n, p)-error.lower}
       fit1=uniroot(eq1, c(0, 1))
       eq2=function(p){pbinom(x, n, p)-error.upper}
       fit2=uniroot(eq2, c(0, 1))
       lower=fit2$root
       upper=fit1$root
       return(list(lower=lower, upper=upper))
       }
       





diff.exact=function(x1,x2,n1,n2,delta.grd, n.grd=15, midp=T)
      {
        fit1=binom.confint(x1,n1,0.9995,method="exact"); 
        fit2=binom.confint(x2,n2,0.9995,method="exact");

        rev=F
        if(fit2$upper-fit2$lower<fit1$upper-fit1$lower)
          {fit1=fit2
           n3=n2
           x3=x2
           n2=n1
           x2=x1                
           n1=n3
           x1=x3

           rev=T}  


        l=fit1$lower; 
        u=fit1$upper


        p1.grd=seq(l, u, length=n.grd); 
        pnull1.tot=matrix(0, n1+1, n.grd)
        for(b in 1:n.grd)
           {p1=p1.grd[b]
            pnull1.tot[,b]=dbinom((0:n1),n1,p1)}
        dfnull=matrix(0, n1+1,n2+1)
        sdnull=matrix(0, n1+1,n2+1)
        for(i in 0:n1)
           {p1=(i+0.5)/(n1+1)
            p2=(0:n2+0.5)/(n2+1)
            dfnull[i+1, ]=(0:n2)/n2-i/n1
            sdnull[i+1, ]=sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
           }
        pv1=pv2=numeric(0)
        for(theta in delta.grd)
           {p1=(x1+0.5)/(n1+1)
            p2=(x2+0.5)/(n2+1)
            t=(x2/n2-x1/n1-theta)/sqrt(p1*(1-p1)/n1+p2*(1-p2)/n2)
            tnull=(dfnull-theta)/sdnull
            pvalue1=pvalue2=rep(0,n.grd); error=1e-6
            if(min(p1.grd)+theta>=1)
                    { pv1=c(pv1, 1)
                      pv2=c(pv2, 0) }
            if(max(p1.grd)+theta<=0)
                    { pv1=c(pv1, 0)
                      pv2=c(pv2, 1) }
            if(min(p1.grd)+theta<1 && max(p1.grd)+theta>0)
            {
            for(b in 1:n.grd)
               {p1=p1.grd[b]
                p2=p1+theta
                if(p2>=0 && p2<=1)
                  {pnull1=pnull1.tot[,b]
                   pnull2=dbinom((0:n2), n2, p2)
                   n1.adj=n1+1-max(c(1,(1:(n1+1))[cumsum(sort(pnull1))<error]))
                   n2.adj=n2+1-max(c(1,(1:(n2+1))[cumsum(sort(pnull2))<error]))
                   id1=order(pnull1)
                   id2=order(pnull2)
                   id1=(id1[(n1+1):1])[1:n1.adj]
                   id2=(id2[(n2+1):1])[1:n2.adj]
                   pnull=pnull1[id1]%*%t(pnull2[id2])
                   if(midp==T)
                      { pvalue1[b]=sum(pnull[tnull[id1,id2]>t])+sum(pnull[tnull[id1,id2]==t])*0.5
                        pvalue2[b]=sum(pnull[tnull[id1,id2]<t])+sum(pnull[tnull[id1,id2]==t])*0.5
                      }else{
                        pvalue1[b]=sum(pnull[tnull[id1,id2]>t])+sum(pnull[tnull[id1,id2]==t])
                        pvalue2[b]=sum(pnull[tnull[id1,id2]<t])+sum(pnull[tnull[id1,id2]==t])
                      }
                  }
    
              }
            pv1=c(pv1, max(pvalue1)+(1-0.9995)); pv2=c(pv2, max(pvalue2)+(1-0.9995))
           }
          }
        
        if(rev==T)
          {pv3=pv2
           pv2=rev(pv1)
           pv1=rev(pv3)
           }
          
        return(list(pv1=pv1, pv2=pv2))
      }



meta.RD=function(data.mi, BB=1000, alpha=seq(0.1,0.95,length=20), B=20000, cov.prob=0.95, MH.imputation=F)
  {
    n=length(data.mi[,1])
    n1=data.mi[,3];             n2=data.mi[,4]
    p1=data.mi[,1]/data.mi[,3]; p2=data.mi[,2]/data.mi[,4]
 
    if(MH.imputation==T)
      {id=(1:n)[p1*p2==0]
       p1[id]=(data.mi[id,1]+0.5)/(data.mi[id,3]+1);  p2[id]=(data.mi[id,2]+0.5)/(data.mi[id,4]+1)
       n1[id]=data.mi[id,3]+1;                        n2[id]=data.mi[id,4]+1
       }

    deltap=p2-p1
    varp=p1*(1-p1)/n1+p2*(1-p2)/n2
    weight=(n1*n2/(n1+n2))/sum(n1*n2/(n1+n2))
    mu.MH=sum(deltap*weight);  sd.MH=sqrt(sum(weight^2*varp))
    ci.MH=c(mu.MH-qnorm((1+cov.prob)/2)*sd.MH, mu.MH+qnorm((1+cov.prob)/2)*sd.MH)
    p.MH=1-pchisq(mu.MH^2/sd.MH^2,1)
 
    d0=max(abs(ci.MH))
    delta.grd=seq(max(-1, -d0*5), min(1, d0*5),length=BB-1); delta.grd=sort(c(0,delta.grd))

    

    pv1.pool=pv2.pool=numeric(0)
    for(kk in 1:n)
      { x1=data.mi[kk,1]
        x2=data.mi[kk,2]
        n1=data.mi[kk,3] 
        n2=data.mi[kk,4]
        fit=diff.exact(x1,x2,n1,n2, delta.grd)
        pv1.pool=rbind(pv1.pool, fit$pv1); pv2.pool=rbind(pv2.pool, fit$pv2)
        print(kk)
      }


    for(i in 1:n)
      { for(j in 1:BB)
          { pv1.pool[i,(BB-j+1)]=max(pv1.pool[i,1:(BB-j+1)]);pv2.pool[i,j]=max(pv2.pool[i,j:BB])
          }
      }
  

    K=length(alpha); sigma0=1/(data.mi[,3]+data.mi[,4])
    weight=matrix(0,K,4); weight[,1]=1; weight[,2]=1/(alpha*(1-alpha)); weight[,3]=1/(1-alpha)


    set.seed(100)
    t0=matrix(0,B,3); c=matrix(0,K,n)
    for(b in 1:B)
      { y=runif(n)
        for(k in 1:K){c[k,]=(y<alpha[k])-alpha[k]}
        t0[b,1]=sum(t(c*weight[,1])/sigma0); t0[b,2]=sum(t(c*weight[,2])/sigma0)
        t0[b,3]=sum(t(c*weight[,3])/sigma0)
      }

    t0=t0/n; alpha0=(1+cov.prob)/2; cut0=rep(0,3)
    for(b in 1:3){cut0[b]=quantile(t0[,b], 1-alpha0)}

    

    t1=t2=matrix(0,BB,3); c1=c2=c

    for(b in 1:BB)
      { for(k in 1:K)    
          { c1[k,]=(pv1.pool[,b]>=(1-alpha[k]))-alpha[k]
            c2[k,]=(pv2.pool[,b]>=(1-alpha[k]))-alpha[k]
          }   
        t1[b,1]=sum(t(c1*weight[,1])/sigma0);   t2[b,1]=sum(t(c2*weight[,1])/sigma0)
        t1[b,2]=sum(t(c1*weight[,2])/sigma0);   t2[b,2]=sum(t(c2*weight[,2])/sigma0)
        t1[b,3]=sum(t(c1*weight[,3])/sigma0);   t2[b,3]=sum(t(c2*weight[,3])/sigma0) 
      }

    t1=t1/n; t2=t2/n

    ci.cons=c(min(delta.grd[t1[,1]>=cut0[1]]),max(delta.grd[t2[,1]>=cut0[1]]))
    ci.iv=c(min(delta.grd[t1[,2]>=cut0[2]]),max(delta.grd[t2[,2]>=cut0[2]]))
    ci.fisher=c(min(delta.grd[t1[,3]>=cut0[3]]),max(delta.grd[t2[,3]>=cut0[3]]))
    ci.range=c(min(delta.grd), max(delta.grd))

    est.cons=delta.grd[abs(t2[,1]-t1[,1])==min(abs(t2[,1]-t1[,1]))][1]    
    est.iv=delta.grd[abs(t2[,2]-t1[,2])==min(abs(t2[,2]-t1[,2]))][1]    
    est.fisher=delta.grd[abs(t2[,3]-t1[,3])==min(abs(t2[,3]-t1[,3]))][1]    
    est.range=99999
    est.MH=mu.MH



    n0=max((1:BB)[delta.grd==0])
    c1=t1[n0,]; c2=t2[n0,]

    p.cons=min(1, 2*min(c(1-mean(t0[,1]>=c1[1]), 1-mean(t0[,1]>=c2[1]))))
    p.iv=min(1, 2*min(c(1-mean(t0[,2]>=c1[2]), 1-mean(t0[,2]>=c2[2]))))
    p.fisher=min(1, 2*min(c(1-mean(t0[,3]>=c1[3]), 1-mean(t0[,3]>=c2[3]))))

    pvalue=c(p.cons, p.iv, p.fisher, p.MH)
    names(pvalue)=c("constant", "inverse-variance", "fisher", "asymptotical-MH")
    ci=cbind(ci.cons, ci.iv, ci.fisher,ci.MH, ci.range)
    ci=rbind(c(est.cons, est.iv, est.fisher, est.MH, est.range), ci)
    rownames(ci)=c("point est", "lower end of CI", "upper end of CI")
    colnames(ci)=c("constant", "inverse-variance", "fisher", "asymptotical-MH", "search range")
  
###################################################################################################  
    n1=data.mi[,3];             n2=data.mi[,4]
    p1=data.mi[,1]/data.mi[,3]; p2=data.mi[,2]/data.mi[,4]
 
    id=(1:n)[p1*p2==0]
    p1[id]=(data.mi[id,1]+0.5)/(data.mi[id,3]+1);  p2[id]=(data.mi[id,2]+0.5)/(data.mi[id,4]+1)
    n1[id]=data.mi[id,3]+1;                        n2[id]=data.mi[id,4]+1       
    deltap=p2-p1
    varp=p1*(1-p1)/n1+p2*(1-p2)/n2

    sigmadelta=sqrt(varp)
    study.ci=matrix(0, n, 4)
    colnames(study.ci)=c("point est", "lower end of CI", "upper end of CI", "search limit")
    for(kk in 1:n) 
       {deltastudy.grd=seq(deltap[kk]-10*sigmadelta[kk], deltap[kk]+10*sigmadelta[kk],length=BB)
        xx1=data.mi[kk,1]
        xx2=data.mi[kk,2] 
        ee1=data.mi[kk,3] 
        ee2=data.mi[kk,4]
        fit=diff.exact(xx1, xx2, ee1, ee2, deltastudy.grd)
        pv1=fit$pv1
        pv2=fit$pv2
        for(j in 1:BB)
           {pv1[BB-j+1]=max(pv1[1:(BB-j+1)])
            pv2[j]=max(pv2[j:BB])
            }
        study.ci[kk,2:3]=range(deltastudy.grd[pv1>(1-cov.prob)/2 & pv2>(1-cov.prob)/2])
        study.ci[kk,1]=deltastudy.grd[abs(pv2-pv1)==min(abs(pv2-pv1))][1]
        study.ci[kk,4]="OK"
        if(min(deltastudy.grd)==study.ci[,2] || max(deltastudy.grd)==study.ci[,3])
          study.ci[kk,4]="Reach Limit"
        }


    plot(range(delta.grd), c(0, n+5), xlab="risk difference", ylab="study", type="n")
    for(i in 1:n)
       {lines( study.ci[i,2:3], rep(n+5-i+1, 2))}
    lines(ci.cons, rep(4, 2), lwd=2, col=2)
    lines(ci.iv, rep(3, 2), lwd=2, col=3)
    lines(ci.fisher, rep(2, 2), lwd=2, col=4)
    lines(ci.MH, rep(1, 2), lwd=2, col=5)
    lines(c(0, 0), c(0, n+6))

    text(min(delta.grd), 4, "const", col=2)
    text(min(delta.grd), 3, "iv", col=3)
    text(min(delta.grd), 2, "fisher", col=4)
    text(min(delta.grd), 1, "MH", col=5)
  
 


    return(list(ci.fixed=ci, study.ci=study.ci, accuracy=(max(delta.grd)-min(delta.grd))/BB))
 }