John
2011-08-04
I've been trying to figure out how JAGS calculates the deviance at each
iteration during the sampling process from the posterior. I used the
parameterization for the normal distribution detailed in the JAGS manual
within my model and calculated the deviance manually for a simple linear
regression. However, what I calculate and what JAGS calculates are slightly
different. Any help on what I'm doing wrong would be greatly appreciated.
Thanks in advance.
{ model{ for(i in 1:n){ y[i] ~ dnorm(mu[i],tauy) mu[i] <- beta[1] + beta[2]*x[i] D[i] <- log(pow((tauy/(2*pi)),0.5)*exp(-(pow((x[i]-mu[i]),2)*tauy))) } for(i in 1:2){ beta[i] ~ dnorm(0.0,1.0E-4) } tauy ~ dgamma(0.1,0.1) Deviance <- -2*sum(D[]) } }[\code]
Martyn Plummer
2011-08-04
Try this:
D[i] <- log(pow((tauy/(2*pi)),0.5)*exp(-[b]0.5*[/b](pow((x[j]-mu[j]),2)*tauy)))
John
2011-08-04
Hey Dr. Plummer,
I tried your suggestion, but there is still a difference in the deviance
calculated by JAGS and what is being calculated manually within the model.
Here is an example of the output:
Deviance
deviance
861.4990
593.4112
844.0658
596.6953
935.9164
593.5914
Deviance calculated manually is in the left column, while deviance calulated
by JAGS is in the right. I'm not sure what is going wrong.
Martyn Plummer
2011-08-05
OK then try this.
D <- log(pow((tauy/(2*pi)),0.5)*exp(-0.5*(pow(([b]y[j][/b]-mu[j]),2)*tauy)))
That's two mistakes in one formula, and now I've wasted half an hour double
checking my code.
You should be working on a log scale exclusively to avoid numerical problems.
I used this
D[j] <- - log(tauy) + log(2*pi) + pow(y[j]-mu[j],2) * tauy Deviance <- sum(D)
and got perfect agreement.
Mauricio Garnier-Villarreal
2014-02-02
Hi
This thread was extremely helpful, currently I am working on the calculation of the WAIC, and it requires the calculation of the log likelihood for each data row. This could also be calculated outside of jags, like in R, but i want to calculate it inside the mcmc
I have tried to adjust this code you included here, but i ahvent been able to adapt it properly because my model is multivariate, I would appreciate your help on adapting the deviance function for a multivariate model
My model is as follow, you can see that I have try to adjust the formula, but it doesnt get to the same deviance as jags
model{ for(i in 1:nRows){ Y[i, 1:nVars] ~ dmnorm(Y.hat[i, 1:nVars], TauY[1:nVars,1:nVars]) Y.hat[i,1:nVars] <- b[id[i],1,1:nVars] + age[i]*b[id[i],2,1:nVars] + (b[id[i], 3, 1:nVars]*((age[i] - b[id[i],4,1:nVars])*step(age[i]-b[id[i],4,1:nVars]))) for(j in 1:nVars){ log_lik01[i,j] <- (- log(TauY[j,j]) + log(2*pi) + pow(Y[i,j]-Y.hat[i,j],2) * TauY[j,j])/(-2) } log_lik[i] <- sum(log_lik01[i,1:nVars]) } for(p in 1:nVars){ for(j in 1:N){ b[j, 1:params,p] ~dmnorm(mub[1:params,p], Taub[1:params, 1:params,p]) } mub[1:params,p] ~ dmnorm(pm[1:params,p], prec[1:params, 1:params]) Taub[1:params, 1:params,p] ~ dwish(R[1:params,1:params], params) Sigma2b[1:params,1:params,p] <- inverse(Taub[1:params,1:params,p]) for (i in 1 : params) {Stdb[i, p] <- sqrt(Sigma2b[i, i, p])} } log_lik0 <- sum(log_lik) dev <- -2*sum(log_lik) for (j in 1:nVars) {stdy[j] <- sqrt(SigmaY[j,j])} TauY[1:nVars,1:nVars] ~ dwish(W[1:nVars,1:nVars], nVars) SigmaY[1:nVars,1:nVars] <- inverse(TauY[1:nVars,1:nVars]) }
Thank you very much
Martyn Plummer
2014-02-07
[Edit: I have edited and expanded this response after mistakes highlighted by Jared below]
The correct formula for the deviancelog density is
log_lik[i] <- logdet(TauY)/2 - nVars * log(2*pi)/2 - t(Y[i,] - Y.hat[i,]) %*% TauY %*% (Y[i,] - Y.hat[i,])/2
The code in the "bugs" module that calculates the deviance uses the non-normalized version of the log density:
logdet(TauY)/2 - t(Y[i,] - Y.hat[i,]) %*% TauY %*% (Y[i,] - Y.hat[i,])/2
i.e. it is missing the term nVars * log(2*pi)/2
. Of course this term is constant so it won't make a difference when you look at the difference between two models.
Strictly speaking, this is not a bug, since the deviance is defined only up to an additive constant. Nevertheless, I have corrected the code in the development version (future JAGS 4.0.0) so that it uses the fully normalized version.
I also added a feature in JAGS 4.0.0 that generates a log-density function for each distribution. So in future you will be able to write
log_lik[i] <- logdensity.mnorm(Y[i,], Y.hat[i,], TauY)
Finally, I note that the formula for the density of the multivariate normal distribution in the JAGS User Manual is also wrong. I have corrected this in the development version.
Jared Harpole
2014-02-10
Dr. Plummer I am having a similar problem and created an example to demonstrate the issue.
Here is the code to for the example:
library("R2jags")
load.module("glm")
load.module("bugs")
load.module("msm")
reps <- 20000
burn <- floor(reps/2)
age <- matrix(c(rep(8, 20), rep(8.5, 20), rep(9, 20), rep(9.5, 20)), ncol = 4)
inits <- list(beta0 = 40, beta1 = 1)
datList <- list(M = 4, N = 20, pi=pi, Y = structure(
.Data = c(47.8, 48.8, 49.0, 49.7,
46.4, 47.3, 47.7, 48.4,
46.3, 46.8, 47.8, 48.5,
45.1, 45.3, 46.1, 47.2,
47.6, 48.5, 48.9, 49.3,
52.5, 53.2, 53.3, 53.7,
51.2, 53.0, 54.3, 54.5,
49.8, 50.0, 50.3, 52.7,
48.1, 50.8, 52.3, 54.4,
45.0, 47.0, 47.3, 48.3,
51.2, 51.4, 51.6, 51.9,
48.5, 49.2, 53.0, 55.5,
52.1, 52.8, 53.7, 55.0,
48.2, 48.9, 49.3, 49.8,
49.6, 50.4, 51.2, 51.8,
50.7, 51.7, 52.7, 53.3,
47.2, 47.7, 48.4, 49.5,
53.3, 54.6, 55.1, 55.3,
46.2, 47.5, 48.1, 48.4,
46.3, 47.6, 51.3, 51.8),
.Dim = c(20, 4)),
age = age,
R = structure(
.Data = c(1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1),
.Dim = c(4, 4)))
param2 <- c("log_lik0","dev")
jagsOUT_ll <- jags(data=datList,inits=NULL,param=param2,
n.chains=2,n.iter=reps,n.burnin=burn,n.thin=1,
model.file="test.txt")
jagsOUT_ll
model
{
beta0 ~ dnorm(0.0, 0.001)
beta1 ~ dnorm(0.0, 0.001)
for (i in 1:N) {
Y[i, 1:M] ~ dmnorm(mu[i,1:M], Omega[ , ])
mu[i,1:M] <- beta0 + beta1*age[i,1:M]
}
Omega[1 : M , 1 : M] ~ dwish(R[ , ], 4) Sigma[1 : M , 1 : M] <- inverse(Omega[ , ]) for(i in 1:N){ log_lik[i] <- - M * log(2*pi)/2 + logdet(Omega)/2 + t(Y[i,] - mu[i,]) %*% Omega %*% (Y[i,] - mu[i,])/2 } log_lik0 <- sum(log_lik[]) dev <- -2*log_lik0 }
First, I change the log_lik[i] <- -logdet(Omega)/2 + M * log(2pi)/2 -
t(Y[i,] - mu[i,]) %% Omega %*% (Y[i,] - mu[i,])/2
(This is from the post above)
to
log_lik[i] <- - M * log(2pi)/2 + logdet(Omega)/2 + t(Y[i,] - mu[i,]) %% Omega %*% (Y[i,] - mu[i,])/2.
When I did it using the first way the Deviance is negative (see Example below) but when I do it the second way everything looks good.
Is the second way the correct way to do it?
Second, when I compute the deviance in the code above (dev) my deviance is still a bit off of what JAGS is doing. See below:
jagsOUT_ll
Inference for Bugs model at "test.txt", fit using jags,
2 chains, each with 20000 iterations (first 10000 discarded)
n.sims = 20000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
dev 167.009 25.480 115.840 150.307 167.526 184.236 215.585 1.001 20000
log_lik0 -83.505 12.740 -107.793 -92.118 -83.763 -75.154 -57.920 1.001 20000
deviance 208.523 5.263 200.904 204.891 207.796 211.398 219.871 1.001 2900
Here the deviance is I calculated is still off by quite a bit and I don't see what I am doing wrong. Any assistance would be helpful.
(NEGATIVE DEVIANCE EXAMPLE)
When I use the log_lik[i] <- -logdet(Omega)/2 + M * log(2pi)/2 -
t(Y[i,] - Y.hat[i,]) %% Omega %*% (Y[i,] - Y.hat[i,])/2
to calculate the deviance I get the following:
jagsOUT_ll
Inference for Bugs model at "test.txt", fit using jags,
2 chains, each with 20000 iterations (first 10000 discarded)
n.sims = 20000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
dev -166.986 25.832 -215.604 -184.451 -167.465 -149.832 -114.730 1.001 5900
log_lik0 83.493 12.916 57.365 74.916 83.733 92.226 107.802 1.001 6900
deviance 208.565 5.237 200.996 204.992 207.832 211.344 219.922 1.001 5600
Which has a negative deviance (I don't think that should be negative).
Martyn Plummer
2014-02-11
Sorry, there was a sign error in my answer above, which I have now corrected by editing my post. You should use
for(i in 1:N){ log_lik[i] <- - M * log(2*pi)/2 + logdet(Omega)/2 - t(Y[i,] - mu[i,]) %*% Omega %*% (Y[i,] - mu[i,])/2 } log_lik0 <- sum(log_lik[]) dev <- -2*log_lik0
If you want to get the same results as the deviance function in JAGS 3.4.0 then drop the term in pi:
for(i in 1:N){ log_lik[i] <- logdet(Omega)/2 - t(Y[i,] - mu[i,]) %*% Omega %*% (Y[i,] - mu[i,])/2 } log_lik0 <- sum(log_lik[]) dev <- -2*log_lik0
John
2011-08-10
Thanks for the help Dr. Plummer and I apologize for causing any unwarranted
frustration.
Best Regards!!!
Jun Xu
2014-02-07
I believe my question is kind of related. So I just put it here. I am trying to get DIC's for several models that I am estimating. I am not sure whether there is any guideline as to how many iterations I should have for dic.samples(). Any suggestion?
Martyn Plummer
2014-02-07
You can set a trace monitor for the deviance and look at its mixing properties. After loading the "dic" module, just monitor "deviance" as if it were a node in your model (NB If you do actually have a node called "deviance" in your model then the monitor will be set to that node, otherwise it will calculate the deviance for you).
Jun Xu
2014-02-07
In addition, I also had the following warning messages after I run dic.samples. Would these be problems? What would be causes of these missing?
Warning messages: 1: In dic.samples(jagsModel, 1000, "pD") : NAs introduced by coercion 2: In dic.samples(jagsModel, 1000, "popt") : NAs introduced by coercion
Martyn Plummer
2014-02-07
I don't know. But I can tell you the message is coming from R not from JAGS.
Jun Xu
2014-02-09
Hi,
My understanding is that dic involves numerical integration of log likelihood (the d bar part). So is my choice of 1000 (number of iteration in the dic.samples function) reasonable regardless of my initial sample size (over 10000)? I did follow your advice to load the dic module (load.module("dic", quiet=T) right after I require the rjags package, and I added a node for deviance to monitor it. Deviance seems to mix well. So is it ok to use 1000 iteration for calculating DIC in the dic.samples function? Sorry for my confusion and thanks a lot for your help!