Effect of smoking on lung cancer incidence with time-weighted pack-years

Question

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?

Answer 1

In Stan, you could start out like

data{
  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.

Answer 2

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.