Archive for May, 2009


Bayesian Hierarchical Modeling

May 15, 2009

I’ve now dug a bit deeper into Carlin & Louis (2008) compliments of the end of the semester and have finally gotten around to Bayesian hierarchical modeling. The extensions of Bayesian models to more levels is very easy to do in WinBUGS (or OpenBUGS). Really the only thing you need to specify are priors and your model. Specification of priors seems to be a relatively simple thing to do in most cases. Priors for most model encountered in social sciences seem often follow Gamma distributions. Adding random effects to a model just requires additional priors for these variables (I believe they’re referred to as hyperpriors). So in the model you’ll have priors for precision (the reciprocal of variance), the random effects (often Gamma or Normal), and the fixed effects (often Normal or you could specify a flat or uniform prior). It’s important to try and pick a conjugate prior. Beyond priors, you just need to specify your model and finally initial values and WinBUGS is ready to roll. Don’t forget you have to burn in before you simulate. The authors seem to always run MCMC 1000 times before actual simulations.

Later today or in the next few days I plan to run a hierarchical model with achievement as my outcome, students and schools as random effects, and risk status and gender as my independent variables. I will post the WinBUGS code and discuss a little of the model selection done with DIC. Also I will play around with prior specification and see how that affects my estimates and Bayesian CIs.  Note: All the data shared here is simulated and does not represent any projects that I am currently involved.


Bayesian Longitudinal Modeling in R

May 1, 2009

A lot of my interests and the data I seem to get my hands on are longitudinal. I frequently use the lme4() package in R to do my analyses (with both discrete and continuous data). However, I recently learned about a package called MCMCglmm() written by Jarrod Hadfield, which allows a full Bayesian analysis on longitudinal discrete or continuous data. I strongly recommend this package and Jarrod is incredibly helpful when you encounter walls.

To use the package:
my.model <- MCMCglmm(y ~ 1 + predictor + time + time*predictor, random= ~ id, family="poisson", data=mydata, verbose=FALSE, DIC=TRUE)

This code will model my dependent variable (y) by time, my predictor, and a time-predictor interaction. I have allowed only random intercepts in my model (random = ~ id), and specified the poisson family. The more crucial bit of code is the DIC=TRUE, if you want to use DIC for model comparison.

To plot the posterior distributions of my parameters:

To see the posterior mode of the parameters:
posterior.mode(my.model$Sol[,1] ## For Intercept
posterior.mode(my.model$Sol[,2] ## For Predictor (lets assume 2 level dummy)
posterior.mode(my.model$Sol[,3] ## For Time
posterior.mode(my.model$Sol[,4] ## For Time*Predictor

Credible intervals can be obtained by the following:
HPDinterval(my.model$Sol[,1], 0.95) ## For Intercept … and so on …

For the DIC:

For more information about this package read the vignette, e.g. specifying priors, examining random effects.