Area above threshold feature

General
1

Let M_{i}(t) be some summary of risk factors at time t for participant i from the longitudinal risk factor model. The hazard sub-model is
h_{i}(t \vert M_{i}(t)) = h_{0}\exp\{\alpha^{'T}w_{i} + \alpha \int_{s=0}^{t} m_{i}(s) I(m_{i}(s) > \gamma)ds
where \gamma is a risk factor threshold and linear predictor of longitudinal sub-model
m_{i}(t)=x_{i}^T(t)\beta+z_{i}^T(t)b_{i}.

Is there any way to specify the area above threshold feature (integration part in hazard model) like this in STAN?

2

Anyone can help me out, please? @maxbiostat @Bob_Carpenter

3

Hi,

This thread might be of value. Just beware that the behaviour of integrate_1d might have changed in recent releases.

1 Like
4

Thank you for your reply! That thread is helpful but I am wondering how I can write the indicator part involving longitudinal parameters in transformed parameters block.

5

Right. Here’s some stan code:

functions{
  real indicator(real x, real c){
    real ans = x > c  ? 1 : 0; 
    return(ans);
  } 
  real phi( real x, real x_c,
  array[] real theta, array[] real x_r, array[] int x_i ) {  
    real ans = x^2 * exp(-x^2/2) * indicator(x, 0);
    return(ans);
    
  }
}
transformed data {
  real x_c[0];
  real x_r[0];
  int x_i[0];
  real theta[2];
}
generated quantities{
  real the_int = integrate_1d (phi,
  negative_infinity(), positive_infinity(), theta, x_r, x_i);
}

I ran with the following R code:

library(cmdstanr)

integrator <- cmdstan_model("indicator_integration.stan")

integrator$sample(fixed_param = TRUE, chains = 1,
                  iter_warmup = 0, iter_sampling = 1)

phi <- function(x) x^2 * exp(-x^2/2) * as.numeric(x >= 0)
integrate(phi, 0, Inf)
sqrt(pi/2)
1 Like
6

Thank you so much! Now, I can write the area above threshold feature in stan. Just one question to confirm that I am doing it correctly. Is it okay to compute this feature (involves indicator function) within transformed parameter block and add it as a covariate in the survival model?

7

I can’t see why not, but then again I’m not an expert on survival models. Maybe if @harrelfe sees this he might be able to advise. I suggest you ask a new, separate question.