Bayesian Tobit Model

ex bayesian tobit.r

#   ------------------------------------------------------------------
#  | PROGRAM NAME: ex_bayesian_tobit
#  | DATE: 9/17/11
#  | CREATED BY:  Matt Bogard
#  | PROJECT FILE: http://econometricsense.blogspot.com/2011/09/bayesian-models-with-censored-data.html             
#  |----------------------------------------------------------------
#  | PURPOSE: comparison of models for censored dependent variables              
#  | 1 - least squares
#  | 2 - tobit model
#  | 3 - bayesian model
#  |------------------------------------------------------------------
#  | REFERENCES: 
#  | UCLA Statistical Computing: http://www.ats.ucla.edu/stat/R/dae/tobit.htm
#  | R Package 'MCMCpack' documentation : # http://mcmcpack.wustl.edu/documentation.html
#  | 
#  | Literature:
#  |   Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. “MCMCpack: Markov Chain Monte Carlo in R.”, 
#  |   Journal of Statistical Software. 42(9): 1-21. http://www.jstatsoft. org/v42/i09/.
#  |   
#  |   Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. Scythe Statistical Library 1.0.
#  |   http://  scythe.wustl.edu.
#  |
#  |   Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002. Output Analysis and Diagnos- tics for 
#  |   MCMC(CODA). http://www-fis.iarc.fr/coda/.
#  |
#  |  Siddhartha Chib. 1992. “Bayes inference in the Tobit censored regression model." Journal of Econometrics. #  |  51:79-99.
#  |
#  |
#  |
#   ------------------------------------------------------------------


# example tobit model
# get data

mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/dae/tobit.csv"))

#explore dataset

names(mydata) # list var names

dim(mydata) # data dimensions

hist(mydata$apt)  # histogram of dependent variable for academic aptitude

#  indcates right or upper bound censoring at 'y' = 800
# run model using standard ols regression

ols <- lm(mydata$apt~mydata$read + mydata$math + as.factor(mydata$prog))

summary(ols)

#  tobit model 
library(VGAM) # load package
tobit <- vglm(mydata$apt ~ mydata$read + mydata$math + as.factor(mydata$prog), tobit(Upper=800))

summary(tobit)

# note the coefficients for the tobit model are larger, indicating the downward bias
# of the OLS estimates

# bayesian model

library(MCMCpack) # load package

bayes.tobit <- MCMCtobit(mydata$apt ~ mydata$read + mydata$math + as.factor(mydata$prog), above = 800, mcmc = 30000, verbose = 1000)
summary(bayes.tobit)

plot(bayes.tobit)

# the empirical (posterior mean) looks very similar to the tobit estimates.