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)

    }

    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

    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

    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.

     

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks