I am trying to calculate the Bayesian predictive probability P(P(Diff>0)>0.95), for two sample means. I assume the following
(1) Standard deviation = 1.44 for both groups
(2) mean 1 = 2.06, mean 2 = 1.34 so we expect to show a mean difference of 0.72
(3) In group 1, n = 80; In group 2, n = 40
(4) we try to borrow information for n = 30 from both group, so those 30 were treated as prior information, and data should be collected for n = 50 from group 1, and n = 10 from group 2.
My OpenBUGS code are
model{
# likelihoods
for (i in 1:n){
x[i] ~ dnorm(mua,taua)
y[i] ~ dnorm(mup,taup)
diff[i] <- x[i] - y[i]
difffl[i] <- step(diff[i])
}
# normal priors
mua ~ dnorm(mu1, inv.sigma)
mup ~ dnorm(mu2, inv.sigma)
inv.sigma <- n3/sigma.squared
taua <- n1/sigma.squared
taup <- n2/sigma.squared
# calculate the probability P(diff > delta)
# delta will be 0 here; # eta will be 0.95 here
# calculate the predictive probability P(P(diff >delta)>0.95)
sdx <- sd(x[])
sdy <- sd(y[])
mdiff <- mean(diff[])
mdifffl <- mean(difffl[])
powerfl <- step(mdifffl - 0.95)
}
list(mu1 = 2.06, mu2 = 1.34, n = 1000, n1 = 50, n2 = 10, n3 = 30,
sigma.squared = 2.0736)
And my results are
mean sd MC_error val2.5pc median val97.5pc start sample
mdiff 0.7152 0.3724 0.003383 -0.01549 0.7143 1.46 1 10000
mdifffl 0.8749 0.1376 0.001318 0.489 0.924 0.999 1 10000
powerfl 0.3932 0.4885 0.005197 0.0 0.0 1.0 1 10000
sdx 0.2036 0.004437 4.722E-5 0.195 0.2036 0.2124 1 10000
sdy 0.4552 0.01011 1.079E-4 0.4357 0.4551 0.4752 1 10000
The results don't look correct, as the standard deviation of mean in group 1 (sdx) should be 1.44/sqrt(80) = 0.161. current value 0.2036 in fact is 1.44/sqrt(50). The power I obtained is also too small (powerfl = 0.3932, which should be around 93%). Seems the code ddin't generate samples from posterior distributions.
It will be great if someone can take a look at my code and find the mistake inside. Much appreciated!
best,
