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:
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"?