# varrd.s : Variance- based method for risk difference # plus cumulative meta analysis shown by "circles" # dat<-scan("c:/splus/DMincidence.txt",list(a=0,n1=0,c=0,n0=0,d=0)) k<-length(dat$a) # # parameters for figures xposi <- 0.2 # x-axis for titles yup <- 10 # upper limit of y-axis ylow <- -5 # lower limit of y-axis xup <- 0.5 # upper limit of x-axis xlow <- -0.5 # lower limit of x-axis label<-rep("",k) label[1]<-"(1) 1997" label[2]<-"(2) 1999" label[3]<-"(3) 2001" label[4]<-"(4) 2002" label[5]<-"(5) 2003" label[6]<-"(6) 2005" label[7]<-"(7) 2006" # par(mar=c(5,4,4,3)) dd<-dat ai<-dat$a bi<-dat$n1-dat$a ci<-dat$c di<-dat$n0-dat$c tn<-dat$n1+dat$n0 n1<-dat$n1 n0<-dat$n0 rd<- (ai/n1) - (ci/n0) se<-sqrt(ai*bi/n1^3 + ci*di/n0^3) low<-rd-1.96*se upp<-rd+1.96*se # ---------- individual graph ------------------ id<-k:1 ind<-1:k res<-0.3 plot(rd,id,ylim=c(ylow, yup), xlim=c(xlow, xup),yaxt="n",pch=" ", ylab="Citation",xlab="Risk difference") title(main=" Variance-based method ( circles are cumulative MA )") symbols(rd,id+res,squares=sqrt(tn), add=T,inches=0.20) abline(v=0) for (i in 1:k){ j<-k-i+1 x<-c(low[i],upp[i]) y<-c(j,j)+0.3 lines(x,y,type="l") text(xposi, i,label[j]) } # ----- cumulative graph ---- hh<-1:k xh<-hh xhl<-hh xhu<-hh for (ir in 1:k){ dat$a <-dd$a[ind<=ir] dat$n1<-dd$n1[ind<=ir] dat$c <-dd$c[ind<=ir] dat$n0<-dd$n0[ind<=ir] ai<-dat$a bi<-dat$n1-dat$a ci<-dat$c di<-dat$n0-dat$c tn<-dat$n1+dat$n0 n1<-dat$n1 n0<-dat$n0 rd<- (ai/n1) - (ci/n0) se<-sqrt(ai*bi/n1^3 + ci*di/n0^3) hh[ir]<-sum(tn) # ----- fixed effects ---- w<-1/se/se sw<-sum(w) varrd<- sum(rd*w)/sw varrdl<- varrd -1.96*sqrt(1/sw) varrdu<- varrd +1.96*sqrt(1/sw) q1<-sum( w*(rd-varrd)^2 ) df1<-k-1 pval1<- 1-pchisq(q1,df1) q2<-varrd^2*sw df2<-1 pval2<- 1-pchisq(q2,df2) # ----- random-effects ---- if (ir == 1){ varrdd<- rd[1] varrddl<-low[1] varrddu<-upp[1] } if (ir > 1){ tau2<-(q1-(k-1))/(sw-sum(w*w)/sw) tau2<-max(0,tau2) wx<-1/(tau2+se*se) swx<-sum(wx) varrdd<-sum(rd*wx)/swx varrddl<- varrdd-1.96*sqrt(1/swx) varrddu<- varrdd+1.96*sqrt(1/swx) qx2<-varrdd^2*swx pvalx2<- 1-pchisq(qx2,df2) } # -- graph of cumulative meta ---- j<- k-ir+1 x<-c(varrddl, varrddu) xh[ir]<- varrdd xhl[ir]<-varrddl xhu[ir]<-varrddu y<-c(j,j) lines(x,y,type="l",lty=3) } # -------- graph ---------------- symbols(xh,id,circles=sqrt(hh), add=T,inches=0.10) x<-c(varrdl,varrdu) y<-c(-1,-1) lines(x,y,type="b") x<-c(varrddl,varrddu) y<-c(-2,-2) lines(x,y,type="b") abline(v=c(varrd,varrdd), lty=2) text(xposi,-1, "Combined : fixed") text(xposi,-2, "Combined : random") # ------ output variables ------- out <- round( c(varrd, varrdl, varrdu), 4) outq1<- round( c(q1, df1, pval1), 6) outq2<- round( c(q2, df2, pval2), 6) outR<- round( c(varrdd, varrddl, varrddu), 4) outq2R<- round( c(qx2, df2, pvalx2), 6) # ----- REML method ---- intau<-tau2 tau<-intau # nrep<-500 newt<-1:nrep for (i in 1:nrep){ wb<-1/(tau+se*se) rdmb<-sum(rd*wb)/sum(wb) #======================================= wb2<-1/(tau+se*se)/(tau+se*se) sb<-sum(wb) sb2<-sum(wb2) dw<- -wb2 sdw<- sum(dw) h<-sb*sum(wb-wb2*(rd- rdmb )^2) - sb2 saa<-sum(dw - 2*wb*dw * (rd - rdmb )^2 ) dh<-sdw*sum(wb-wb2*(rd - rdmb )^2) + sb*saa - sum(2*wb*dw) #======================================= newt[i]<-tau - h/dh rel<-abs((newt[i]-tau)/tau) if (rel < 0.0000001) break tau<-newt[i] } if (i == 500) tau<-NaN # wg<-1/(tau+se*se) swg<-sum(wg) rdRM<- sum(rd*wg)/swg rdRMl<- rdRM -1.96*sqrt(1/swg) rdRMu<- rdRM +1.96*sqrt(1/swg) qx2RM<- rdRM^2*swg pvalx2RM<- 1-pchisq(qx2RM,df2) x<-c(rdRMl,rdRMu) y<-c(-3, -3) lines(x,y,type="b") text(xposi,-3, "Combined : REML") tau2<-tau outRM<- round( c(rdRM, rdRMl, rdRMu), 4) outq2RM<- round( c(qx2RM, df2, pvalx2RM), 6) #====================================== # Noma(2011)'s Bartlett corrected LR 95% CI # we need six parameters (y, var, muR, tauR, clow, cup) # y<-rd # data (effect size, e.g., # difference in means, log OR ) var<-se*se # variance muR <- varrdd # pooled effect size (random) tauR <- intau # estimated tau^2 (random) clow <- varrddl # lower limit of 95% CI (random) cup <- varrddu # upper limit of 95% CI (random) # outBC<-noma(y, var, muR, tauR, clow, cup) # x<-c(outBC[2], outBC[3]) y<-c(-4, -4) lines(x,y,type="b") text(xposi, -4, "Combined :random ") text(xposi, -4.6, "(Bartlett corrected LR)") outBLR<-round( outBC, 4)