Many divergent transitions with truncated data model

[Eric Mittman]

I am building a model for lifetime data that assumes a proportion, p, of units are defective and subject to an additional failure mode. The distribution of X is

Pr(X <= t) = F(t) = 1 - (1 - p * F1(t)) * ( 1 - F2(t))

where X is the failure time for the unit. I assume F1 and F2 are Weibull distributions, and I actually log the data. Y ~ Weibull(shape, scale) implies log(Y) ~ gumbel(-mu, sigma), where mu = log(scale) and sigma = 1/shape.

I simulated data from the model and left-truncated the data at randomly chosen time points. I sampled 4000 iterations from each of 8 chains (with 2000 warmup) setting init=0 ( after running into some problems with init="random".) The samples I obtained look reasonable (pairs_stan.png, attached), but these divergent transitions are concerning. This problem was not evident in earlier fits that did not involve truncation.

The warnings I see suggest that I might constrain mu, but it is ordered type for identifiability. I am hoping someone might suggest improvements to the Stan code (mcr_trunc_obs.stan, attached).

Here is the call to sampling:

s <- sampling(m2, data = data, init = 0, iter = 2000, chains=8, thin=2, control = list(stepsize = .1, adapt_delta=.95))

Here is print(s):

Inference for Stan model: mcr_trunc_obs.

8 chains, each with iter=2000; warmup=1000; thin=2; 

post-warmup draws per chain=500, total post-warmup draws=4000.

                 mean se_mean   sd     2.5%      25%      50%      75%

mu[1]            7.12    0.00 0.05     7.01     7.08     7.12     7.15

mu[2]            9.99    0.00 0.02     9.96     9.98     9.99    10.01

log_sigma[1]    -0.52    0.00 0.04    -0.57    -0.55    -0.53    -0.50

log_sigma[2]    -0.05    0.00 0.02    -0.09    -0.07    -0.05    -0.04

p                0.12    0.00 0.01     0.10     0.12     0.12     0.13

lp__         -7423.71    0.04 1.62 -7427.70 -7424.59 -7423.38 -7422.51

                97.5% n_eff Rhat

mu[1]            7.22   861 1.00

mu[2]           10.03  1700 1.00

log_sigma[1]    -0.42   721 1.02

log_sigma[2]    -0.02  1231 1.00

p                0.14  1117 1.00

lp__         -7421.48  1349 1.01

The true values are: mu = c(7, 10), log_sigma = c(-.69, 0), p = .1

The warnings I see:

The following numerical problems occured the indicated number of times after warmup on chain 8

                                                                                                     count

Exception thrown at line 40: Exception thrown at line 6: gumbel_log: Location parameter is -inf, but     9

When a numerical problem occurs, the Hamiltonian proposal gets rejected.

If the number in the 'count' column is small, do not ask about this message on stan-users.

Warning messages:

1: There were 2787 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help.

2: Examine the pairs() plot to diagnose sampling problems

SessionInfo:

R version 3.3.1 (2016-06-21)

Platform: x86_64-w64-mingw32/x64 (64-bit)

Running under: Windows >= 8 x64 (build 9200)

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] rstan_2.11.1       StanHeaders_2.11.0 ggplot2_2.1.0     

loaded via a namespace (and not attached):

 [1] colorspace_1.2-6   scales_0.4.0       plyr_1.8.4         tools_3.3.1       

 [5] inline_0.3.14      gtable_0.2.0       gridExtra_2.2.1    Rcpp_0.12.6       

 [9] KernSmooth_2.23-15 grid_3.3.1         munsell_0.4.3      stats4_3.3.1 

[Bob Carpenter]

You don't need to assign before returning:

real sev_logpdf(real y, real mu, real sigma){
real out;
out = gumbel_lpdf(-y| -mu, sigma);
return out;
}

That can just be return gumbel_lpdf(...).

This prior for p is saying that you expect the
value to be really close to 1 most of the time, but
sometimes very close to 0 and almost never in between:

p ~ beta(.1, .9);

I'd be inclined to go with a more regularizing prior.
It can really help with computation.

It looks like the divergences (I think that's what's
printed in red?) may all be due to numerical issues in
evaluation.

For efficiency, you don't want to keep taking those logs
and exps of the same values inside the loops --- compute them
on the outside and store in a local variable.

You might want to try coding up the Gumbel directly with
a log scale parameter so you don't need the log to exp to log
transform, which can be unstable. And then reducing to the point
where you can

And I'd recommend dropping the interval constraints on
log_sigma and just using the prior as is --- normal(0, 1)
should keep it in (-5, 5) unless there's a ton of data
pulling it out, in which case you probably don't want that
constraint anyway.

Also, there's a log1m_exp(u) which is equal to
log_diff_exp(0, u).

- Bob



> On Aug 26, 2016, at 6:30 PM, Eric Mittman <mathteach...@gmail.com> wrote:
>
> I am using RStan version

>
> I am building a model for lifetime data that assumes a proportion, p, of units are defective and subject to an additional failure mode. The distribution of X is
>
> Pr(X <= t) = F(t) = 1 - (1 - p * F1(t)) * ( 1 - F2(t))
>

> where X is the failure time for the unit. I assume F1 and F2 are Weibull distributions, and I actually log the data. This implies that log(X) ~ gumbel(-mu, sigma), where mu = log(scale) and sigma = 1/shape of the Weibull.

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <mcr_trunc_obs.stan><pairs_stan.png>

[Eric Mittman]

Bob,

Thanks for the tips. I no longer get divergent transitions and now my code is cleaner. Revised code attached.

Follow up question: Suppose I wanted to put a muti_normal prior on (mu[1], log_sigma[1]). Can I do this:

parameters{

    ordered[2] mu;

    vector[2] log_sigma;

    ...

}

transformed parameters{

    vector[2] theta;

    theta[1] = mu[1];

    theta[2] = log_sigma[1];

}

model{

    ...

    theta ~ multi_normal(m, C);

}

?

Thanks for your time,

Eric

[Bob Carpenter]

If you want to put a multivariate prior on theta, that'll
work. You can transform log_sigma to sigma via exp() and
use them as parameters.

We recommend the LKJ prior on a correlation matrix (can concentrate
prior around unit correlation matrix or make it uniform) combined
with an independent prior on the scales. There's a discussion in
the regression chapter of the manual.

- Bob

> On Aug 26, 2016, at 10:46 PM, Eric Mittman <mathteach...@gmail.com> wrote:
>
> Bob,
>
> No more divergent transitions and cleaner code. Revised code attached.
>
> Follow up question: Suppose I wanted to put a muti_normal prior on (mu[1], log_sigma[1]). Do you think I can get away with constructing such a vector in the transformed parameters block, i.e.
> tranformed parameters{

> vector[2] theta;
> theta[1] = mu[1];
> theta[2] = log_sigma[1];
> }
>

> Thanks for your time,
>
> Eric
>