Portfolio 2 - Introduction to Survival Analysis in R; A beginners tutorial
2025-10-13
This tutorial aims to teach you how to use survival analysis to in sports data. The tutorial will allow you to:
Understand survival analysis, key principles and when to use survival analysis in sports analytics
Handling of censored data
Kaplan Meier Survival curves
Log Rank tests to compare survival between groups
Cox Proportion Hazards models
Interpretation of hazard ratios
Make predictions
Survival analysis investigates time to event data. Modelling the duration until an event occurs. In sport science this could be time until an athlete’s re injury, time to performance decline (below a certain threshold), career longevity etc.
Key concepts:
Time (t) – duration from study start until the event occurs (T>0)
Event: Binary outcome of interest. For example injury occurs/does not occur
Status (d) – Event indicator. Where d=1 “failed”, event occurred, where d=0 event not observed “censored”.
Censoring
Censoring refers to incomplete data, resulting in an unknown survival time. This will be due to: Study ending prior to event occurring, or if a participant withdraws from the study.
For example if we were investigating injury of cricket players over 10 months. Dan is re injured at 3 months; D=1, t=3 (failed), Mike remains injury free; D=0, t=10 (censored),and Susan retires at month 5; D=0, t=5(censored).
When to use
Time to until event analysis
Censored observations are present
Time intervals between observations may or may not be unequal
Analysing survival rates between groups
When multiple predictors affecting survival time requires adjusting
Do not use when all event times are observed, no censoring occurs, not time based or only concerned about an event occurring (logistic).
Load required packages
library(tidyverse) ## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(survival)
library(survminer) ## Loading required package: ggpubr
##
## Attaching package: 'survminer'
##
## The following object is masked from 'package:survival':
##
## myelomalibrary(ggfortify) Don’t forget to install any packages that aren’t already installed through install.packages(“package_name”)
Import Dataset
# created example dataset of a 20 week study
cricket <- data.frame(
player_id = 1:20,
time = c(8, 12, 6, 15, 3, 9, 14, 7, 11, 5,
18, 10, 4, 13, 16, 2, 8, 12, 9, 14),
status = c(1, 0, 1, 0, 1, 1, 0, 1, 0, 1,
0, 1, 1, 0, 0, 1, 1, 0, 1, 0),
age = c(25, 30, 22, 28, 24, 27, 29, 23, 31, 26,
32, 21, 24, 30, 28, 22, 25, 27, 29, 26),
training_load = c(rep("High", 10), rep("Low", 10))
)Dataset Variables:
player_id: Unique player identifiertime: Weeks until injury (or end of study)status: 1 = injury occurred, 0 = censored (no injury observed)age: Player age in yearstraining_load: Training intensity (High or Low)
Exploratory Data Analysis
#load the data and examine.
head(cricket)## player_id time status age training_load
## 1 1 8 1 25 High
## 2 2 12 0 30 High
## 3 3 6 1 22 High
## 4 4 15 0 28 High
## 5 5 3 1 24 High
## 6 6 9 1 27 High# Summary statistics
summary(cricket)## player_id time status age
## Min. : 1.00 Min. : 2.00 Min. :0.00 Min. :21.00
## 1st Qu.: 5.75 1st Qu.: 6.75 1st Qu.:0.00 1st Qu.:24.00
## Median :10.50 Median : 9.50 Median :1.00 Median :26.50
## Mean :10.50 Mean : 9.80 Mean :0.55 Mean :26.45
## 3rd Qu.:15.25 3rd Qu.:13.25 3rd Qu.:1.00 3rd Qu.:29.00
## Max. :20.00 Max. :18.00 Max. :1.00 Max. :32.00
## training_load
## Length:20
## Class :character
## Mode :character
##
##
## # Count events vs censored
table(cricket$status)##
## 0 1
## 9 11# Distribution by group
table(cricket$training_load, cricket$status)##
## 0 1
## High 4 6
## Low 5 5# Summary by training load
cricket %>%
group_by(training_load) %>%
summarise(
n = n(),
n_events = sum(status),
n_censored = sum(status == 0),
mean_time = mean(time)
)## # A tibble: 2 × 5
## training_load n n_events n_censored mean_time
## <chr> <int> <dbl> <int> <dbl>
## 1 High 10 6 4 9
## 2 Low 10 5 5 10.6Through these EDA’s we can see the dataset contains 20 cricket
players with an average time to injury of 9.8 weeks. Out of 20 players,
11 (55%) experienced injury within the study period, whilst the
remaining 45% were censored.
When comparing training load groups, high training load encountered 6
out of 10 players were injured whilst the low training group had 5
injuries out of 10 players.
High load players also show a shorter meant time till injury (9 weeks vs
10.6 weeks).
These reported statistics dictate our use of survival analysis to test whether training load significantly affects injury risk
Visualisations
# Plot 1: Distribution of survival times
ggplot(cricket, aes(x = time, fill = factor(status))) +
geom_histogram(bins = 10, alpha = 0.7, position = "identity") +
scale_fill_manual(values = c("0" = "blue", "1" = "red"),
labels = c("Censored", "Injured")) +
labs(title = "Distribution of Time to Injury",
x = "Time (weeks)",
y = "Count",
fill = "Status") +
theme_minimal()
The histogram shows that injuries (red) and censored observations (blue) occurred throughout the study period validating our use of survival analysis to appropriately handle the incomplete observations.
Survival Analysis
The Kaplan-Meier method estimates the probability of subjects survival probability at each event time. The curve begins at 100% event-free and decreases in steps steeper drops indicate higher event rate. This is known as the step function where probability drops when events occur, and stays flat between events. . Confidence intervals widen over time as fewer subjects remain at risk.
Understanding the output:
time: when events occurred
n.risk: Number still at risk at that time
n.event: Number of events at that time
survival: Cumulative probability of remaining event-free
95% CI: Confidence interval around the survival estimate
# Kaplan-Meier Survival ====
# KM for all players
km_overall <- survfit(Surv(time, status) ~ 1, data = cricket)
# Summary
summary(km_overall)## Call: survfit(formula = Surv(time, status) ~ 1, data = cricket)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 20 1 0.95 0.0487 0.859 1.000
## 3 19 1 0.90 0.0671 0.778 1.000
## 4 18 1 0.85 0.0798 0.707 1.000
## 5 17 1 0.80 0.0894 0.643 0.996
## 6 16 1 0.75 0.0968 0.582 0.966
## 7 15 1 0.70 0.1025 0.525 0.933
## 8 14 2 0.60 0.1095 0.420 0.858
## 9 12 2 0.50 0.1118 0.323 0.775
## 10 10 1 0.45 0.1112 0.277 0.731# KM Plot overall survival
autoplot(km_overall) +
labs(title = "Overall Kaplan-Meier Survival Curve",
x = "Time (weeks)",
y = "Probability of Remaining Injury-Free") +
theme_minimal()
# KM by training load group
km_group <- survfit(Surv(time, status) ~ training_load, data = cricket)
# Plot stratified survival curves
ggsurvplot(km_group,
data = cricket,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
xlab = "Time (weeks)",
ylab = "Probability of Remaining Injury-Free",
title = "Survival Curves by Training Load",
legend.title = "Training Load",
legend.labs = c("High", "Low"))
All players started injury free (survival =1.0). The first injury occurred at week two and survival dropped to 95%, weeks 8 and 9 experienced steeper drops due to two injuries in each of those weeks and by week 10, 45% remained injury free. Remaining observations were censored after week 10.
The comparison plot shows similar curves throughout the study period for both the low and high training group. this is confirmed by a p value of 0.58, indicating no statistical significance between the two groups.
Log Rank Test - Comparing Groups
Log rank test compare the survival probability difference between two or more groups by investigating the statistical significance between the survival curves
logrank_test <- survdiff(Surv(time, status) ~ training_load,
data = cricket)
print(logrank_test)## Call:
## survdiff(formula = Surv(time, status) ~ training_load, data = cricket)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## training_load=High 10 6 5.11 0.154 0.299
## training_load=Low 10 5 5.89 0.134 0.299
##
## Chisq= 0.3 on 1 degrees of freedom, p= 0.6Understanding the output
N: Sample size per group (10 players each)
Observed: Number of injuries in each group (High = 6, Low = 5)
Expected: Expected number of injuries per group (High = 5.11, Low = 5.89
(O-E)²/E: Contribution of each group to the chi-square statistic
Chi-square = 0.3: Test statistic (measures overall difference between groups). Calculated by comparing observed vs expected values across both groups
p = 0.6: Probability that differences are due to chance, statistical significance (p <0.05)
Interpretation:
High training load group encountered 6 injuries (expected5.11) while the low load group had 5 injuries (expected 5.89). P value of 0.6 indicate no statistical significant difference between groups (p>0.0.5).
Cox Proportional Hazards Model
Cox proportional hazards regression models the relationship between the predictors and the time to event, estimating hazards ratio (HR). CPH regression is able to adjust for multiple variables and quantify their effects
Understanding Hazard Ratios:
HR = 1: No effect on injury risk
HR > 1: Increased hazard (higher injury risk)
HR < 1: Decreased hazard (lower injury risk)
# Fit Cox model with training load
cox1 <- coxph(Surv(time, status) ~ training_load, data = cricket)
summary(cox1)## Call:
## coxph(formula = Surv(time, status) ~ training_load, data = cricket)
##
## n= 20, number of events= 11
##
## coef exp(coef) se(coef) z Pr(>|z|)
## training_loadLow -0.3311 0.7182 0.6082 -0.544 0.586
##
## exp(coef) exp(-coef) lower .95 upper .95
## training_loadLow 0.7182 1.392 0.218 2.366
##
## Concordance= 0.546 (se = 0.082 )
## Likelihood ratio test= 0.3 on 1 df, p=0.6
## Wald test = 0.3 on 1 df, p=0.6
## Score (logrank) test = 0.3 on 1 df, p=0.6# Extract hazard ratio
exp(coef(cox1))## training_loadLow
## 0.7181626# Add age as a covariate
cox2 <- coxph(Surv(time, status) ~ training_load + age,
data = cricket)
summary(cox2)## Call:
## coxph(formula = Surv(time, status) ~ training_load + age, data = cricket)
##
## n= 20, number of events= 11
##
## coef exp(coef) se(coef) z Pr(>|z|)
## training_loadLow -0.9300 0.3945 0.7130 -1.304 0.19213
## age -0.3555 0.7008 0.1095 -3.247 0.00117 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## training_loadLow 0.3945 2.535 0.09754 1.5960
## age 0.7008 1.427 0.56549 0.8686
##
## Concordance= 0.816 (se = 0.054 )
## Likelihood ratio test= 12.33 on 2 df, p=0.002
## Wald test = 10.55 on 2 df, p=0.005
## Score (logrank) test = 12.74 on 2 df, p=0.002Model 1: Training Load Only
Key Results:
HR = 0.72 (95% CI: 0.22 - 2.37, p = 0.59)
Low training load has 28% lower injury hazard than high load, but this is not statistically significant
Concordance = 0.55: Poor predictive accuracy (barely better than chance)
Interpretation: Training load alone does not significantly predict injury risk in this dataset.
Model 2: Training Load + Age
Key Results:
Training Load: HR = 0.39 (95% CI: 0.10 - 1.60, p = 0.19) - Not significant
Age: HR = 0.70 (95% CI: 0.57 - 0.87, p = 0.001*) - Highly significant
Each additional year of age reduces injury hazard by 30%
Concordance = 0.82: Good predictive accuracy
Interpretation: Age is a strong predictor of injury risk. After adjusting for age, the training load effect becomes stronger but remains non-significant. The counter-intuitive age effect (older = safer) is a reminder that this is a sample dataset
Model Comparison
The model with the lowest AIC value is considered the best, as it balances goodness-of-fit with parsimony, avoiding the risk of over fitting.
Likelihood ratio test - compare how well nested models fit the data by measuring the difference in their log-likelihoods
# Compare models using AIC
AIC(cox1, cox2)## df AIC
## cox1 1 60.76971
## cox2 2 50.74019# Likelihood ratio test
anova(cox1, cox2)## Analysis of Deviance Table
## Cox model: response is Surv(time, status)
## Model 1: ~ training_load
## Model 2: ~ training_load + age
## loglik Chisq Df Pr(>|Chi|)
## 1 -29.385
## 2 -23.370 12.03 1 0.0005236 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Forest plot of hazard ratios
ggforest(cox1, data = cricket)
ggforest(cox2, data = cricket)
Which model is better?
AIC Comparison:
Model 1 (training load only): AIC = 60.77
Model 2 (training load + age): AIC = 50.74
Lower AIC = better model. Model 2 is preferred (difference > 10 indicates strong improvement).
Likelihood Ratio Test:
Chi-square = 12.03, p = 0.0005 (highly significant)
This confirms that adding age significantly improves model fit compared to training load alone.
Interpretation: Model 2 is substantially better. Age is an important predictor of injury risk and should be included when analysing training load effects.
Forest plots visualises the hazard ratios with confidence intervals.
Prediction
# Predict survival for a new player with High training load
new_player_high<- data.frame(training_load = "High")
new_player_low <- data.frame(training_load = "Low")
# Get survival probability at 12 weeks
summary(survfit(cox1, newdata = new_player_high), times = 12)## Call: survfit(formula = cox1, newdata = new_player_high)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 12 8 11 0.403 0.154 0.191 0.851summary(survfit(cox1, newdata = new_player_low), times = 12)## Call: survfit(formula = cox1, newdata = new_player_low)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 12 8 11 0.521 0.153 0.293 0.925Use of the cox model predicts the probability that a new player remains injury-free at 12 weeks, based on their training load.
Results:
High Training Load: 40.3% probability of remaining injury-free at 12 weeks (95% CI: 0.19 - 0.85)
Low Training Load: 52.1% probability of remaining injury-free at 12 weeks (95% CI: 0.29 - 0.93)
Summary Section
In this tutorial you have learnt how to understand and implement survival analysis in r. How to build and interpret models, visuals and statistic results.
Key Takeaways:
Survival analysis handles time-to-event data with censoring
Kaplan-Meier curves visualises survival probabilities over time
Log-rank tests compare survival between groups statistically using chi square statistics and p values
Cox regression models the effects of predictors on survival
Hazard ratios quantify risk (HR > 1 = increased risk, HR < 1 = decreased risk, HR = 1 = no effect)
Always check if results are statistically significant (p<0.05) and practically meaningful (effect size)