beginner requests help: censored regression

[vit...@osu.edu]

For my first stan analysis, I wanted to extend the
censored regression model decribed in the stan_2.9.0
manual, pp. 108-109, to allow for covariates, ie to
implement what econometricians call the tobit model.
This is the model version that does not impute the
censored values.

My stan results are so different from the ML results
using the tobit estimator from the AER package that
I think I must be doing something wrong. (See
particularly the mean of the distribution of beta[2]).
With stan using by default uninformative priors,
one would expect the results to be close, no?

Could I  persuade someone to take a look at this
and see if they can diagnose what's going on?

Thanks!!

##### R session     (results below)
library(rstan)
library(AER)

load("c:/work/limdep/aasho/r/aasho.rdata")

U<-max(rigid25$ltrig)

model<-ltrig~lslabdd+laxles+lwta

# fit the tobit model
(summary(tobit(model,data=rigid25,right=U)))

# now set up the data for stan
id_cens<-rigid25$ltrig==U 
mmx<-model.matrix(model,data=rigid25)
y_obs<-rigid25$ltrig[!id_cens]
x_obs<-mmx[!id_cens,]
x_cens<-mmx[id_cens,]
K<-dim(mmx)[2]
N_obs<-length(y_obs)
N_cens<-sum(id_cens)

# lets recheck tobit using these constructed data
dd<-cbind(y_obs,x_obs)
dd<-rbind(dd,cbind(U,x_cens))
dd<-dd[,c(1,3,4,5)] # remove the intercept column
dd<-data.frame(dd) # convert to data frame
colnames(dd)[1]<-"ltrig" # name of dependent variable

(summary(tobit(model, data=dd,right=U)))
# this agrees with the earlier version

# now for stan: form the data list
data<-list(y_obs=y_obs, x_cens=x_cens, x_obs=x_obs,
                     K=K, N_obs=N_obs, N_cens=N_cens, U=U)

# re-running with iter=10000, warmup=2000
# doesn't change the restults signifcantly                     
fit<-stan(file="aasho2.stan",data=data,
                            iter = 5000 , warmup=1000 , chains = 1 )
print(fit)
get_stanmodel(fit)



##### results

# tobit results

> (summary(tobit(model,data=rigid25,right=U)))

Call:
tobit(formula = ltrig ~ lslabdd + laxles + lwta, right = U, data = rigid25)

Observations:
         Total  Left-censored     Uncensored Right-censored
           264              0             73            191

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) 13.50466    0.30659  44.049   <2e-16 ***
lslabdd      5.03874    0.32869  15.330   <2e-16 ***
laxles       2.26879    0.24178   9.384   <2e-16 ***
lwta        -3.23973    0.25971 -12.475   <2e-16 ***
Log(scale)  -1.00280    0.08374 -11.975   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Scale: 0.3668

Gaussian distribution
Number of Newton-Raphson Iterations: 8
Log-likelihood: -55.1 on 5 Df
Wald-statistic: 258.3 on 3 Df, p-value: < 2.22e-16


### repeat of tobit model (here I just show the coefficient block)
## the point is, they're the same as before
> (summary(tobit(model, data=dd,right=U)))
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) 13.50466    0.30659  44.049   <2e-16 ***
lslabdd      5.03874    0.32869  15.330   <2e-16 ***
laxles       2.26879    0.24178   9.384   <2e-16 ***
lwta        -3.23973    0.25971 -12.475   <2e-16 ***
Log(scale)  -1.00280    0.08374 -11.975   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Scale: 0.3668


## now stan

> fit<-stan(file="aasho2.stan",data=data,
+                             iter = 5000 , warmup=1000 , chains = 1 )

SAMPLING FOR MODEL 'aasho2' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)#
#  Elapsed Time: 41.016 seconds (Warm-up)
#                125.062 seconds (Sampling)
#                166.078 seconds (Total)
#
> print(fit)
Inference for Stan model: aasho2.
1 chains, each with iter=5000; warmup=1000; thin=1;
post-warmup draws per chain=4000, total post-warmup draws=4000.

           mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
