Or more elegantly the function below where a and b are the parameters of the beta prior, xa and xb are the current number of events in group A and B respectively; na and nb are the current total number of subjects in group A and B respectively; Na and Nb are the final total number of subject in group A and B respectively; and alpha is the alpha level (from 0 to 1). Don't forget to empty your folder: "DATA/predictivepower" (that you create in advance of course) after each computation.
----------------------------------------------------------------------------- ##Predictive Probability Interim Analysis, John Cook, 2006 predictive.power=function (a,b,xa,xb,na,nb,Na,Nb,alpha){ setwd("C:/R-2.4.1/DATA/predictivepower") 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)} Nta=Na-na for(i in 0:Nta) {sink("pw.power1.txt",append=TRUE); v=prepower(a,b,xa,(na-xa),i,Nta-i) dput(v) sink() } Ntb=Nb-nb for(i in 0:Ntb) {sink("pw.power2.txt",append=TRUE); v=prepower(a,b,xb,(nb-xb),i,Ntb-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/Nb))) z } for(i in 0:Nta){sink("pw.zval.txt",append=TRUE) zval=(zth(xa,i,Na,xb,0:Ntb,Nb)) cat(zval,sep="\n") sink() } zz=scan(file="pw.zval.txt") for(i in 1:(Nta+1)){sink("pw.predict.txt",append=TRUE) pw=((x[i]*y[1:(Ntb+1)])) cat(pw,sep="\n") sink() } p=scan(file="pw.predict.txt") z=-qnorm(alpha) mdf <- cbind(zz, p) d1=subset(mdf,apply(mdf>(z),1,any)) d2=subset(mdf,apply(mdf<(-z),1,any)) s1=sum(d1[,2]) s2=sum(d2[,2]) s3=1-(s1+s2) print(c("p A>B"=s1,"p B>A"=s2, "p A=B"=s3)) } -- View this message in context: http://www.nabble.com/Conditional-power%2C-predictive-power-tf3603396.html#a10094068 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.