Brian Neelon R Programs
This page contains examples of R programs for running standard Bayesian regression models (linear regression, logistic/probit regression,
random effects models, latent class models, etc.). A basic working knowledge of R and Bayesian computation is
Note: Several of the programs below require you to download R packages such as mvtnorm, MCMCpack, etc.
General Disclaimer: These programs are meant as heuristic aids for newcomers to Bayesian inference,
not as exemplars of great code!
Please e-mail me if you find errors or have questions...thanks!
Bayesian Data Analysis by Andrew Gelman et al.
A First Course in Bayesian Statistical Methods by Peter Hoff
Bayesian Computation with R by Jim Albert
Finite Mixture and Markov Switching Models by Sylvia Fruhwirth-Schnatter
Normal and Linear Regression Models
Logistic, Multinomial and Probit Regression Models
- Normal Mean Gibbs.r: Gibbs sampler
for normal model with unknown mean and variance
- Normal Mean Metropolis.r:
Normal mean model with Metropolis update
- Multivariate Normal Mean Gibbs.r: Gibbs sampler for MVN model with unknown mean and variance
- Bivariate Normal.zip:
Two approaches to fitting a bivariate normal linear regression model. BivNorm condtional.r updates regression parameters from their univariate full conditionals. BinNorm Joint.r upates Y1,Y2 jointly -- this is easier code to follow, but much slower.
- Random Intercept Mixed Model.r:
Gibbs sampler for balanced-design mixed model with random intercept
- Random Slope Mixed Model.r:
Gibbs sampler for random slope model with unbalanced design. Centers out the random effects to improve
mixing of fixed effects parameters, beta.
- Here is another version which updates beta
after marginalizing out the random effects to further improve mixing.
See Chib and Carlin (1999) Statistics and Computing for more information.
- Note: both programs are slow as n increases, due to
inversion of large block diagonals. The programs invoke the spam package to try to improve computation a bit. More crudely, you can simply update
random effects one-at-a-time via a loop, as in this program.
This is acutally a bit faster than the matrix-inversion programs for n>200 or so. But it's not
the best solution.
- Random Slope Using MH: Uses MH for updating the random effects and is faster than the previous programs (although still not great). MH Function.r is the Metropolis function for udpating the random slopes, and Random Slope MH.r generates data and runs the MCMC.
- Skew Normal.r: Generates and implements Gibbs sampler for univariate skew-normal data. Full conditionals are here.
- Probit Model for Binary Data.r: Gibbs sampler for binary probit regression using Albert and Chib (1993) data-augmentation method.
- Probit Model for Ordered Categorical Data.r: Gibbs sampler for cumulative probit regression with 3-category ordered response variable.
- t-Link Logistic.r: Gibbs sampler for (approximate) logistic regression via t-link.
- Fixed Effects Logistic Regression.zip: Files for running
basic logistic regression using Metropolis step. Contains two files:
- Update Beta Fixed.r: function for updating beta (run this first)
- Log1.r: Generates data and runs MCMC (either source in "Update Beta" function or run "Update Beta" first)
- Multinomial Logit Model.r: Gibbs sampler
for multinomial logit model using Metropolis algorithm.
- Probit Random Intercept Model.r: Gibbs sampler
for probit random intercept model using data-augmentation approach
- Random Intercept Logistic Regression.zip: Files for running random intercept logistic regression using Gibbs and Metropolis steps. Contains three files:
- Update Beta Random.r: function for updating beta (run this first)
- Update b.r: function for updating random intercept b (run this second)
- Log2.r: Generates data and runs MCMC (either source in "Update Beta" and "Update b" functions or run those first)
- Bivariate Probit Model.zip: Various MCMC programs for fixed-effects bivariate (correlated) probit regression model. See "Read Me" file for details.
- Bivariate Probit Model with Random Effects.r: MCMC for random-intercept bivariate probit regression model.
Bayesian Two-Part Latent Class Model
- Zero-Inflated Poisson.r: Fixed effects zero-inflated Poisson (ZIP) regression using MH.
- Spatial Poisson Hurdle Model.zip: R2WinBUGS and BUGS code for simulating and fitting the spatial Poisson hurdle model described here. See README.txt for file discriptions.
Finite Mixture and Dirichlet Process Models
- Bayesian Two-Part.zip: R code for running
latent-class two-part model as described in
Neelon et al. (2011), with three latent classes and random intercepts within classes. Contains 8 files:
- Three Class MCMC.r: This is the main program for implementing the MCMC algorithm. It calls the other programs,
but, of course, you'll need to change the file path names;
- Update Gamma.r, Update Alpha.r, etc: These are the five functions used in the Gibbs and MH steps;
- simdata.txt: Data file
- Make simdata.r: Creates simdata.txt
- Growth Mixture Model for Blood Pressure and Birth Outcomes: Contains R code for simulating and fitting the Bayesian growth mixture model (GMM) desribed in Neelon et al. (2011) A Bayesian Growth Mixture Model to Examine Maternal Hypertension and Birth Outcomes. See README.txt for file discriptions.
Latent Factor Model
- Two-Component Finite Mixture.r: Gibbs sampler for
two-component (or "two-class") finite mixture normal model with class-specific means and variances.
Dirichlet prior on class-membership probabilities.
The program also illustrates the so-called "label-switching" problem and uses order-constraint solution to resolve
the issue. See Fruhwirth-Schnatter (2006) for a more general discussion.
- Two-Class Bivariate Normal Mixture.zip: Zip file contains 5 programs: 1) makedata.r creates two-class bivariate normal data with a class-level covariate ("male"); 2) data.txt is the data file; 3) MCMC.r runs the mcmc using the conditional udpating scheme described above in BivNorm.zip; programs 4 and 5 are Update C.r and Update Gamma.r---these are functions for updating the latent class indicators and the class-specific regression parameters (gamma's) via MH.
- Three-Component Dirichlet Process Mixture Model.r: Uses a stick-breaking DP mixture model to estimate a three-component mixture of normals. For a nice review of DPs, and Bayesian nonparametrics in general, see this
paper by David Dunson.
Bayesian Isotonic Regression
- Bayesian Latent Factor Model.r: Three-factor linear model with M=6 manifest variables. See Lopes and West (Statistica Sinica, 2004) for details.
Spatial Analysis for Areal Data
- Bayesian Order-Restricted Inference for Categorical Predictors: Contains files for implemented Bayesian order-restricted inference on categorical predictors, as described in Dunson and Neelon (2003). The ZIP file contains two R programs: the first is the minmax function used for imposing the order restriction, and the second simulates a linear model with three categorical covariates and applies the minmax monotonicity restriction as part of a Gibbs sampling algorithm.
- Bayesian Isotonic Regression for Continuous Predictors: The program first
generates normal data with a montone-increasing sigmoidal mean function. It then fits a Bayesian, piecewise-linear
isotonic regression model with a fixed number of "free" knots (i.e., unknown knot locations).
This approach provides a reasonable estimate the mean function, but it assumes that the number of knots is known;
to account for uncertainty in the number of knots (K), one can re-fit the model under different choices of K
and use model averaging. Alternatively, one can use reversible jump MCMC or a model comparison criterion,
such as BIC or DIC, to choose the "optimal" number of knots.
The method is outlined in this working paper.
See Denison, Mallick and Smith (1998) and Holmes and Heard (2003) for a related approaches.
- Normal ICAR Model.r: Generates and fits data for
spatially varying normal data under intrinsic CAR (ICAR) model. To generate data, we set spatial smoothing parameter close
to 1 (but not equal to 1) to ensure that the joint covariance of the spatial random effects is non-singular. The full condtionals are given here.
- CAR_Spam.r: ICAR Model with sparse matrix (SPAM) routines. Much faster!! Thanks to Howard Chang for pointing this out.
- Bivariate CAR (biCAR) Model.zip: Generates and fits data for spatially varying bivariate normal data assuming bivariate intrinsic CAR (biCAR) random effects. The Gibbs sampling algorithm and data generation are in bicar.r, and the function for updating the spatial effects phi are in update phi.r. The full condtionals are given here.
- Spatial Bivariate Probit (Biprobit) Model.zip: Generates and fits data for correlated binary data as described in Neelon et al. (2012) A spatial bivariate probit model for correlated binary data with application to adverse birth outcomes. See Read Me.txt for details.
- Spatial Mixture.zip: Generates and fits data for spatial mixture model as described in Neelon et al. (2014) A Multivariate Spatial Mixture Model for Areal Data: Examining Regional Differences in Standardized Test Scores. To appear in JRSS-C. See READ ME.txt for details.
More to come...
- Linear Quantile Regression.r: Quantile regression for simple linear regerssion model. For more details, see Reed and Yu "A Partially Collapsed Gibbs Sampler for Bayesian Quantile Regression".
Return to Brian's home page