Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo


Bad mixing with measurement error

  • 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.

    Simulate data in R:

    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)

    Results without measurement error

    All looks pretty good. Not perfect mixing on effect size and sd, but good enough.
    Good mixing

    Result with measurement error

    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.
    Bad mixing

    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!

    Last edit: Jonas Lindeløv 2014-04-03
  • Martyn Plummer
    Martyn Plummer

    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.

  • 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.