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.