## pD half the effective number of parameters?

Help
2011-03-22
2012-09-01
• Tim Handley - 2011-03-22

In the DIC, the penalty term pD is generally described as being an 'effective
number of parameters.' However, in my recent modelling, the value of pD has
generally come out half as large as I expected. As an experiment, I ran the
following simple bilinear normal JAGS model using Rjags. (model text and R
code below). I'm pretty confident that this model has four parameters (one
intercept, two slopes, and a precision). Oddly, the call to DIC.out reports a
penalty of 2.0, suggesting that this model effectively has just two
parameters. Is there a missing factor of two in dic.samples? Am I
misinterpreting the meaning of pD? Are there really correlations among these
parameters such that the four actual parameters only contain the same
information as two independent parameters? Any comments would be greatly
appreciated. Thanks,

## JAGS model file

model{

b ~ dunif(-100,100)

m1 ~ dunif(-100,100)

m2 ~ dunif(-100,100)

tau ~ dgamma(.01,.01)

for(i in 1:Npoints)

{

y_~dnorm(b+m1x1_+m2x2_, tau)

}

}##End model

# RCode

library(rjags)

b_true = 20

m1_true=5

m2_true=12

tau_true = 1

Npoints=1e2

Nchains=3

x1=runif(Npoints,-100,100)

x2=runif(Npoints,-100,100)

y = b_true+m1_truex1+m2_truex2 + rnorm(Npoints, mean=0, sd=1/sqrt(tau_true))

data_list = list(x1,x2,y,Npoints)

names(data_list)<-c("x1","x2","y","Npoints")

jagsmodel= jags.model(

file = ".\SimpleBilinear.txt",

data = data_list,

n.chains=Nchains,

coda.out = coda.samples(jagsmodel, c("b","m1","m2","tau"), n.iter=20e3,
thin=1)

DIC.out = dic.samples(jagsmodel, n.iter=40e3, thin=1, type="pD")

• Tim Handley - 2011-03-22

Update: I ran this same model in OpenBUGS, which reported the same deviance
numbers, but a different pD and DIC. Since OpenBUGS reported the expected
value of pD, this suggests that there is indeed an error in the JAGS
calculation of pD and DIC.

JAGS results:

Mean deviance: 295.4

penalty 2.045

Penalized deviance: 297.5

OpenBUGS results:

Dbar: 295.4

Dhat:291.4

pD: 4.0

DIC: 299.4

• Martyn Plummer - 2011-03-23

Oh dear. This bug will occur in any model with normal, Poisson, or binomial
observed nodes. I put in a fast method of calculating pD for these
distributions, but only did half the calculations. (The value of pD is
calculated from the undirected Kullback-Liebler divergence between the
parallel chains, which is the sum of the two directed Kullback-Leibler
divergences: only one divergence is currently included)

This is why we need more regression tests.

• Tim Handley - 2011-03-23

Again, thanks. I feel better knowing that I wasn't imagining things.

• Tim Handley - 2011-03-24

Martyn, I'm not sure how to proceed from here. I have some fairly different
models which I need to compare I think the most straightforward way to do that
is to use DIC, and the most straightforward way to get DIC values would be to
use dic.samples with a pD penalty. However, depending on time, I may need to
work out something else. Do you have a rough idea of when you might release an
update which would fix this issue? Alternatively, does this issue affect the
calculation of the optimism (popt) penalty?

• Martyn Plummer - 2011-03-28

I generally release a new version o f JAGS after each release of R. The next
scheduled release of R is 13 April (R 2.3.0), so you shouldn't have to wait
long.

• Anonymous - 2011-05-06

I hate to be this forward because I know you probably have lots of other stuff
on your plate, but I was wondering if you had any idea of the timeline on
which this issue is likely to be fixed. I'm about to embark on some long-
running estimations and would rather launch them with a version of JAGS that
has a working dic.samples.

Thanks!