I'm trying to specify a model for predicting lung cancer incidence based on a person's smoking history. Most studies simply use pack-years, a measure of the total cigarettes smoked in one's lifetime. I'd like to incorporate time weights into the calculation of pack-years, where the weights are estimated as part of the regression equation. What this adds is the ability to say, for example, that the effect of smoking a cigarette today is 3 times more influential than a cigarette smoked 10 years ago. I'm thinking of a Cox proportional hazards model (semi-parametric) like the one below, or maybe an accelerated failure time (parametric) model that might be more appropriate for prediction. How can I specify this model, preferably using an R package?

$\lambda(t|X_i) = \lambda_0(t)\exp(\beta_1{Age}_i + \beta_2{Gender}_i + \beta_3 WPY_{i})$

The hazard rate at time $t$ is a function of age, gender and the weighted pack-year covariate.

$WPY_i ~ \sim \sum_{j=1}^{m} \omega_j \space PY_{ij}$

Weighted pack-years are calculated as, summing over each previous month $m$, the number of pack-years smoked in that month $PY_{ij}$ multiplied by the weight for that month $\omega_j$. Note that month $m$ means "months since now", not the calendar month. I expect that the weights $w_j$ follow something like a Gamma distribution and collectively sum to 1.

I attempted this already by wrapping a Cox model inside of optim() to get values of $\omega_j$ that maximize the likelihood of the observed data. However, I'm not sure that I'm setting that up correctly and there's probably a less hacky way. I'm thinking this is possible with Stan. Any ideas?


2 Answers 2


In Stan, you could start out like

  int<lower=1> N;
  int<lower=1> M;
  matrix<lower=0>[N, M] PY;
  vector<lower=0>[N] age;
  vector<lower=0,upper=1> gender;
parameters {
  simplex[M] omega;
  vector[3] beta;
model {
  vector[N] WPY = PY * omega;
  vector[N] eta = beta[1] * age + beta[2] * gender + beta[3] * WPY; 

But how exactly you specify the rest of the model block is up to you. Technically, it would be insufficient in Stan to allow $\lambda_0\left(t\right)$ to be any function (without some prior over all functions) because that would yield an improper posterior distribution. However, $\lambda_0\left(t\right)$ could be any smooth function, in which case you could use a Gaussian process prior or a spline or something like that. Just Google for Bayesian Cox model to find various alternatives. There are plenty of examples of parametric survival models in SurvivalStan.

  • $\begingroup$ Good point about lambda_0. This is a good start, although I'm pretty new to Stan. Could you give another example of, rather than having a separate coefficient for each month, have them follow a Weibull distribution (only 2 parameters)? $\endgroup$
    – rsoren
    Jan 8, 2018 at 17:42
  • $\begingroup$ In that case, change simplex[M] omega; to vector<lower=0>[M] omega; and add target += weibull_lpdf(omega | shape, scale); to the model block, where shape and scale would typically be passed in through the data block. $\endgroup$ Jan 9, 2018 at 17:56

Insert WPY from your second equation into the first, multiply beta into the sum, and note that you now have a standard cox ph regression where the coefficients on the PY are each equal to beta3omegaj.

So in practice, just do a cox proportional hazards regression including all the PY as regressors. If there is actually a decaying effect, you will see that effect in the coefficients. You can also test for the decaying effect by testing equality of coefficients.

I don’t think you should require the omegas to sum to one, because the coefficients are scaled by beta3 anyway. (I.e. if you were to put in a restraint forcing the weights to sum to one, that would get absorbed by the estimate on beta3).

Update: R code for estimating the model

model1 <- coxph(Surv(time, cancer) ~ age + gender + PY1 + PY2 + PY3 +  PY4 + PY5, data =  mydata)

PY1 would be the persons pack-years of smoking in the month before the event time, PY2 would be two months ago etc up to the number of months that are relevant (I think you would probably want to estimate per year or quarter when you exceed some small number of months, to reduce the dimensionality of the mode).

If your data are longitudinal (i.e. one row per person per month) you need to aggregate them up so you have one row per person with their longest survival time and one variable/column per month of smoking.

  • $\begingroup$ Could you give an example of how to do this with coxph() from the survival package in R? $\endgroup$
    – rsoren
    Jan 8, 2018 at 17:40
  • $\begingroup$ I'll update the answer with a couple of lines of syntax, but really the point is that the model is a perfectly standard multivariate cox ph regression, so you can use all the usual methods for analysing, estimating etc. It's only the input variables that are different from what is usual. $\endgroup$
    – BjaRule
    Jan 11, 2018 at 7:23
  • $\begingroup$ I'm pretty sure this would suffer from collinearity issues. That's why I was looking for a multilevel model. Thanks though $\endgroup$
    – rsoren
    Jan 14, 2018 at 1:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.