Ok I am replying to my own message! I wrote a "function", it works well but it's a bit twisted because you will have to edit the last file in excel or other. This is to analyze the bayesian predictive power in an analysis where treatment x is compared to treatment y. Example: Total final subject number is 100 (randomization 1:1, so 50 per group). We do an interim analysis when there are 25 subjects in group x and 25 subjects in group y. Number of successes in x=10; number of successes in y=16. Prior=beta(0.6,0.4). ----------------------------------------------------------------------------------------------------------------
##Predictive Probability Interim Analysis, John Cook, 2006 prepower=function(a,b,s0,f0,s1,f1){ factorial(s1+f1)/ (factorial(s1)*factorial(f1))* beta(s0+s1+a,f0+f1+b)/beta(s0+a,f0+b)} for(i in 0:25) {sink("pw.power1.txt",append=TRUE); v=prepower(0.6,0.4,10,15,i,25-i) dput(v) sink() } for(i in 0:25) {sink("pw.power2.txt",append=TRUE); v=prepower(0.6,0.4,16,9,i,25-i) dput(v) sink() } x=scan(file="pw.power1.txt") y=scan(file="pw.power2.txt") zth=function(xa,ya,na,xb,yb,nb){ tha=(xa+ya)/na thb=(xb+yb)/nb th=(xa+xb+ya+yb)/(na+nb) z=(tha-thb)/sqrt(th*(1-th)*((1/na)+(1/na))) z } for(i in 0:25){sink("pw.zval.txt",append=TRUE) zval=(zth(10,i,50,16,0:25,50)) cat(zval,sep="\n") sink() } z=scan(file="pw.zval.txt") for(i in 1:26){sink("pw.predict.txt",append=TRUE) pw=((x[i]*y[1:26])) cat(pw,sep="\n") sink() } p=scan(file="pw.predict.txt") mdf <- data.frame(z, p) pw.bind=mdf[order(z,p),] write.table(pw.bind,file="pw.bind.txt",sep="\t") ##edit pw.bind file in excel for z>1.96, z<-1.96 or in between (or any other value); save as same file predictive=read.table(file="pw.bind.txt",header=TRUE) sum(predictive$p) ## This gives the bayesian predictive power -- View this message in context: http://www.nabble.com/Conditional-power%2C-predictive-power-tf3603396.html#a10079648 Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.