Cohort Regression Models vs Individual Regression Models

Question

I have always had the following question: Were regression models originally designed and intended to be used for predictions/inference only at the "cohort level" (i.e. for large groups of people), or were they also intended for prediction/inference at the individual level (provided the data is large enough and "well behaved" i.e. moderate levels of variance)?

To illustrate my example, I use regression models from the domain of "survival analysis". As a refresher, survival analysis (also called "time to event" modelling) is interested in modelling the "time at which some "event" happens" (e.g. medical death, failure, mortgage defaulting, etc.), with specific focus to accommodate "incomplete" data (i.e. censored data - but this is not overly important for my question). Historically, survival analysis models were useful in estimating the "survival rates" in clinical studies, and estimating the effect of different covariates (i.e. predictor variables) on the survival rate. Using the R programming language, I will demonstrate applications of survival regression models at the cohort level and at the individual level.

1) Survival Analysis Regression at the Cohort Level

Historically, survival regression models were mainly used for analysis at the cohort level - for example, comparing the survival rates between two different groups in a medical study (e.g. survival rate for men who smoke cigarettes vs women who smoke cigarettes).

In the below example, survival analysis is performed to estimate the overall survival rate for all army veterans, and then to compare the survival rate for veterans who received "Treatment A" vs. "Treatment B" in this case, the "treatment" is the cohort, i.e. there are 2 cohorts):

#load libraries
library(survival)
library(ranger)
library(ggplot2)
library(dplyr)
library(ggfortify)

#load data
data(veteran)

#analysis on all veterans 
km <- with(veteran, Surv(time, status))

#cohort level analysis
km_trt_fit <- survfit(Surv(time, status) ~ trt, data=veteran)


#visualize results
autoplot(km_fit)
autoplot(km_trt_fit)

We can also fit a survival regression model (the "Cox Proportional Hazards Model") to estimate the effect of individual variables on the survival rate:

#fit model
cox <- coxph(Surv(time, status) ~ trt + celltype + karno                   + diagtime + age + prior , data = vet)

#view results : summary of variables on survival rate
summary(cox)

 Call:
 coxph(formula = Surv(time, status) ~ trt + celltype + karno + 
     diagtime + age + prior, data = vet)
 
   n= 137, number of events= 128 
 
                         coef  exp(coef)   se(coef)      z Pr(>|z|)    
 trttest            2.946e-01  1.343e+00  2.075e-01  1.419  0.15577    
 celltypesmallcell  8.616e-01  2.367e+00  2.753e-01  3.130  0.00175 ** 
 celltypeadeno      1.196e+00  3.307e+00  3.009e-01  3.975 7.05e-05 ***
 celltypelarge      4.013e-01  1.494e+00  2.827e-01  1.420  0.15574    
 karno             -3.282e-02  9.677e-01  5.508e-03 -5.958 2.55e-09 ***
 diagtime           8.132e-05  1.000e+00  9.136e-03  0.009  0.99290    
 age               -8.706e-03  9.913e-01  9.300e-03 -0.936  0.34920    
 priorYes           7.159e-02  1.074e+00  2.323e-01  0.308  0.75794    

#model fit statistics (good sign: p-value is less than alpha at 0.05)
 Concordance= 0.736  (se = 0.03 )
 Rsquare= 0.364   (max possible= 0.999 )
 Likelihood ratio test= 62.1  on 8 df,   p=2e-10
 Wald test            = 62.37  on 8 df,   p=2e-10
 Score (logrank) test = 66.74  on 8 df,   p=2e-11 

We can also visualize this cohort level survival regression model:

#plot 
autoplot(cox_fit)

These are some of the main applications of survival regression models at the cohort level.

2) Survival Regression Analysis at the Individual Level

I have heard that traditionally, survival analysis regression was not extended to the individual level and was only meant for the cohort level. This is because when these models were developed (e.g. 1950's- 1970's), datasets were very small - thus, survival models applied to small data and the individual level were almost destined to fail, because they data (when analyzed at the individual level) contained insufficient information to generalize to new data (in statistics, this is called the "ecological fallacy" and the "reference class problem"). As a quick example, imagine you have a dataset with 50 "healthy patients" and 50 "sick patients" with the goal of studying survival rates. The "sick patients" are made up of 25 patients with diabetes, 22 patients with kidney failure and 3 patients with liver failure. Assuming that "healthy patients" on average live longer than "sick patients" - if you wanted to make a regression model for only patients with liver failure, you would be forced to make a model with only 3 data points: apart from having too few observations for a model, these 3 patients might contain biases/outliers and likely are not representative of the true population of patients with liver failure, and thus this model is likely to perform poorly in the real world on new patients. However, if you decide to aggregate all the "sick patients" together - you now have a lot more data points and you can expose your model to a larger variety of data, making your model likely to make more accurate predictions about the "general" survival rates of "sick patients".

However, in recent years - researchers now have access to far larger datasets (e.g. 100,000 rows) compared to the data available when the original survival models were developed (e.g. if you have data on 100,000 patients with liver disease, its likely that this contains enough diversity of information such that a model might be able to generalize to new patients). This lead to the intersection of survival analysis, big data and machine learning. For example, the "random forest" algorithm has been extended (Ishwaran et al., 2008) to the survival analysis context, where it can be used to make predictions on the individual level:

#load libraries
library(survival)
library(dplyr)
library(ranger)
library(data.table)
library(ggplot2)

#use the built in "lung" data set
#remove missing values (dataset is called "a")

a = na.omit(lung)

#create id variable
a$ID <- seq_along(a[,1])

#create test set with only the first 3 rows
new = a[1:3,]

#create a training set by removing first three rows
a = a[-c(1:3),]

