Overview

This talk draws heavily on Professor Galit Shmueli's 2010 paper of the same name:

Shmueli, G. (2010), To Explain or To Predict?, Statistical Science, vol 25 no 3, pp. 289-310.

"Since all models are wrong, the scientist must be alert to what is importantly wrong. It is inappropriate to be concerned about mice when there are tigers abroad" Box (1976)

The two broad classes of DS/modelling question:

Explain

• How does an explanatory variable(s) relate to an outcome?
• We might consider two main groups:
  • Description
  • Causal Inference / Counter factual analysis

Predict

• How can predictor variable(s) predict an outcome?
• At future time points or in different context?

You can use many of the same models to fit in either context, but how you do it is different!

Grounding

$$\mathcal{Y} = \mathcal{F(X)}$$

• Imagine that $\mathcal{X}$ causes $\mathcal{Y}$ through some function: $\mathcal{F}$
• $\mathcal{F}$ is a theoretical model, set of statements, path model etc.
• To model with data, we need to use measurable variables

$$E(Y) = f(X)$$

Our modelling goals are

• Explanatory: using various $f$ to estimate $\mathcal{F}$, using $X,Y$
• Predictive: estimate new values of $Y$, using $f(X)$

Example of a Causal Explanatory model

• Directed Acyclic graph (DAG)

Resued from: Arnold et al. Int J Epidemiol, Volume 49, Issue 6, December 2020, Pages 2074–2082, https://doi.org/10.1093/ije/dyaa049

Simple Example:

We will use an example I sourced from Kaggle, related to the paper:

Chicco, D., Jurman, G. Machine learning can predict survival of patients with heart failure from serum creatinine and ejection fraction alone. BMC Med Inform Decis Mak 20, 16 (2020). https://doi.org/10.1186/s12911-020-1023-5

https://www.kaggle.com/datasets/andrewmvd/heart-failure-clinical-data/data

• Data related to heart failure
• Paper tests various prediction methods to see if albumin and serum creatinine alone can predict death.
• Any modelling approach would be more more extensive than the example here.

I will be using logistic regression in both context

Regression in 2 minutes...

Regression models (1)

$$y= \alpha + \beta x + \epsilon$$

Regression models (2)

Regression equation

$$\large{y= \alpha + \beta_{i} x_{i} + \epsilon}$$

How can we use regression on other distributions / data types?

We can 'zoom out' to a more the 'Generalized Linear Model (GLM):

For distributions in the exponential family, GLM allows the linear model to relate to response through a function:

$$\large{g(\mu)= \alpha + \beta_{i} x_{i}}$$

• Where $\mu$ is the expectation of $Y$,
• $g$ is the link function - related to each
  

So, for logistic regression, the link function is 'logit' or the log odds of the event.

Explanatory model - R

r_model_exp <- glm(DEATH_EVENT ~ serum_creatinine + ejection_fraction
              , data=heart_failure_dt
              , family = "binomial")
summary(r_model_exp)

## 
## Call:
## glm(formula = DEATH_EVENT ~ serum_creatinine + ejection_fraction, 
##     family = "binomial", data = heart_failure_dt)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.37769    0.54343   0.695    0.487    
## serum_creatinine   0.74987    0.17932   4.182 2.89e-05 ***
## ejection_fraction -0.05986    0.01350  -4.435 9.19e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 375.35  on 298  degrees of freedom
## Residual deviance: 324.32  on 296  degrees of freedom
## AIC: 330.32
## 
## Number of Fisher Scoring iterations: 4

Explanatory model - Python

import statsmodels.formula.api as smf
py_model1_exp = smf.logit("DEATH_EVENT ~ serum_creatinine + ejection_fraction", data=heart_failure_pd).fit()

print(py_model1_exp_summary)

##                            Logit Regression Results                           
## ==============================================================================
## Dep. Variable:            DEATH_EVENT   No. Observations:                  299
## Model:                          Logit   Df Residuals:                      296
## Method:                           MLE   Df Model:                            2
## Date:                Wed, 20 Nov 2024   Pseudo R-squ.:                  0.1359
## Time:                        13:39:05   Log-Likelihood:                -162.16
## converged:                       True   LL-Null:                       -187.67
## Covariance Type:            nonrobust   LLR p-value:                 8.308e-12
## =====================================================================================
##                         coef    std err          z      P>|z|      [0.025      0.975]
## -------------------------------------------------------------------------------------
## Intercept             0.3777      0.543      0.695      0.487      -0.687       1.443
## serum_creatinine      0.7499      0.179      4.181      0.000       0.398       1.101
## ejection_fraction    -0.0599      0.013     -4.435      0.000      -0.086      -0.033
## =====================================================================================

Testing Fit

• Significance of coefficients in our model summaries
• Assumption of regression being met - a topic for another day

library(ModelMetrics)
auc(r_model_exp)

