# Cohort Regression Models vs Individual Regression Models

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)

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!

Note: I have seen some attempts of estimating individual survival curves over here (https://arxiv.org/pdf/1811.11347.pdf)

References:

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.

• 2nd question: have you heard of the random survival forest?arxiv.org/abs/0811.1645 supposedly, the random survival forest is able to effectively use "bootstrap aggregation" and better estimate individual survival times? Oct 22, 2021 at 23:28
• @stats555 if there were no $\sigma W$ term, a survival function would be a step function dropping from 100% to 0% survival at the time predicted by $x’ \beta$. How many survival functions have you seen that look like such step functions? Random forests and boosted survival trees have strengths in terms of incorporating multiple predictors and interactions without overfitting, but I haven’t seen any that produce anything like such a step-shaped survival function.
– EdM
Oct 23, 2021 at 2:28
• @ EdM : thank you for your reply! I am not denying the fact that sigmaw exists - I am just curious to know its size relative to x-transposebeta. How do we know that x-transpose*beta is much smaller than sigma *w? Thank you! Oct 23, 2021 at 7:14
• Side question: do you think any survival models in the future might be able to solve the individual inference problem? Maybe when RNN (recurrent neural networks) are used in the survival context? Thanks! Oct 23, 2021 at 7:16

In a non-deterministic system, we cannot know with absolute certainty the future state of any individual. Increasing the number of individuals, measuring more of their characteristics, precisely measuring outcomes, or using different methods to discern the association between characteristics and outcomes do not have the emergent property of allowing us to clearly known the exact future state of individuals. These approaches only allow us to asymptotically approach individual prediction by creating more accurate cohorts; we never have access to an individual's future outcome.

Semantics are inescapable here, and I fully accept that based on your definition of "individual prediction" I may have answered an unrelated question. A prediction is a testable claim about an exact future state (eg, "you will have a heart attack"), while risk describes the likelihood of an outcome (eg, there is a 24% chance you will have a heart attack). Testable predictions and risk estimates should be bounded by time.

For example, we would never say to a patient "If you undergo this treatment, you will be cancer-free on July 12,2022." We prefer to say "In patients similar to you who undergo this treatment, 50% were cancer-free within 6 months." We simply do not have the information to provide the certainty required for the individual-level prediction, while we have more confidence in the accuracy of the cohort-derived risk estimate.

When considering your survival curves, these curves are estimates of risk that, by definition, are derived from and pertain to groups of patients. We may be able to draw more accurate curves from larger samples and better models, but individuals cannot have a risk curve; they either had or did not have an event by a given time.

In short:

Can we make predictions about the future state of an individual? Yes, but this is an approach fraught with problems of certainty.

Can we estimate risk more precisely for an individual? Yes, and this leaves appropriate space for uncertainty.

Your question is not new. In 1814, Simon Laplace wrote:

We may regard the present state of the universe as the effect of its past and the cause of its future. An intellect which at a certain moment would know all forces that set nature in motion, and all positions of all items of which nature is composed, if this intellect were also vast enough to submit these data to analysis, it would embrace in a single formula the movements of the greatest bodies of the universe and those of the tiniest atom; for such an intellect nothing would be uncertain and the future just like the past would be present before its eyes.

However, the advent of quantum mechanics years later supports the notion that the universe is largely probabilistic and not deterministic. Machine learning is unlikely to be the "intellect" that Laplace describes.

Your question, as I've interpreted it, may be more philosophic than methodologic, and accordingly, definitive sources are scarce. Thus, there are two books that may help you with the distinction between experiments with one trial or observation in time (eg, an individual) and those with multiple trials or observations over time (eg, aggregates of individuals):

Tetlock, P. E., & Gardner, D. (2015). Superforecasting: The art and science of prediction.

Lewis, M. (2003). Moneyball: The art of winning an unfair game.