2

I'm doing a network meta-analysis including several clinical trials. The response is binomial. Each trial contains several treatments.

When I do a random effects model, the output from JAGS and WinBUGS is similar. When I do a fixed effects model, the DIC and pD components are way out, though the posteriors of the parameters I'm interested in are similar.

I've got similar models that have Gaussian response, not binomial, and JAGS and WinBUGS are in agreement.

The BUGS/JAGS code for the fixed effects model is lifted from page 61 of this and appears below. However, the same code runs and produces similar posteriors using WinBUGS and JAGS, it's just DIC and pD that differ markedly. So I don't think this code is the problem.

for(i in 1:ns){  # Loop over studies
  mu[i] ~ dnorm(0, .0001)
    # Vague priors for all trial baselines
  for (k in 1:na[i]) {  # Loop over arms
    r[i, k] ~ dbin(p[i, k], n[i, k])
      # binomial likelihood
    logit(p[i, k]) <- mu[i] + d[t[i, k]] - d[t[i, 1]]
      # model for linear predictor
    rhat[i, k] <- p[i, k] * n[i, k]
      # expected value of the numerators
    dev[i, k] <- 
      2 * (r[i, k] * (log(r[i, k]) - log(rhat[i, k])) +     
      (n[i, k] - r[i, k]) * (log(n[i, k] - r[i, k]) +
      - log(n[i, k] - rhat[i, k])                         ))
      # Deviance contribution
  }
  resdev[i] <- sum(dev[i, 1:na[i]])
    # summed residual deviance contribution for this trial
}
totresdev <- sum(resdev[])
  # Total Residual Deviance

d[1] <- 0
  # treatment effect is zero for reference treatment
for (k in 2:nt){
  d[k] ~ dnorm(0, .0001)
}  # vague priors for treatment effects

I found an old post describing a known issue, but that's too old for me to think it's the same problem.

Are there any known issues with JAGS reporting the wrong DIC and pD? (Searching for "JAGS bugs" isn't all that helpful.)

I'd be grateful of any pointers.

Marcelo Ventura
  • 581
  • 7
  • 19
harrys
  • 179
  • 2
  • 9

1 Answers1

3

There are a number of different ways to calculate pD, and the method used by JAGS differs to that used by WinBUGS. See the Details section of the following help file:

?rjags::dic

Specifically:

DIC (Spiegelhalter et al 2002) is calculated by adding the “effective number of parameters” (pD) to the expected deviance. The definition of pD used by dic.samples is the one proposed by Plummer (2002) and requires two or more parallel chains in the model.

The details are in the (lengthy, but well worth reading) discussion of the following paper:

Spiegelhalter, D., N. Best, B. Carlin, and A. van der Linde (2002), Bayesian measures of model complexity and fit (with discussion). Journal of the Royal Statistical Society Series B 64, 583-639.

In general: all methods to calculate pD (of which I am aware) are approximations, and if two such methods disagree then it is probably because the assumptions behind for one (or both) approximation(s) are not being met.

Matt

Matt Denwood
  • 2,537
  • 1
  • 12
  • 16