# peto.s: Peto's method # dat<-scan("c:/splus/beta2.txt",list(a=0,n1=0,c=0,n0=0,d=0)) k<-length(dat$a) # # parameters for figures xposi <- 4 # x-axis for titles yup <- 20 # upper limit of y-axis ylow <- -2 # lower limit of y-axis # # par(mar=c(5,4,4,3)) 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 m1<-ai+ci m0<-bi+di o<-ai e<-n1*m1/tn v<-n1*n0*m1*m0/tn/tn/(tn-1) lgor<-(o-e)/v se<-1/sqrt(v) w<-1/se/se low<-exp(lgor-1.96*se) upp<-exp(lgor+1.96*se) # ---------- individual graph ------------------ id<-k:1 plot(exp(lgor),id,ylim=c(ylow, yup), log="x", xlim=c(0.1,10),yaxt="n",pch=" ", ylab="Citation",xlab="Odds ratio") title(main=" Peto's method ") symbols(exp(lgor),id,squares=sqrt(tn), add=T,inches=0.20) 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(xposi,i,j) } # ----- fixed effects ---- sw<-sum(w) peto <- exp( sum( lgor*w )/sw ) petol<- exp( log(peto)-1.96*sqrt(1/sw) ) petou<- exp( log(peto)+1.96*sqrt(1/sw) ) q1<-sum( w*(lgor-log(peto))^2 ) df1<-k-1 pval1<- 1-pchisq(q1,df1) q2<-(abs( sum(o-e) ) -0.5)^2/sum(v) df2<-1 pval2<- 1-pchisq(q2,df2) # -------- graph ---------------- x<-c(petol,petou) y<-c(-1, -1) lines(x,y,type="b",lty=1) abline(v=c(peto),lty=2) abline(v=1) text(xposi,-1, "Combined : fixed") # ------ output variables ------- out <- round( c(peto, petol, petou), 4) outq1<- round( c(q1, df1, pval1), 6) outq2<- round( c(q2, df2, pval2), 6)