Manually Calculating Deviance

Help
John
2011-08-04
2014-02-11
  • John
    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
    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
    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
    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.

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

         
        Last edit: Martyn Plummer 2014-02-11
        • Jared Harpole
          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:

          ####Beginning of R Syntax

          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

          #############END of R Syntax######################
          #############Beginning of JAGS Syntax

          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
          }
          
          #############End of JAGS Syntax

          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
            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
    John
    2011-08-10

    Thanks for the help Dr. Plummer and I apologize for causing any unwarranted
    frustration.

    Best Regards!!!

     
  • Jun Xu
    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
      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
    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
      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
    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!