Jonas Lindeløv
2014-04-03
Here's a very simple model which infers a population sd and effect size from one-sample data with known population mean. The effect size is the deviation of the sample from this population mean, standardized to SD's. I've used JZS-priors.
The problem is that "version 1" (with a small measurement error) somehow has very bad mixing while "version 2" (without measurement error) does not. If anyone can explain why this arises in version 1 but not version 2 - and maybe even fix it - that would be much appreciated.
# Priors on population parameters pop.pr ~ dgamma(0.00001, 0.00001) # approximation to Jeffrey's prior on variance pop.sd <- sqrt(1 / pop.pr) # convert to SD # Priors on (population) effect size effect.d ~ dt(0, 1, 1) # cauchy prior on standardized effect size (Cohen's d) effect.abs <- effect.d * pop.sd # convert to absolute effect size for(i in 1:length(x1)) { # Version with Bad mixing: with measurement error (SD=1) ability[i, 1] ~ dnorm(100, pop.pr) x2[i] ~ dnorm(ability[i, 1] + effect.abs, pow(1, -2)) #Version with relatively good mixing: without measurement error #x2[i] ~ dnorm(100 + effect.abs, pop.pr) }
I ran it on simulated data with mu=100, sd=15, d=2 (large effect) and n=100. With measurement error, there's a clear inverse relationship between estimated population sd and Cohen's d.
We're just going to use x2 for this simple simulation. It's a measurement after an intervention.
n = 100 pop.sd = 15 pop.mean = 100 effect.d = 2 x1 = as.vector(scale(rnorm(n))*pop.sd + pop.mean) # pretest (exact mean and sd) x2 = x1 + effect.d * pop.sd # posttest (noise-free effect)
All looks pretty good. Not perfect mixing on effect size and sd, but good enough.
Horrible mixing. It is even worse if the measurement error is decreased. That makes no sense to me since that should approximate the results from above.
By the way, my earlier post on another simple mixing problem remains unanswered: https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/aae85a9e/
Thanks for any help!
Martyn Plummer
2014-04-03
This is a parameterization problem. If you rewrite the model in the form below you get much better mixing:
ability[i, 1] ~ dnorm(100 + effect.abs, pop.pr) x2[i] ~ dnorm(ability[i, 1], pow(1, -2))
This is a weakness of sampling node-by-node. The sampling efficiency depends on how you parameterize your model.
Your example would work fine if you were using the glm module and if effect.abs
had a normal distribution. Then effect.abs
would be updated in a block with ability[,1]
and you would get good mixing, but the glm model cannot (yet) handle variables with a t distribution so this opportunity is missed.
Jonas Lindeløv
2014-04-04
Thanks, Martyn! That solved it for this particular problem, but mixing remained a problem when I tried to extend the model to repeated measures (x2 = x1 + effect). The whole purpose was to put the prior on Cohen's d rather than the absolute effect size, in order to do a JZS-prior. I figured out a solution which is mathematically equivalent to the script above but yields good mixing.
Instead of
# Priors on (population) effect size effect.d ~ dt(0, 1, 1) # cauchy prior on standardized effect size (Cohen's d) effect.abs <- effect.d * pop.sd # convert to absolute effect size
I did
# Priors on (population) effect size effect.abs ~ dt(0, pop.sd, 1) # cauchy prior on cohen's d of this effect effect.d <- effect.abs / pop.sd
Exactly why this works, I'm not sure. But the new model is simpler (cohen's d is just computed post-hoc), which is a nice thing.
Best,
Jonas