# admeanRM.s : admean.s + REML for difference in means # d<-scan("c:/MetaAnalysis/Rcode/normal.txt",list(n1=0,m1=0,s1=0,n0=0,m0=0,s0=0)) k<-length(d$n1) # # parameters for figures xposi <- 60 # x-axis for titles yup <- 11 # upper limit of y-axis ylow <- -4 # lower limit of y-axis xup <- 100 # upper limit of x-axis xlow <- -100 # lower limit of x-axis # # par(mar=c(5,4,4,3)) ad<-d$m1-d$m0 s<-( (d$n1-1)*d$s1^2 + (d$n0-1)*d$s0^2 )/(d$n1+d$n0-2) se<-sqrt((1/d$n1+1/d$n0)*s) tn<-d$n1+d$n0 low<-ad-1.96*se upp<-ad+1.96*se # ---------- individual graph ------------------ id<-k:1 plot(ad,id,ylim=c(ylow, yup), xlim=c(xlow, xup),yaxt="n",pch=" ", ylab="Citation",xlab="Difference in means") title(main=" Mean difference model ") symbols(ad,id,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) lines(x,y,type="l") text(60 ,i,j) } # ----- fixed effects ---- w<-1/se/se sw<-sum(w) adm<-sum(ad*w)/sw adl<-adm-1.96*sqrt(1/sw) adu<- adm+1.96*sqrt(1/sw) q1<-sum(w*(ad-adm)^2) df1<-k-1 pval1<- 1-pchisq(q1,df1) # I2 <- (q1-df1)/q1 * 100 # q2<-adm^2*sw df2<-1 pval2<- 1-pchisq(q2,df2) # ----- random-effects ---- tau2<-(q1-(k-1))/(sw-sum(w*w)/sw) tau2<-max(0,tau2) wx<-1/(tau2+se*se) swx<-sum(wx) admx<-sum(ad*wx)/swx adxl<- admx-1.96*sqrt(1/swx) adxu<- admx+1.96*sqrt(1/swx) qx2<-admx^2*swx pvalx2<- 1-pchisq(qx2,df2) # -------- graph ---------------- x<-c(adl,adu) y<-c(-1, -1) lines(x,y,type="b") x<-c(adxl,adxu) y<-c(-2, -2) lines(x,y,type="b") abline(v=c(adm, admx), lty=2) text(xposi,-1, "Combined : fixed") text(xposi,-2, "Combined : random") # ------ output variables ------- out <- round( c(adm, adl, adu), 4) outq1<- round( c(q1, df1, pval1), 6) outq2<- round( c(q2, df2, pval2), 6) outR<- round( c(admx, adxl, adxu), 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) #======================================= wb2<-1/(tau+se*se)/(tau+se*se) sb<-sum(wb) sb2<-sum(wb2) dw<- -wb2 sdw<- sum(dw) admb<-sum(ad*wb)/sum(wb) h<-sb*sum(wb-wb2*(ad-admb )^2) - sb2 saa<-sum(dw - 2*wb*dw * (ad-admb)^2 ) dh<-sdw*sum(wb-wb2*(ad-admb )^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) adRM<-sum(ad*wg)/swg adRMl<- adRM-1.96*sqrt(1/swg) adRMu<- adRM+1.96*sqrt(1/swg) qx2RM<-adRM^2*swg pvalx2RM<- 1-pchisq(qx2RM,df2) x<-c(adRMl,adRMu) y<-c(-3, -3) lines(x,y,type="b") tau2 <- tau outRM<- round( c(adRM, adRMl, adRMu), 4) outq2RM<- round( c(qx2RM, df2, pvalx2RM), 6) text(xposi,-3, "Combined : REML")