Hi all,

I am trying to set up simple tobit model that should evolve toward a

hierarchical structure later. Unfortunately, even in its basic incarnation

something is wrong. Since I am fairly new to jags so that I wouldn't exclude

the possibility of very basic error. As one can see in the bellow code the

coefficient for "a" and "b" are both wrong. Any help or advice is really

appreciated. Best, Antonio.

# jags model

model {

for (i in 1:N){

y.censor_ ~ dnorm(y.hat_,tau)T(0,)

y.hat_ <- a + b * x_

}

a ~ dnorm(0, .0001)

b ~ dnorm(0, .0001)

tau <- pow(sigma, -2)

sigma ~ dunif(0, 100)

}

# R code

## simulated data for the left tobit

N <- 1000

x <- seq(0,18, length.out = N)-8

epsilon <- rnorm(N, 20, 10)

y <- 3*x + epsilon

hist(y)

y.censor <- ifelse(y>=0,y,0)

## from survival

tobit <-survreg(Surv(y.censor, y > 0, type ='left')~ x, dist = 'gaussian')

print(tobit)

Call:

survreg(formula = Surv(y.censor, y > 0, type = "left") ~ x, dist = "gaussian")

Coefficients:

(Intercept) x

20.007908 2.870581

Scale= 10.13069

Loglik(model)= -3434.3 Loglik(intercept only)= -4004.7

Chisq= 1140.86 on 1 degrees of freedom, p= 0

n= 1000

## from MCMCpack

tobit.mcmc <-MCMCtobit(y.censor ~ x,below=0)

summary(tobit.mcmc)

Iterations = 1001:11000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 10000

- Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

(Intercept) 20.008 0.34147 0.0034147 0.003437

x 2.871 0.06465 0.0006465 0.000762

sigma2 103.088 5.00353 0.0500353 0.053112

- Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

(Intercept) 19.342 19.778 20.012 20.243 20.671

x 2.744 2.827 2.870 2.915 2.999

sigma2 93.754 99.637 102.878 106.345 113.251

tobit.inits <-

function(){list(a=rnorm(1),b=rnorm(1),sigma=runif(1),p=runif(1))}## ind <-ifelse(y.censor==0,1,0) # index for the censoring

tobit.parameters <-c("a","b","sigma")

tobit.data <- list("x","y.censor","N")

tobit.jags <- jags(tobit.data,tobit.inits,

tobit.parameters,model.file="tobit_1left.bug",

- n.chains =3,n.iter=1000,n.burnin=100,n.thin=2)

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%

|******************** ****| 100%

tobit.jags

Inference for Bugs model at "tobit_1left.bug", fit using jags,

3 chains, each with 1000 iterations (first 100 discarded), n.thin = 2

n.sims = 1350 iterations saved

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff

a 16.7 0.5 15.5 16.3 16.7 17.1 17.6 1 420

b 3.5 0.1 3.3 3.4 3.5 3.5 3.7 1 860

deviance 7011.2 2.5 7008.4 7009.4 7010.5 7012.3 7017.9 1 1400

sigma 10.7 0.3 10.1 10.5 10.7 10.9 11.3 1 1300

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 3.1 and DIC = 7014.4

DIC is an estimate of expected predictive error (lower deviance is better).