beta[1]   16.32    0.02 0.48   15.48   15.97   16.28   16.61   17.34   815    1
beta[2]   13.43    0.05 1.37   11.03   12.48   13.37   14.31   16.33   627    1
beta[3]    5.62    0.03 0.64    4.48    5.17    5.59    6.03    7.00   654    1
beta[4]   -8.70    0.04 0.91  -10.66   -9.29   -8.66   -8.05   -7.12   607    1
sigma      1.19    0.00 0.13    0.97    1.10    1.18    1.27    1.47   667    1
lp__    -238.28    0.06 1.63 -242.23 -239.14 -237.93 -237.08 -236.13   805    1

Samples were drawn using NUTS(diag_e) at Wed Apr 20 09:26:54 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

> get_stanmodel(fit)
S4 class stanmodel 'aasho2' coded as follows:
# this allows for mu to be determined by data    as x * beta
# this version does not impute the censored values


data {
  int<lower=0> N_obs ;
  int<lower=0> N_cens ;
  int<lower=0> K ; // number of predictors
  real y_obs[N_obs]    ;
  real<lower=max(y_obs)> U;
  matrix[N_obs,K] x_obs; // uncensored predictors data mx, with 1s in col 1
  matrix[N_cens,K] x_cens; // censored predictors data mx
}

parameters {
  vector[K] beta ;
  real<lower=0> sigma ;
}

model {
  y_obs ~ normal(x_obs * beta,sigma);
  increment_log_prob(N_cens * normal_ccdf_log( U , x_cens * beta , sigma)) ;
}



### session info

> sessionInfo()
R version 3.2.4 (2016-03-10)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows XP (build 2600) Service Pack 3

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252  
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                         
[5] LC_TIME=English_United States.1252   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base    

other attached packages:
[1] AER_1.2-4       survival_2.38-3 sandwich_2.3-4  lmtest_0.9-34   zoo_1.7-12    
[6] car_2.1-1       rstan_2.9.0-3   ggplot2_2.1.0 

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.3        Formula_1.2-1      splines_3.2.4      MASS_7.3-45      
 [5] munsell_0.4.3      colorspace_1.2-6   lattice_0.20-33    minqa_1.2.4      
 [9] plyr_1.8.3         tools_3.2.4        nnet_7.3-12        pbkrtest_0.4-6   
[13] parallel_3.2.4     grid_3.2.4         gtable_0.2.0       nlme_3.1-126     
[17] mgcv_1.8-12        quantreg_5.21      MatrixModels_0.4-1 lme4_1.1-11      
[21] Matrix_1.2-4       gridExtra_2.2.1    nloptr_1.0.4       codetools_0.2-14 
[25] inline_0.3.14      scales_0.4.0       stats4_3.2.4       SparseM_1.7

[Jonah Gabry]

The Stan code from the manual at the bottom of pg. 108 (and top of 109) is just for right censoring so it wouldn't be equivalent to your use of the tobit function, which I think defaults to left censoring at 0. That is, I'm pretty sure the default in AER::tobit is left = 0 so your formula

tobit(model,data=rigid25,right=U)

is really equivalent to 

tobit(model,data=rigid25,right=U, left=0)

If you only want right censoring then I think you need to set left = -Inf.

Hope that helps, 

Jonah

[Viton, Philip]

Thanks you very much!  You’re correct that AER:tobit defaults to left=0, but in my case the actual data meant that specifying left=0 wasn’t needed. As it happens, the tobit results in my note agree with a result in the published literature, so the tobit results aren’t in question.

 

So my real question is: what’s wrong with my stan implementation such that it doesn’t agree with the tobit results? I’m sorry if that wasn’t clear in my original posting.

 

 

 

 

From: stan-...@googlegroups.com [mailto:stan-...@googlegroups.com] On Behalf Of Jonah Gabry
Sent: Wednesday, April 20, 2016 2:57 PM
To: Stan users mailing list <stan-...@googlegroups.com>
Subject: [stan-users] Re: beginner requests help: censored regression

 

The Stan code from the manual at the bottom of pg. 108 (and top of 109) is just for right censoring so it wouldn't be equivalent to your use of the tobit function, which I think defaults to left censoring at 0. That is, I'm pretty sure the default in AER::tobit is left = 0 so your formula

 

tobit(model,data=rigid25,right=U))

is really equivalent to 

 

tobit(model,data=rigid25,right=U, left=0)

