R‎ > ‎R‎ > ‎

Bayesian Tobit Model

#The following R/OpenBUGS code runs a Bayesian Tobit-like model, with an arbitrary lower threshold
#Unlike a standard Tobit, this model outputs two separate sets of coefficients,
#the logit coefficients predict y>=Lower Threshold
#the normal coefficients predict y when y>Lower Threshold
#the second set of coefficients must be modified to reflect only the mean of the truncated portion of the model
#Like the Tobit, the coefficients of each part of the model cannot have opposite signs.

Code


tobitmodel <- function(){
  for (i in 1:J) {  #step through all the data
    z[i]<-step(Y[i]-LL) #LL is the lower limit, the truncation value
    z[i]~dbern(p[i])
    logit(p[i]) <- alpha0*X0[i] + alpha1*X1[i]  #Can change probit
   for (i in N+1:J){ 
    Y[i] ~ dnorm(mu[i],tau)%_%I(LL,) #This truncating adds a warning on the OpenBUGS module when the model is running
    mu[i] <- beta0*X0[i] + beta1*X1[i] #mu is normal mean, not truncated normal mean
  }
  
  #Prior on tau and regression coefficients
  tau ~ dgamma(0.001,0.001)
  alpha0 ~ dnorm(0,0.001)
  alpha1 ~ dnorm(0,0.001)
  beta0 ~ dnorm(0,0.001)
  beta1~dnorm(0,0.001)
    
  mzp <-1-mean(p[]) # Mean  zero prob
  sigma<-1/sqrt(tau)
}

# write the model code out to a file
setwd("~/Desktop/Size")
modelfile=function(){
  write.model(tobitmodel, "model1w.txt")
  model.file = paste(getwd(),"model1w.txt", sep="/")
  return(model.file)
}

#Initialize near the know values
inits=function(){
  list(alpha0=.9,alpha1=.7,beta0=3,tau=.2,beta1=2)
}

#Simulated data
Y=c(rep(0,1000),rtnorm(1000,6,1.79,upper=Inf,lower=LL),rtnorm(1500,4.5,1.1,upper=Inf,lower=LL))
X0=rep(1,length(Y)) #Intercept
X1=c(rep(0,150),rep(1,850),rep(1,1000),rep(0,1500))
J=length(Y)

#Input data into OpenBUGS
data=function(){
  list(Y=Y,X0=X0,J=J,LL=LL,X1=X1,N=N)
}

parameters.to.save = list("alpha0","alpha1","beta0","beta1","mzp","sigma") 

OpenBUGS.pgm="/Users/[accountname]/.wine/drive_c/Program Files/OpenBUGS/OpenBUGS322/OpenBUGS.exe"

library(arm)
WINE="/opt/local/bin/wine"
WINEPATH="/opt/local/bin/winepath"

Fit<- bugs(data(), inits, parameters.to.save=parameters.to.save, model.file = modelfile(),
           debug = T, DIC = F, digits = 5,n.chains=2,n.iter=2000,n.burnin=200,useWINE=T,
           OpenBUGS.pgm=OpenBUGS.pgm,WINE=WINE,WINEPATH=WINEPATH)

#mu is the mean of the UNTRUNCATED distribution
#to obtain the truncated mean, it needs to be adjusted by the following formula:
#Left censorsed std. normal, E(X) = density((a-mu)/sigma)/(1-cumdist((a-mu)/sigma))
#Compute the adjusted mean
mu.func <- function(mu,LL,sigma){
  value=mu+(dnorm((LL-mu)/sigma))/ (1 - pnorm((LL-mu)/sigma))*sigma
  return(value)
}

Comments