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/e20) { 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))t])+sum(pnull[tnull[id1,id2]==t])*0.5 pvalue2[b]=sum(pnull[tnull[id1,id2]t])+sum(pnull[tnull[id1,id2]==t]) pvalue2[b]=sum(pnull[tnull[id1,id2]=(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=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))t])+sum(pnull[tnull[id1,id2]==t])*0.5 pvalue2[b]=sum(pnull[tnull[id1,id2]t])+sum(pnull[tnull[id1,id2]==t]) pvalue2[b]=sum(pnull[tnull[id1,id2]=(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)) }