RAPID GLM — Tutorial Notebook¶
This notebook walks you through using the RAPID glm pipeline from end to end.
By the end of this tutorial you will know how to:
- Instantiate a glm model using the RAPID Pipeline Factory
- Fit the model and customise variable labels
- Interpret every assumption test, performance metric, and plot
- Choose the right distributional family and link function for your outcome
This pipeline is designed for use in epidemiological and clinical research contexts, where outcome variables are continuous — such as length of hospital stay, biomarker levels, physiological measurements, or symptom scores.
1. Setup and Installation¶
Before starting, ensure the RAPID package is installed in your environment.
pip install isaric-0.1.0-py3-none-any.whl
Then import the factory and any other dependencies you need.
!pip install gdown
import gdown
url = "https://drive.google.com/uc?export=download&id=1OpP-r3YzKGTKYHijHE-yAgFNYLG-5vj-"
gdown.download(url, 'isaric-0.1.0-py3-none-any.whl', quiet=False)
!pip install isaric-0.1.0-py3-none-any.whl
from isaric.pipelines.pipeline_factory import RAPID_PipelineFactory
import pandas as pd
2. The Pipeline Factory¶
All RAPID pipelines are created through the RAPID_PipelineFactory. Rather than importing and instantiating regression classes directly, the factory provides a single, consistent entry point for creating any supported pipeline by name.
This keeps your code cleaner and ensures you are always using a correctly configured pipeline.
# Instantiate the factory — no arguments needed
factory = RAPID_PipelineFactory()
# See all pipelines available out of the box
print(factory.available())
To create a glm model, call factory.create() with "glm" as the first argument, followed by your data and modelling parameters.
3. Preparing Your Data¶
The pipeline accepts a pandas DataFrame. Rows with missing values are dropped automatically before fitting, so no manual imputation is required unless you want more control over that process.
Your data should have:
- One column representing the continuous outcome you want to model (e.g. length of stay, CRP level)
- One or more columns representing the independent_vars (e.g. age, sex, comorbidity score)
Below is an example using a simulated clinical dataset.
import numpy as np
np.random.seed(42)
n = 500
df = pd.DataFrame({
"length_of_stay": np.random.gamma(shape=3, scale=2, size=n),
"age": np.random.randint(18, 90, size=n),
"sex": np.random.choice([0, 1], size=n),
"comorbidity_score": np.random.poisson(lam=2, size=n),
"icu_admission": np.random.choice([0, 1], size=n, p=[0.7, 0.3])
})
df.head()
4. Creating and Fitting the Model¶
Pass your DataFrame, outcome variable name, and list of independent_vars to factory.create(). The regression_type parameter controls whether this is a univariable ("Uni") or multivariable ("Multi") regression — this affects the column names in the results table but not the underlying model.
model = factory.create(
"glm",
data=df,
dependent_var="length_of_stay",
independent_vars=["age", "sex", "comorbidity_score", "icu_admission"],
regression_type="Multi"
)
Fitting the Model¶
Call .fit() to run the regression. You can optionally pass a labels dictionary to map raw column names to human-readable display names — these appear in all result tables and plots.
model.fit(
labels={
"age": "Age (years)",
"sex": "Sex (Male=1)",
"comorbidity_score": "Comorbidity Score",
"icu_admission": "ICU Admission"
},
cross_val=True,
n_splits=5
)
During fit(), the pipeline automatically:
- Fits the GLM
- Runs all assumption tests
- Computes all performance metrics
- Runs k-fold cross-validation (if
cross_val=True)
5. Results Table¶
The coefficient table is stored in model.summary_df. Each row represents one predictor. The columns show the regression coefficient, its 95% confidence interval, and the p-value.
model.summary_df
How to interpret the results table¶
- Coefficient: The estimated change in the outcome for a one-unit increase in the predictor, holding all other independent_vars constant. A positive value means the outcome increases; a negative value means it decreases.
- LowerCI / UpperCI: The 95% confidence interval around the coefficient. If this interval does not cross zero, the association is statistically significant at p < 0.05.
- p-value: The probability of observing a coefficient this large (or larger) if the true coefficient were zero. Conventionally, values below 0.05 are considered statistically significant — but always interpret in clinical context.
6. Assumption Tests¶
glm relies on several assumptions about the data and residuals. Violating these assumptions can invalidate your results. The pipeline tests all of them automatically.
Display all assumption results with:
model.summary(assumptions="all")
Or access the underlying dataframe directly:
model.assumption_metrics_df
6.1 Durbin-Watson — Independence of Errors¶
What it tests: Whether the residuals (errors) from the model are independent of each other, or whether there is autocorrelation — i.e. the error for one observation predicts the error for another.
Why it matters in epidemiology: Autocorrelation commonly arises in longitudinal data, time-series data, or clustered data (e.g. patients within hospitals). If residuals are correlated, your standard errors will be underestimated and your confidence intervals will be too narrow.
How to interpret:
| Value | Interpretation |
|---|---|
| Close to 2 | Residuals are independent — assumption is met |
| Below 1.5 | Positive autocorrelation — residuals in one direction tend to be followed by residuals in the same direction |
| Above 2.5 | Negative autocorrelation — residuals alternate direction |
What to do if violated: Consider whether your data has a natural ordering (time, geography, clustering). If so, you may need a mixed-effects model or GEE (Generalised Estimating Equations) to account for correlation structure.
# Access the Durbin-Watson statistic directly
print(f"Durbin-Watson: {model.dw:.3f}")
# Or display only this test
model.summary(assumptions=["Durbin-Watson"])
6.2 Shapiro-Wilk — Normality of Residuals¶
What it tests: Whether the residuals from the model follow a normal distribution.
Why it matters in epidemiology: glm assumes that the errors are normally distributed. If this assumption is violated, p-values and confidence intervals may be unreliable, particularly in small samples. In large samples (n > 200), the central limit theorem means this assumption becomes less critical.
How to interpret:
| p-value | Interpretation |
|---|---|
| > 0.05 | Fail to reject H₀ — residuals are consistent with normality |
| ≤ 0.05 | Reject H₀ — residuals are not normally distributed |
Important caveat: In large epidemiological datasets, the Shapiro-Wilk test is very sensitive and will often reject normality even for minor, inconsequential deviations. Always inspect the Q-Q plot alongside this test result.
What to do if violated: Consider transforming the outcome variable (e.g. log transformation for right-skewed data), or switching to a GLM family that better matches your outcome's distribution (see Section 9).
print(f"Shapiro-Wilk statistic: {model.shapiro_wilk_test_statistic:.4f}")
print(f"Shapiro-Wilk p-value: {model.shapiro_wilk_p_value:.4f}")
model.summary(assumptions=["Shapiro-Wilk Statistic", "Shapiro-Wilk p-value"])
6.3 VIF — Multicollinearity¶
What it tests: Whether predictor variables are highly correlated with each other (multicollinearity). The Variance Inflation Factor (VIF) measures how much the variance of a coefficient is inflated due to its correlation with other independent_vars.
Why it matters in epidemiology: High multicollinearity makes it difficult to isolate the individual effect of each predictor. Coefficients become unstable — small changes in the data can produce large changes in the estimates. This is a common problem when modelling correlated risk factors (e.g. BMI, weight, and waist circumference together).
How to interpret:
| VIF | Interpretation |
|---|---|
| 1–5 | Acceptable — low to moderate correlation |
| 5–10 | Potential concern — consider whether both variables are needed |
| > 10 | Serious multicollinearity — coefficients are unreliable |
What to do if violated: Remove one of the correlated independent_vars, combine them into a composite score, or use dimensionality reduction.
# View VIF for all independent_vars
model.vif_df
# Or display via summary with a custom threshold
model.summary(assumptions=["VIF"], vif_threshold=5.0)
6.4 Influential Outliers — Cook's Distance¶
What it tests: Whether any individual observations have a disproportionate influence on the model's coefficient estimates. Cook's distance summarises how much all fitted values change when a given observation is removed.
Why it matters in epidemiology: A single unusual patient — for example, an extreme outlier in length of stay due to a rare complication — could distort the estimated associations for the entire cohort. Identifying such observations is essential before drawing clinical conclusions.
How to interpret:
The threshold is set at 4/n (where n is the sample size). Observations with Cook's distance above this threshold are flagged as potentially influential.
This does not mean these observations should be automatically removed. Instead:
- Investigate whether they represent data entry errors
- Consider whether they are clinically meaningful extreme cases
- Re-run the model without them to assess sensitivity
What to do if influential points are found: Report them transparently. If they represent genuine data, consider robust regression methods.
print(f"Threshold (4/n): {model.influential_outliers_threshold:.4f}")
print(f"Number of influential points: {len(model.influential_points)}")
print(f"Indices: {model.influential_points}")
model.summary(assumptions=["Influential Outliers", "Number of Influential Points"])
7. Diagnostic Plots¶
Visual diagnostics are an essential complement to the numerical assumption tests above. Always inspect these plots — numbers alone can miss important patterns.
7.1 Residuals vs Fitted Values¶
What it shows: The residuals (observed minus predicted values) plotted against the fitted (predicted) values.
What to look for:
- Ideal: Points scattered randomly around the horizontal line at zero, with no discernible pattern
- Funnel shape (heteroscedasticity): Residuals spread out as fitted values increase — this violates the assumption of constant variance and suggests a transformation or a different GLM family may be needed
- Curved pattern (non-linearity): A U-shape or arch suggests the relationship between predictor and outcome is not linear — consider adding polynomial terms or transforming a predictor
In epidemiological practice: A funnel shape is very common with outcomes like length of stay or biomarker levels, which are right-skewed and have variance that increases with the mean. This is a strong signal to consider the Gamma or Inverse Gaussian family (see Section 9).
model.summary(plots=["residuals_vs_fitted"])
7.2 Q-Q Plot (Quantile-Quantile)¶
What it shows: The quantiles of the model's residuals plotted against the theoretical quantiles of a normal distribution. If residuals are perfectly normal, all points fall on the diagonal reference line.
What to look for:
- Points following the diagonal: Normality assumption is met
- S-shaped curve: The distribution is heavier-tailed than normal (common with skewed clinical outcomes)
- Points deviating at the upper tail: Right skew — typical of outcomes like length of stay or viral load
- Points deviating at both tails: Heavy tails — the distribution has more extreme values than a normal distribution would predict
In epidemiological practice: Some deviation in the tails is very common with real-world clinical data and does not necessarily invalidate your results, especially in larger samples. Use this plot alongside the Shapiro-Wilk test.
model.summary(plots=["qq_plot"])
7.3 Forest Plot¶
What it shows: Coefficient estimates and their 95% confidence intervals displayed as horizontal lines with central point estimates. A vertical reference line is drawn at zero (the null value — no association).
What to look for:
- Confidence interval entirely to the right of zero: Significant positive association
- Confidence interval entirely to the left of zero: Significant negative association
- Confidence interval crossing zero: No statistically significant association at p < 0.05
- Width of the interval: Narrow intervals indicate more precise estimates (larger sample or stronger signal); wide intervals indicate uncertainty
In epidemiological practice: The forest plot gives an immediate visual overview of all effect sizes in the model, making it easy to compare the relative magnitude and direction of different independent_vars at a glance.
model.summary(plots=["forest_plot"])
8. Performance Metrics¶
Performance metrics quantify how well the model fits the data. Display all of them with:
model.summary(performance="all")
# Or access the dataframe directly
model.performance_metrics_df
8.1 Error Metrics — MSE, RMSE, MAE¶
These three metrics all measure prediction error in terms of the outcome's original scale.
Mean Squared Error (MSE): The average of squared differences between observed and predicted values. Squaring penalises large errors more heavily. Useful for comparing models, but not directly interpretable in outcome units.
Root Mean Squared Error (RMSE): The square root of MSE. This is in the same units as your outcome variable, making it directly interpretable. For example, if your outcome is length of stay in days, an RMSE of 3.2 means your model's predictions are off by approximately 3.2 days on average.
Mean Absolute Error (MAE): The average of absolute differences between observed and predicted values. Less sensitive to large outliers than RMSE. In clinical contexts where extreme cases are common, MAE often gives a more representative picture of typical prediction error.
print(f"MSE: {model.mse:.4f}")
print(f"RMSE: {model.rmse:.4f}")
print(f"MAE: {model.mae:.4f}")
model.summary(performance=["MSE", "RMSE", "MAE"])
8.2 R² and Adjusted R²¶
R² (Coefficient of Determination): The proportion of variance in the outcome explained by the model. A value of 0.40 means the model explains 40% of the variability in the outcome.
Adjusted R²: R² penalised for the number of independent_vars in the model. Adding more independent_vars always increases R², even if they have no real relationship with the outcome. Adjusted R² corrects for this and should always be preferred over R² when comparing models with different numbers of independent_vars.
Interpreting R² in epidemiology: R² values in clinical and epidemiological research are often lower than in other fields — values of 0.10–0.30 are common and can still represent meaningful findings. Human health outcomes are inherently variable and influenced by many unmeasured factors. A low R² does not mean your independent_vars are unimportant.
print(f"R²: {model.r2:.4f}")
print(f"Adjusted R²: {model.adjusted_r2:.4f}")
8.3 Pseudo R² Metrics — McFadden, Efron¶
Because RAPID uses a GLM framework, it also computes pseudo R² metrics derived from log-likelihoods. These are especially relevant when using non-Gaussian families.
McFadden R²: Compares the log-likelihood of your fitted model to the log-likelihood of a null model (intercept only). Values between 0.2 and 0.4 are generally considered good fit in epidemiological models.
Adjusted McFadden R²: Penalises for model complexity — equivalent to Adjusted R² but in the likelihood framework.
Efron R²: Computed from residuals relative to the outcome mean. For a Gaussian model with identity link, this equals the standard R². It is more informative when using non-Gaussian families.
print(f"McFadden R²: {model.mcfadden_r2:.4f}")
print(f"Adjusted McFadden R²: {model.mcfadden_adj_r2:.4f}")
print(f"Efron R²: {model.efron_r2:.4f}")
8.4 Information Criteria — AIC and BIC¶
AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) are used primarily for model comparison, not for evaluating a single model in isolation. Lower values indicate a better trade-off between model fit and complexity.
AIC penalises model complexity lightly and is preferred when prediction is the main goal.
BIC penalises complexity more heavily (especially in large samples) and is preferred when parsimony and interpretability are priorities — which is common in epidemiological modelling.
Rule of thumb: A difference in AIC or BIC of more than 10 between two models is considered strong evidence in favour of the model with the lower value.
print(f"AIC: {model.aic:.4f}")
print(f"BIC: {model.bic:.4f}")
print(f"LLF: {model.llf:.4f}")
9. Cross-Validation¶
K-fold cross-validation provides an estimate of how well the model generalises to new, unseen data. The dataset is split into k folds; the model is trained on k-1 folds and evaluated on the held-out fold, repeating until each fold has been used as the test set once.
Why this matters in epidemiology: A model can appear to fit the training data very well but perform poorly on new patients. Cross-validation is particularly important when your dataset is small or when you plan to deploy the model for prediction in a clinical setting.
model.summary(cross_val="all")
# Or access directly
model.cv_df
print(f"Mean CV MSE: {model.cv_mse_scores.mean():.4f}")
print(f"Std CV MSE: {model.cv_mse_scores.std():.4f}")
Interpreting Cross-Validation Results¶
- Mean CV MSE: The average prediction error across all folds. Compare this to the training MSE — if the CV MSE is much higher, the model may be overfitting.
- Standard Deviation of CV MSE: Measures stability across folds. A high standard deviation suggests the model's performance varies considerably depending on which data it sees, which may indicate instability — particularly in smaller datasets.
- Individual Fold MSEs: Inspect these to check whether any single fold is driving a high mean or standard deviation.
10. Choosing a Family and Link Function¶
The default glm (Gaussian family, identity link) assumes your outcome is normally distributed with constant variance. In practice, many continuous outcomes in clinical research do not meet these assumptions.
RAPID supports several GLM families that better reflect the true distribution of common epidemiological outcomes. Choosing the right family can substantially improve model fit, interpretability, and the validity of your inferences.
Pass the desired family and link when creating the model:
model = factory.create(
"glm",
data=df,
dependent_var="length_of_stay",
independent_vars=["age", "sex", "comorbidity_score"],
family="gamma",
link="log"
)
model.fit()
10.1 Gaussian + Identity (Default)¶
Family: "gaussian" | Link: "identity"
The standard glm model. Assumes the outcome is normally distributed around the predicted mean, and that the relationship between independent_vars and outcome is linear and additive.
When to use:
- Outcome is approximately normally distributed (e.g. blood pressure, height, certain physiological measurements)
- Residuals from an initial model appear normally distributed on a Q-Q plot
- Outcome can take negative values
Epidemiological examples: Systolic blood pressure, body temperature, spirometry values (FEV1), cognitive test scores.
Caution: Do not use when your outcome is strictly positive and right-skewed (e.g. length of stay, costs, biomarker levels) — the Gaussian family can predict negative values, which is nonsensical for such outcomes.
model_gaussian = factory.create(
"glm",
data=df,
dependent_var="length_of_stay",
independent_vars=["age", "sex", "comorbidity_score"],
family="gaussian",
link="identity"
)
model_gaussian.fit()
10.2 Gamma + Log¶
Family: "gamma" | Link: "log"
The Gamma family models outcomes that are strictly positive and right-skewed, where the variance increases with the mean. The log link means the model is multiplicative — a one-unit increase in a predictor multiplies the outcome by a factor of exp(coefficient).
When to use:
- Outcome is always positive and positively skewed
- Variance appears to increase with the mean (visible as a funnel shape in your residuals vs fitted plot)
- You want a multiplicative interpretation of effects
Epidemiological examples: Length of hospital stay, healthcare costs, viral load, CRP levels, time-to-event outcomes treated as continuous.
Interpreting coefficients: Exponentiate the coefficient to get the multiplicative effect. A coefficient of 0.05 means a one-unit increase in the predictor multiplies the outcome by exp(0.05) ≈ 1.05, i.e. a 5% increase.
model_gamma_log = factory.create(
"glm",
data=df,
dependent_var="length_of_stay",
independent_vars=["age", "sex", "comorbidity_score"],
family="gamma",
link="log"
)
model_gamma_log.fit()
10.3 Gamma + Inverse¶
Family: "gamma" | Link: "inverse"
Also models strictly positive, right-skewed outcomes, but the inverse link means the model is formulated in terms of the reciprocal of the mean. This is the canonical link for the Gamma family and is mathematically natural, though coefficients are harder to interpret directly.
When to use:
- Same distributional assumptions as Gamma + Log
- Historically preferred in some fields as the canonical Gamma link
- When comparing models: if AIC/BIC is lower than Gamma + Log, prefer this
Epidemiological examples: Same as Gamma + Log. The choice between log and inverse links is often made empirically by comparing model fit using AIC/BIC.
Interpreting coefficients: A positive coefficient means an increase in the predictor is associated with a decrease in the outcome (because the model is on the inverse scale). This unintuitive direction is one reason Gamma + Log is more commonly preferred in applied epidemiological work.
model_gamma_inv = factory.create(
"glm",
data=df,
dependent_var="length_of_stay",
independent_vars=["age", "sex", "comorbidity_score"],
family="gamma",
link="inverse"
)
model_gamma_inv.fit()
10.4 Inverse Gaussian + Inverse¶
Family: "inv_gaussian" | Link: "inverse"
The Inverse Gaussian family models outcomes that are strictly positive and even more heavily right-skewed than the Gamma. The variance increases faster with the mean than in the Gamma family, making it suitable for outcomes with a very long right tail.
When to use:
- Outcome is strictly positive with an extreme right skew
- Gamma family has been tried but residual diagnostics still show poor fit
- Outcomes where very large values (extreme cases) are a genuine feature of the distribution, not anomalies
Epidemiological examples: ICU length of stay (which can be extremely long for a small number of patients), total healthcare costs over a lifetime, time to a rare event.
Practical note: This family is less commonly used than Gamma and can be harder to fit — convergence issues are more likely with small datasets or highly influential observations. Always inspect assumption diagnostics carefully.
model_inv_gauss = factory.create(
"glm",
data=df,
dependent_var="length_of_stay",
independent_vars=["age", "sex", "comorbidity_score"],
family="inv_gaussian",
link="inverse"
)
model_inv_gauss.fit()
10.5 Tweedie + Log¶
Family: "tweedie" | Link: "log"
The Tweedie family is a flexible generalisation that encompasses several other distributions. With a log link, it models outcomes that are non-negative, may include exact zeros, and have a right-skewed positive component. This makes it uniquely suited to outcomes that are zero for some participants and positive for others.
When to use:
- Outcome contains a meaningful proportion of exact zeros alongside positive values
- Examples include healthcare resource utilisation, where many patients use zero resources and a smaller number have very high utilisation
- Modelling costs where a large proportion of patients incur no cost
Epidemiological examples: Number of GP visits (many patients have zero), total medication costs (zero for those on no medications), days of oxygen supplementation (zero for mild cases).
Why not use zero-inflated models? In many epidemiological contexts, the zeros are not a separate process but part of the same underlying distribution. The Tweedie family handles this as a single unified model, which is often more parsimonious.
Interpreting coefficients: As with Gamma + Log, exponentiate coefficients for a multiplicative interpretation.
model_tweedie = factory.create(
"glm",
data=df,
dependent_var="length_of_stay",
independent_vars=["age", "sex", "comorbidity_score"],
family="tweedie",
link="log"
)
model_tweedie.fit()
10.6 Comparing Families¶
When you are unsure which family is most appropriate, fit several models and compare their AIC and BIC. Lower values indicate a better fit, penalised for complexity.
Use the RAPID_ModelComparator to compare performance across models side by side:
from isaric.pipelines.model_comparison import RAPID_ModelComparator
comparator = RAPID_ModelComparator({
"Gaussian (identity)": model_gaussian,
"Gamma (log)": model_gamma_log,
"Gamma (inverse)": model_gamma_inv,
"Tweedie (log)": model_tweedie,
})
comparator.report()
In addition to AIC and BIC, inspect the residuals vs fitted and Q-Q plots for each model. The best-fitting family will show:
- Randomly scattered residuals with no funnel shape
- Points closely following the diagonal on the Q-Q plot
11. Displaying Selected Results¶
You do not need to display everything at once. The summary() method accepts lists of specific metrics so you can focus on what is relevant for your analysis.
# Show only the metrics most relevant for a publication-ready report
model.summary(
assumptions=["Durbin-Watson", "Shapiro-Wilk p-value", "VIF", "Influential Outliers"],
performance=["RMSE", "MAE", "Adjusted R2", "AIC", "BIC"],
cross_val=["Mean CV MSE", "Standard Deviation of CV MSE"],
plots=["forest_plot", "residuals_vs_fitted", "qq_plot"],
vif_threshold=5.0
)
12. Accessing Results Programmatically¶
All results are stored as attributes on the model object and can be accessed directly for downstream use — for example, to export to a report, feed into another analysis, or display in a custom format.
# Coefficient table
model.summary_df
# Performance metrics
model.performance_metrics_df
# Assumption tests
model.assumption_metrics_df
# VIF table
model.vif_df
# Cross-validation results
model.cv_df
# Individual metrics
print(model.r2)
print(model.adjusted_r2)
print(model.aic)
print(model.dw)
print(model.shapiro_wilk_p_value)
print(model.influential_points)
print(model.cv_mse_scores)
Summary¶
This tutorial has covered the full RAPID glm workflow:
| Step | What you did |
|---|---|
| Factory | Created the model via RAPID_PipelineFactory |
| fit() | Fitted the model with labels and cross-validation |
| Results | Interpreted the coefficient table |
| Assumptions | Tested and interpreted Durbin-Watson, Shapiro-Wilk, VIF, and Cook's distance |
| Plots | Interpreted residuals vs fitted, Q-Q plot, and forest plot |
| Performance | Understood MSE, RMSE, MAE, R², AIC, BIC, and pseudo R² metrics |
| Cross-validation | Assessed generalisation using k-fold CV |
| Families | Chose between Gaussian, Gamma, Inverse Gaussian, and Tweedie families |
For the logistic regression equivalent of this tutorial, see the Logistic Regression Tutorial notebook.