Brian Neelon R Programs

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

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!

Useful References
    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

R Programs

Normal and Linear Regression Models
  1. Normal Mean Gibbs.r: Gibbs sampler for normal model with unknown mean and variance

  2. Normal Mean Metropolis.r: Normal mean model with Metropolis update

  3. Multivariate Normal Mean Gibbs.r: Gibbs sampler for MVN model with unknown mean and variance

  4. Bivariate 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.

  5. Random Intercept Mixed Model.r: Gibbs sampler for random intercept mixed model with unbalanced design. Here's a slightly revised version that allows for larger n by avoiding random effect design matrix.

  6. 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's a much faster version that uses a conditional update of the random effects: b1|b2 and b2|b1. This is computationally more efficient since it relies on univariate updates of the random effects.This is to appear in a forthcoming paper on improved MCMC for multivariate random effects models.
    • Here's 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.
    • Random Slope Using MH: Uses MH for updating the random effects. MH Function.r is the Metropolis function for udpating the random slopes, and Random Slope MH.r generates data and runs the MCMC.

  7. Skew Normal.r: Generates and implements Gibbs sampler for univariate skew-normal data. Full conditionals are here.

Logistic, Multinomial and Probit Regression Models
  1. Probit Model for Binary Data.r: Gibbs sampler for binary probit regression using Albert and Chib (1993) data-augmentation method.

  2. Logistic Regression for Binary Data: Gibbs sampler for (approximate) logistic regression via t-link.

  3. MH Logistic Files for running basic logistic regression using Metropolis step. Contains two files:
    1. Update Beta Fixed.r: function for updating beta (run this first)
    2. Log1.r: Generates data and runs MCMC (either source in "Update Beta" function or run "Update Beta" first)

  4. Probit Model for Ordered Categorical Data.r: Gibbs sampler for cumulative probit regression with 3-category ordered response variable.

  5. Multinomial Logit Model.r: Gibbs sampler for multinomial logit model using Metropolis algorithm.

  6. Probit Random Intercept Model.r: Gibbs sampler for probit random intercept model using data-augmentation approach.

  7. Logistic Random Intercept Model.r: Gibbs sampler for logit random intercept model using data-augmentation approach with t-link.

  8. MH Logistic Random Intercept Files for running random intercept logistic regression using Gibbs and Metropolis steps. Contains three files:
    1. Update Beta Random.r: function for updating beta (run this first)
    2. Update b.r: function for updating random intercept b (run this second)
    3. Log2.r: Generates data and runs MCMC (either source in "Update Beta" and "Update b" functions or run those first)

  9. Bivariate Probit Various MCMC programs for fixed-effects bivariate (correlated) probit regression model. See "Read Me" file for details.

  10. Bivariate Probit Model with Random Effects.r: MCMC for random-intercept bivariate probit regression model.

Count Models
  1. Zero-Inflated Poisson.r: Fixed effects zero-inflated Poisson (ZIP) regression using MH.

  2. Spatial Poisson Hurdle R2WinBUGS and BUGS code for simulating and fitting the spatial Poisson hurdle model described here. See README.txt for file discriptions.

Bayesian Two-Part Latent Class Model
  1. Bayesian 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:
    1. 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;
    2. Update Gamma.r, Update Alpha.r, etc: These are the five functions used in the Gibbs and MH steps;
    3. simdata.txt: Data file
    4. Make simdata.r: Creates simdata.txt

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

Finite Mixture and Dirichlet Process Models
  1. 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.

  2. Two-Class Bivariate Normal 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; 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.

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

Latent Factor Model
  1. Bayesian Latent Factor Model.r: Three-factor linear model with M=6 manifest variables. See Lopes and West (Statistica Sinica, 2004) for details.

Bayesian Isotonic Regression
  1. 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.

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

Spatial Analysis for Areal Data
  1. 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.

  2. CAR_Spam.r: ICAR Model with sparse matrix (SPAM) routines. Much faster!! Thanks to Howard Chang for pointing this out.

  3. Bivariate CAR (biCAR) Generates and fits data for spatially varying bivariate normal data assuming bivariate intrinsic CAR (biCAR) random effects. The full condtionals are given here. Note that this program uses conditional updates for the bivariate spatial effects: phi1|phi2 and phi2|phi1. This permits univariate updates of the spatial effects, thus substantially increasing computational speed. This will appear in a forthcoming paper on improved MCMC for multivariate random effect models.

  4. Spatial Bivariate Probit (Biprobit) 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.

  5. Spatial 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. A new version with computationally efficient univariate updates for the spatial effects is forthcoming.

Quantile Regression
  1. 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".

More to come...

Return to Brian's home page