How do define a hazard function in rstan

Hi everyone,

I am trying to write a rstan code for parameter calibration. In my calibration model, there is a hazard function h(t) which uses to derive the age t. This model comes from the paper MILC: A microsimulation model of the natural history of lung cancer | International Journal of Microsimulation.
This is how I tried to define the function

  real h_integral(real t, real X, real mu, real nu, real alpha1, real alpha2, real gamma_ns, real alpha_ns, real q, real[] x_r, int[] x_i){
    real alpha = alpha_ns * (1 + alpha1* pow(q,alpha2));
    real gamma = gamma_ns * (1 + alpha1* pow(q,alpha2)); 
    real B = (-gamma + sqrt(gamma^2 + 4 * alpha * mu))/2;
    return (nu * mu * X * (exp((gamma+2*B)*t)-1))/(gamma+B*exp((gamma+2*B)*t)+1);

Can someone please guide me on a way of defining this function in order to return age t ?

(In R, I defined it as

hr[[i]] <- function(t) (nu * mu * X * (exp((gamma[i]+2*B[i])*t)-1))/(gamma[i]+B[i]*exp((gamma[i]+2*B[i])*t)+1) # Hazard rate[[i]] <- function(t) 1-exp(-(integrate(hr[[i]], lower = 0, upper = t)$value)) # CDF function)

Thank you

@jonah @Bob_Carpenter Please have a look at this. Thank you all.

I want to make sure I understand your question, because there are a few things that I find confusing here:

  • you say that you want the function to return t, but you provide t as an input to the function.
  • you already have what looks like a syntactically valid Stan function, but it has a bunch of unused inputs. Does this function work, or does it error?
  • what is going on with the inputs x_r and x_i`, which never show up anywhere in the Stan or R code except as inputs?

My guess is that you want to modify your function so that its new output will be one minus the exponential of the definite integral from t to 0 of your current function.

The thing is, your R function hr[[i]] has an antiderivative that according to Wolfram is readily expressed in the Stan language. So if I’ve understood correctly what your question is, you should skip the numerical integration and just plug in the exact function that you seek.

Does that help?