## [1] 0.7614173

from sklearn import metrics
py_auc = metrics.roc_auc_score(heart_failure_pd['DEATH_EVENT'], py_model1_exp.fittedvalues)
print(py_auc)

## 0.7614172824302136

Is over-fitting an issue?

No, not if your goal is explanatory

Why statsmodels ?

https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression

Predictive Model - R

heart_failure_dt$sc_serum_creatinine <- scale(heart_failure_dt$serum_creatinine)
heart_failure_dt$sc_ejection_fraction <- scale(heart_failure_dt$ejection_fraction)
trainIndex <- caret::createDataPartition(heart_failure_dt$DEATH_EVENT
                                         , p = .8
                                         , list = FALSE
                                         , times = 1)
Train <- heart_failure_dt[ trainIndex,]
Test  <- heart_failure_dt[-trainIndex,]
r_model_pred <- glm(DEATH_EVENT ~ sc_serum_creatinine + sc_ejection_fraction
              ,  data=Train
              , family = "binomial")
predictions <- predict(r_model_pred, newdata=Test, type="response")
# Model performance metrics
auc(Test$DEATH_EVENT, predictions)

## [1] 0.8339599

Predictive model - Python

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
sc= StandardScaler()
X = sc.fit_transform(heart_failure_pd[['serum_creatinine', 'ejection_fraction']])
y = heart_failure_pd[['DEATH_EVENT']]
X_train, X_test, y_train, y_test = train_test_split(X,y , test_size = 0.2)
log_reg = LogisticRegression(penalty = 'None').fit(X_train,y_train)

y_pred = log_reg.predict_proba(X_test)
py_auc = metrics.roc_auc_score(y_test, y_pred[:,1])
print(py_auc)

## 0.7957351290684623

What is important in assessing the fit:

Explain or predict Bingo (1):

Forecasting attendances at an Emergency Department

Predict!

Explain or predict Bingo (2):

What drives people to attend an Emergency Department?

Explain!

Explain or predict Bingo (3):

What is a person's risk of attending an Emergency Department?

It depends:... is it about the person's individual risk based on explanatory factors, the best prediction you can make, or is it for risk-adjustment?

Explain or predict Bingo (4):

Did our new UTC pathway decrease how often people attend the Emergency Department?

Explain!

Explain or predict Bingo (5):

Building a Large Language Model to answer questions as a chatbot

Predict!

Explain or predict Bingo (6):

Modelling long-term population health-state changes

It depends:... are you testing what causes it, or predicting future states of the population?

Summary

References

Arnold, K.F. et al. (2020) ‘Reflection on modern methods: generalized linear models for prognosis and intervention—theory, practice and implications for machine learning’, International Journal of Epidemiology, 49(6), pp. 2074–2082. Available at: https://doi.org/10.1093/ije/dyaa049.

Box, G.E.P. (1976) “Science and Statistics.” Journal of the American Statistical Association 71, no. 356: 791–99. https://doi.org/10.2307/2286841.

Chicco, D. and Jurman, G. (2020) ‘Machine learning can predict survival of patients with heart failure from serum creatinine and ejection fraction alone’, BMC Medical Informatics and Decision Making, 20(1), p. 16. Available at: https://doi.org/10.1186/s12911-020-1023-5.

Hernán, M. A., Hsu, J. and Healy, B. (2019) ‘A Second Chance to Get Causal Inference Right: A Classification of Data Science Tasks’, CHANCE, 32(1), pp. 42–49. doi: 10.1080/09332480.2019.1579578.

Shmueli, G. (2010) 'To Explain or to Predict?' Statistical Science 25, no. 3 : 289–310. http://www.jstor.org/stable/41058949.

Predictive Model - R bonus (like Scikit learn assumes you want...)

heart_failure_dt$sc_ejection_fraction <- scale(heart_failure_dt$ejection_fraction)
trainIndex <- caret::createDataPartition(heart_failure_dt$DEATH_EVENT
                                         , p = .8
                                         , list = FALSE
                                         , times = 1)
Train <- heart_failure_dt[ trainIndex,]
Test  <- heart_failure_dt[-trainIndex,]
x <- model.matrix(DEATH_EVENT~log(serum_creatinine)+sc_ejection_fraction, Train)[,-1]
y <- Train$DEATH_EVENT
library(glmnet)
# Cross validate to get best lambda (shrinkage)
cv <- cv.glmnet(x, y, alpha = 0, family="binomial")
ridge1<-glmnet(x,y, alpha=0, lamda=cv$lambda.1se, family="binomial")
# Make predictions on the test data
x.test <- model.matrix(DEATH_EVENT~log(serum_creatinine)+sc_ejection_fraction, Test)[,-1]
predictions <- predict(ridge1, newx=x.test, type="response") |> as.vector()
ModelMetrics::auc(Test$DEATH_EVENT, predictions)

## [1] 71.20213