If you only want right censoring then I think you need to set left = -Inf.

 

Hope that helps, 

 

Jonah

 

 

 

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/1tnXA0WiSx0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jonah Sol Gabry]

I'll take a look at your Stan code, but in case it's not immediate apparent what the problem is, is the data you're using available somewhere (or small enough to share) or could you reproduce the issue with simulated data? It's typically easier to diagnose issues when we can run the model.

Jonah

[vit...@osu.edu]

Thanks!  I'm attaching the data file (aasho.rdata) in case the stand code's problem prove hard to unravel. Thanks again for all your trouble!

Regards,

On Wednesday, April 20, 2016 at 3:40:56 PM UTC-4, Jonah Gabry wrote:

I'll take a look at your Stan code, but in case it's not immediate apparent what the problem is, is the data you're using available somewhere (or small enough to share) or could you reproduce the issue with simulated data? It's typically easier to diagnose issues when we can run the model.

Jonah

On Wednesday, April 20, 2016, Viton, Philip <vit...@osu.edu> wrote:

Thanks you very much!  You’re correct that AER:tobit defaults to left=0, but in my case the actual data meant that specifying left=0 wasn’t needed. As it happens, the tobit results in my note agree with a result in the published literature, so the tobit results aren’t in question.

 

So my real question is: what’s wrong with my stan implementation such that it doesn’t agree with the tobit results? I’m sorry if that wasn’t clear in my original posting.

 

 

 

 

From: stan-...@googlegroups.com [mailto:stan-users@googlegroups.com] On Behalf Of Jonah Gabry
Sent: Wednesday, April 20, 2016 2:57 PM
To: Stan users mailing list <stan-...@googlegroups.com>
Subject: [stan-users] Re: beginner requests help: censored regression

 

The Stan code from the manual at the bottom of pg. 108 (and top of 109) is just for right censoring so it wouldn't be equivalent to your use of the tobit function, which I think defaults to left censoring at 0. That is, I'm pretty sure the default in AER::tobit is left = 0 so your formula

 

tobit(model,data=rigid25,right=U))

is really equivalent to 

 

tobit(model,data=rigid25,right=U, left=0)

If you only want right censoring then I think you need to set left = -Inf.

 

Hope that helps, 

 

Jonah

 

 

 

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/1tnXA0WiSx0/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/1tnXA0WiSx0/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Jonah Gabry]

Ok I've figured out the issue. It's one of those annoying ones that's so simple but takes a while to actually spot. You just need to get rid of the multiplication by N_cens in the increment_log_prob line. Here's why:

In the example in the Stan manual mu is a scalar parameter so, since there are N_cens censored observations, in the line 

increment_log_prob(N_cens * normal_ccdf_log(U, mu, sigma))

you need the multiplication by N_cens. However, when you go from mu to x_cens * beta you have a vector instead of a scalar and the call to normal_ccdf_log becomes vectorized. The vectorized call to normal_ccdf_log obviates the need to multiply by N_cens since x_cens * beta contains a value for each of the N_cens censored observations. So you just want 

increment_log_prob(normal_ccdf_log(U, x_cens * beta, sigma)) ;

Anyway, when you make that change you should get the correct estimates. 

On another note, I would run at least 4 chains. You also shouldn't need more than the default number of iterations per chain (2000) for a model like this, and you probably don't even need close to that many (I would double check that though).

Best,

Jonah

[vit...@osu.edu]

Thanks very much! (And the explanation wqs a bit of an embarrassment for me - I really should have noticed).  Your advice as to chains is well taken, though at this stage, given that this was my first attempt, I wanted to be confident that the model was right before working on production runs.

[Jonah Sol Gabry]

On Thursday, April 21, 2016, <vit...@osu.edu> wrote:

Thanks very much! (And the explanation wqs a bit of an embarrassment for me - I really should have noticed). 

It's one of those things that seems obvious in hindsight but is easy to overlook. Took me a while to see it too. Also, vectorization is great but in this case we would have probably noticed the problem very quickly if normal_ccdf_log didn't accept a vector and threw an error. Anyway, glad we got it sorted out.

Jonah

[Bob Carpenter]

The reason you want to run multiple chains is to make
sure that diffuse initializations converge to the same
region of the posterior. So it's important to do it
before you do a final run.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.