Numerical problem in fitting gamma and Weibull distributions

Modeling
1

I am a beginner for RStan and I am trying to fit gamma and Weibull distributions to some failure data. When fitting the models they produce following error.

[1] “The following numerical problems occured the indicated number of times after warmup on chain 1”
count

Exception thrown at line 18: weibull_log: Shape parameter is 0, but must be > 0! 7

I get similar error messages for all the chains (4 chains) and for the gamma distributions.

Here is my attempt to fit the Weibull model.

library (rstan)

modelString = "data {
int<lower=0> J; // number of observations
real y[J]; // failure in hours
}
parameters {
real<lower=0> alpha; // shape
real<lower=0> sigma; // scale
}

model {
//define priors of alpha and beta
alpha ~ uniform(0, 20);
sigma ~ inv_gamma(0.1, 0.1);

//Likelihood of the data
y ~ weibull(alpha, sigma);
}"

mod_weib <- stan_model(model_code = modelString)

failure_data <- list (J = 88,
                  y = rep (c (8, 16, 32, 40, 56, 60, 64, 72, 80, 96, 104, 108, 112, 114, 
                              120, 128, 136, 152, 156, 160, 168, 176, 184, 194, 208, 216, 
                              224, 232, 240, 246, 256, 264, 272, 280, 288, 304, 308, 328,
                              340, 352, 358, 360, 384, 392, 400, 424, 438, 448, 464, 480, 
                              536, 552, 576, 608, 656, 716), times = c (1, 4, 2, 4, 3, 1, 1, 
                                                                        5, 4, 2, 1, 1, 2, 1, 1, 1, 1, 3, 1, 1, 5, 1, 3, 1, 2, 1, 4, 1,
                                                                        1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                                                                        1, 1, 1, 1, 1, 1, 1)) 
)

fit_weib <- sampling(mod_weib, data = failure_data, chains = 4)
fit_weib

Any suggestions to solve this problem will be highly appreciated.

Thank you in advance!

2

You’re only getting this a couple times per chain, right? And then the inference finishes? That can happen. It’s fine to ignore these errors if so.

If you really have a problem you’ll get tons and tons and tons of errors and the sampler won’t run though. Is this happening? I don’t think it should but I could always be missing something.

Your model looks right. Putting a zero avoiding prior on alpha should also help you avoid these errors. With a uniform prior like that the sampler is allowed to walk right up to zero. Maybe use a gamma of some sort?

3

Thank you Ben for the quick response. Yes, I can proceed with the posterior inferences. Anyway, I will try using a different prior.

4

Good morning,
I have an issue with ffiting a Weibull with lognormal priors on the parameters.

here is my code


model <- 'data {
          int<lower=0> N;// Number of observations
          vector[N] time; //predictor
          vector[N] concentration;  //response
          
          real lambda_pop;
          real beta_pop;
          real<lower=0> omega_lambda;
          real<lower=0> omega_beta;
        }
        parameters {
          vector<lower=0>[2] param;
        }

        model {
          //Priors
          param[1] ~ lognormal( lambda_pop , omega_lambda);
          param[2] ~ lognormal( beta_pop , omega_beta);
          concentration ~ weibull(param[1], param[2]);
        }'



modelstan <- stan_model(model_name = "rtte",model_code = model)

stan_data <- list(N = length(obs),concentration = obs
					,time = time,
					lambda_pop=lambda,beta_pop=beta,
					omega_lambda=omega.eta[1,1],omega_beta=omega.eta[2,2])

	warmup <- 1000
	fit <- sampling(modelstan, data = stan_data, iter = 100000+warmup,
		warmup = warmup,chains = 1,algorithm = "NUTS") #can try "HMC", "Fixed_param"

and i get this error

Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
error occurred during calling the sampler; sampling not done

My values for lambda_pop and beta_pop are positive so I don’t understand the issue.
Thanks for the help.

Belhal

5

Actually since my vector of observations include some 0 and 1, I believe I need ```weibull_ccdf```` to deal with the zeros (it is a right censored weibull model).

Another question would be: my hazard function is equal to h(t,param) = beta/lambda*(t/lambda)^{beta-1}
How would you proceed to include beta, lambda and time t into your abstraction weibull(.)?

Thanks a lot

6

I’ve been encountering similar problems in my Weibull models. At least for my data (which is simulated to have shape \approx 1), the Weibull is finicky about initial values for the shape parameter. The uniform(-2, 2) for log(shape) (log(param[1] for you) is just too dispersed. Does it work better if you run with init = "0"? If so, you may want to change the initialization range.

I’m not sure I understand your second question about including beta, lambda, and t. Every Weibull pdf parameterization will have a corresponding hazard parameterization.

(Also, are you saying that some elements of the concentration vector will have a zero value? How does a zero value correspond to right censoring?)

7

Thanks for your reply.
I am predicting a binary variable. 1 if the event occurred and 0 if the right censoring time is reached (and 0 at the initial time 0).
So my vector is actually composed of a 0, many 1’s and a final 0 (giving the information that the censoring time is reached).

Also, regarding my particular hazard function of parameter beta and lambda, and regressor time t, how would you define the model using the weibull function implemented in stan?

Thanks

8

Are you saying concentration is 0/1? That’s not what it looks like from your code. If X is Weibull-distributed, its support should be \mathbb{R}^+, and a right-censored version would have support (0, U) for some upper limit U > 0. The typical approach to right-censoring has a variable T_i=\min(X_i, C_i) that contains the value of X_i (if observed) or the right censoring value C_i. There is also a separate binary variable \delta_i that indicates whether the observation was censored or not. If you have censoring, you cannot just write concentration ~ weibull(param[1], param[2]).

It would be helpful if you explained your input data more. Are these event times? Chemical concentrations with an upper limit of quantification?

9

You are right. Here is my problem: I am considering a Weibull model for right-censored repeated time to event data. My observations are thus the Times at which those repeated events occurred.
My hazard function is defined as h(t, \psi_i) = \frac{\beta_i}{\lambda_i}\left(\frac{t}{\lambda_i}\right)^{\beta_i-1}

I put lognormal priors on the parameters \lambda and \beta.

I would like to use rstan to sample from the conditional distribution of the latent variables (\lambda and \beta) given the observations (that are the times of those events).

How would you proceed?

10

Ah, okay. This is a lot more like a classical use of the Weibull than I thought! I am somewhat confused why your hazard function has \beta and \lambda with i subscripts. Your code has only two parameters, not 2N. (It is probably best to change this after you get a simpler model working.)

Another question is how you are dealing with repeated events within an observation unit. Let’s say observation unit i has two events: one at t_1 and one at t_2. With a Weibull model, the hazard changes over time as a function of the shape parameter. You can choose to model gap times (t_1 - 0) and (t_2 - t_1) if you think that the hazard “resets” after each event occurrence. Alternatively, you can model the actual event times t_1 and t_2, accounting for the fact that the at-risk period for the second event did not start until t_1. Which of these approaches describes your scenario best?

At this point I suggest making a new post about recurrent event Weibull models so that others can find this topic in the future.

11

You are right the presence of the subscript i is a typo.

Thanks a lot for your help. I will indeed create a new topic on design on weibull model on rstan with right censored repeated time to event data.

Thanks again!
Belhal