Hi, Kalinda!

I'm having some trouble understanding the question. You said you get different DIC values depending on the level you calculate it. But, according to Gelman's et al. "Bayesian Data Analysis", de Deviance Information Criterion is calculated as

DIC = 2* D^avg(y) - Dpars_hat(y)

where the first term in the left side is the mean deviance, -2 * log p(y|pars) averaged over the posterior distribution; and the second term is the deviance at the expected parameter value. This is the same as

DIC = D^avg(y) + pd

where pd is the effective number of parameters, which takes into account the partial pooling of estimates in hierarchical models with random coefficients.

As the DIC depends on the log-likelihood of the outcome variable, I don't understand how you are computing it at various levels of the model. Could you explain how you are computing the DICs for your model?

Hi Erikson,

Thanks for your response.  For my random coefficients probit (hierarchical) model, I calculated the DIC using Spiegelhalter's approach, but my question would be the same even with Gelman's .  With Spiegelhalter's, I used p_d=D^avg(y)-- Dpars_hat(y).   What I mean by levels is that I calculated at two different "levels" of the model.  

My utility function in my model is more or less structured as follows:

U_I = a_i + b_iX + e (Level 1)

b_i = g_0 +g_1Z+u  (Level 2)

For level 1, I calculated the DIC for each of the 404 individuals that I have in my hierarchical data, using each individual's own estimated coefficients (b_i) to calculate log-likelihood and the deviance statistics, and summed these individual DICs up to represent the DIC for the overall hierarchical model.   For the second approach, I calculated one DIC from the coefficients estimated from my Z matrix (level 2 predictor) which would represent the coefficients for the average individual in my data.  This would be similar to how it would be calculated for  a model without hierarchy, right?

Based on the first calculation, it suggests that my hierarchical model is much worse than the plain probit model (which doesn't make sense to me why this is.  My hit rate is much better with the hierarchical model).   Based on the second approach, it suggests that my hierarchical model is much better.   

I don't know which approach is the right way to calculate DIC for a hierarchical model.  

I hope this clarifies what I am asking for guidance on.

Kalinda

Kalinda,

You shouldn't compute the DIC for each outcome in your dataset. The DIC is a measure of model fit based on the full y vector. If you compute the DIC for each outcome and then add them together, your pd will be multiplied by N! This might explain why the DIC for the hierarchical model is much higher.

This is how I understand you are calculating the DIC for the first level (it might change depending which formula you are using, but the idea remains the same):

DIC = sumN(-2 log p(yi|parsi^) + 2pd)

       = 2Npd + sumN(-2 log p(yi|parsi^))

This is how it should be done:

DIC =  sumN(-2 log p(yi|parsi^)) + 2pd

Which means you should add together the deviance for each outcome and only then sum the effective number of parameters.

I think that the DIC calculated at the second level doesn't make much sense - I still don't understand how you are doing it. What might make sense is to compute the likelihood based on the second level parameters, changing the "focus" of the model, as Spiegelhalter puts it.

What are you using to fit your model? BUGS and JAGS have built-in DIC computation.

Hi Erikson,

I built my model in R (wanted more control than BUGS/JAGS allows).   So, I am calculating DIC directly in the R code.  

I'm unclear on how you would calculate 2pd in the response you gave.  This is one of the sources of my confusion.  Should not 2pd = 2*(D^avg(yi)-- Dpars_hat(yi)) for all i because each individual i has their own set of parameters, so should not pd be summed across all i?

For the second level, could you elaborate on how you "compute the likelihood based on the second level parameters, changing the "focus" of the model"?   Specifically, what parameters are included in this?  Should it include the parameters for the interaction terms, or just the main effects?  This is the second source of my confusion.

Kalinda,

I misunderstood your calculations earlier. If you computed pd(i) for each outcome i and then summed them together, you obtained the correct value for pd.

But there's an important issue here: pd is not a measure for each outcome, but a property of the model. It's a measure of effective number of parameters, taking into account the partial pooling of coefficient estimates in multilevel models.

Using a notation similar to R functions:

pd = meanS(sumN(-2 * log p(yi | parss))) - sumN(-2 log p(yi|parsi^)

Where:

• parss is a vector of simulated values of pars (i.e., in a MCMC chain);
• sumN(-2 * log p(|)) is the model deviance, computed for each simulated parameter value and averaged (first term in the right side) and for the average parameter value (second term in the right side).

The problem might be the choice of the likelihood function, p(y|pars). What you choose as pars in the likelihood function changes the "focus" of the model. Spielgelhalter has some nice slides with examples of model focus in the attached link.

Simply put, the model focus is the level of the model for which you want to predict new cases. In your case, for the first level (responses within individuals?):

p(yi|ai, bi, se) = dnorm(yi, mean=ai+bi*xi, sd=se)

Where each parameter is either the posterior expected value (par^) or a simulated posterior value from the MCMC chain (pars). In case you want to focus on the second level (new individuals?), you must integrate over the first level parameters that depend on the second level. This can be achieved with:

p(yi|ai, g0, g1, se, sb) = dnorm(yi, mean=ai+(g0+g1*zi)*xi, sd=sqrt(se^2 + sb^2))

http://www.mrc-bsu.cam.ac.uk/wp-content/uploads/DIC-slides.pdf