#fit survival model (random survival forest)
r_fit <- ranger(Surv(time,status) ~ age + sex + ph.ecog + ph.karno + pat.karno + meal.cal + wt.loss, data = a, mtry = 4, importance = "permutation", splitrule = "extratrees", verbose = TRUE)

#create new intermediate variables required for the survival curves

death_times <- r_fit$unique.death.times
surv_prob <-data.frame(r_fit$survival)
avg_prob <- sapply(surv_prob, mean)

#use survival model to produce estimated survival curves for the first three observations

pred <- predict(r_fit, new, type = 'response')$survival
pred <- data.table(pred)
colnames(pred) <- as.character(r_fit$unique.death.times)

#plot the results for these 3 patients

plot(r_fit$unique.death.times, pred[1,], type = "l", col = "red")
lines(r_fit$unique.death.times, pred[2,], type = "l", col = "green")
lines(r_fit$unique.death.times, pred[3,], type = "l", col = "blue")

Optional: Variable Importance (not the same as Variable "Effect") and Model Fit

 vi <- data.frame(sort(round(r_fit$variable.importance, 4), decreasing = TRUE))
 names(vi) <- "importance"

 head(vi)
          importance
ph.ecog       0.0296
sex           0.0205
pat.karno     0.0123
ph.karno      0.0075
wt.loss       0.0022
age          -0.0012

cat("Prediction Error = 1 - Harrell's c-index = ", r_fit$prediction.error)

#note: in this particular example the model happens to perform poorly (c-index = 0 is bad, c-index = 1 is good)

Prediction Error = 1 - Harrell's c-index =  0.3800807

Conclusion: Thus, when "enough" data is available (and if you are ready to assume a certain amount of "risk") : Can regression models be used to make predictions and inferences at the individual level? In the above picture, assuming we trust our model (e.g. if the c-value was higher), could we infer that the patient with the "Red Curve" is likely to survive longer than the patients with the "Blue Curve" and "Green Curve" - therefore, in some scenarios, it might be more advisable to prioritize the treatment of the patients corresponding to the "Blue Curve" and "Green Curve"?

Thanks!

References:

• https://rviews.rstudio.com/2017/09/25/survival-analysis-with-r/
• https://en.wikipedia.org/wiki/Reference_class_problem
• https://en.wikipedia.org/wiki/Ecological_fallacy

Answer 1

The main distinction between "individual" and "cohort" predictions has more to do with what you do after you build the model. Survival models are necessarily built on data from cohorts. They can nevertheless make predictions for an individual's survival function based on that individual's set of covariate values. See any survival software's predict() or similar function.

A major goal of clinical survival modeling is to try to estimate survival for individual patient counseling or to pick the therapy that has the best chance of extending survival for that particular patient. Consider TNM (tumor-node-metastasis) staging, which might distinguish 4 groups with progressively shorter overall survival after diagnosis of a particular type of cancer. That's based on very large cohorts, and regularly re-evaluated by clinical experts. It's not very precise for any individual patient, but at least it provides a rough answer to the patient's question: "How long do I have to live?" From that perspective, survival modeling has always been about individual survival--at least in terms of trying to figure out which cohort best describes an individual's survival probability.

Chapter 10 of Therneau and Grambsch, now 21 years old, discusses "Expected Survival" estimates at both the individual and the cohort level, with individual survival estimates discussed first. For example, although you will sometimes hear it claimed that Cox models aren't intended for individual predictions, the combination of a model-estimated baseline hazard with regression coefficients does provide predictions of survival curves based on individual sets of covariates, as discussed in that chapter. For "cohort" analysis, as explained there, the problem is perhaps more difficult, as you need to make a decision about what type of "average" to make among members of the cohort.

Larger data sets might provide more refined estimates of regression coefficients or allow models to incorporate larger numbers of predictors. That might in principle make estimates of things like median survival more precise than previously.

Larger data sets, however, do not overcome the fundamental problem with making individual predictions in survival analysis: the typically wide distribution of survival times around the median or other central measure. Survival models effectively start with some baseline form of a survival function and then examine the associations of covariates with properties of that function, for example its position in time, its steepness or shape. You model the entire survival function. A highly precise estimate of an individual's survival would require something close to a step function in the predicted survival function for that individual. In clinical survival functions, at least, you don't typically find step functions, except perhaps for those with very poor prognosis.

The variability in clinical survival outcomes at any given set of covariate values is inherent in proportional-hazards (PH) models, like Cox regressions. The baseline survival function for a PH model runs over the entire period from time = 0 to the last event time in the data set. In clinical survival studies that period can be years or decades. An individual with a "favorable" set of covariate values could in principle, under the PH assumption, experience the event at any time during that long period, just with a lower risk of event at any time than an individual with "less-favorable" covariate values. You thus don't get a precise individual survival prediction from a PH model for an individual with "favorable" covariate values, even though you might have a model that describes the survival characteristics of a cohort of such individuals quite well.

You can see the issues also in the formula that describes parametric accelerated failure time (AFT) models in survival analysis,

$$\log T = -x'\beta + \sigma W, $$

where $T$ is survival time, $x$ is the vector of covariates with corresponding coefficients $\beta$, $\sigma$ is a scale factor and $W$ is an error term with a specified distribution (e.g., extreme-value, normal, logistic). The distribution of errors $\sigma W$ limits the precision of the estimated survival time $T$ based on $x' \beta$ for an individual. In practice, that distribution is generally quite wide.

The term $\sigma W$ captures the un-modeled variability in outcomes, similar to the error term in a linear regression. I suppose that if one could incorporate information about all the interactions among 20,000 genes and environmental, social and economic conditions one might be able to model more of the variability. Yes, that would require very large data sets and would probably benefit from "machine learning" approaches. That still would not, however, capture variability arising from unpredictable accidents or infections--as the last 2 years should make very clear.