## tobit model

Help
Anonymous
2011-03-07
2012-09-01
• Anonymous - 2011-03-07

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)

}

# 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

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

1. 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).

• Martyn Plummer - 2011-03-14

You are confusing censoring and truncation. Please have a look at the JAGS
user manual.

To prepare your data for JAGS in R, use this

```y.censored <- ifelse(y>=0, y, NA)
y.ind <- as.numeric(y>=0)
```

ao that censored observations are set to missing, not zero. Then the JAGS
model is:

```model {
for (j in 1:N){
y.ind[j] ~ dinterval(y.censored[j], 0)
y.censored[j] ~ dnorm(y.hat[j],tau)
y.hat[j] <- a + b * x[j]
}
a ~ dnorm(0, .0001)
b ~ dnorm(0, .0001)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)
}
```

We need to explicitly set initial values for the missing y.censored values so
that they are consistent with the a posteriori constraint imposed by y.ind.
It is easiest to generate starting values with a=0, b=0.

```tobit.inits <- function() {
sigma <- runif(1)
list(a=0, b=0, sigma=sigma, p=runif(1),
y.censored=ifelse(y.ind, NA, -abs(rnorm(N, 0, sigma))))
}
```

Finally you can run this through JAGS using the rjags package

```m <- jags.model("tobit.bug",
data=list(y.censored=y.censored, y.ind=y.ind, x=x, N=N),
inits=tobit.inits)
m.out <- coda.samples(m, c("a", "b", "sigma"), n.iter=10000)
```

• Anonymous - 2011-03-23

Thank you very much for your reply. I guess the main issue for me is that I
still don't fully understand the logic of dinterval(,). Now I am trying to
work with double censoring but it is unclear how to specify dinterval(,) for
those situations as I am getting errors all over the place.