Package Comparisons
Compare native Aalen-Johansen outputs from polarstate, lifelines, and tidycmprsk.
This page uses the same competing-risks data in three implementations. As in the original polarstate documentation, each tab shows the packageβs native output, not only a normalized summary.
Shared example
Each row is one observation. Code 0 means right-censored, 1 is the event of interest, and 2 is the competing event.
times |
24.1 | 9.7 | 49.9 | 18.6 | 34.8 | 14.2 | 39.2 | 46.0 | 31.5 | 4.3 |
|---|---|---|---|---|---|---|---|---|---|---|
reals |
1 | 1 | 1 | 1 | 0 | 2 | 1 | 2 | 0 | 1 |
polarstate first returns the complete Polars event table, exposing the risk set, cause-specific hazards, survival, transition increments, and state occupancy probabilities.
import polars as pl
from polarstate import prepare_event_table, predict_aj_estimates
data = pl.DataFrame({
"times": [24.1, 9.7, 49.9, 18.6, 34.8, 14.2, 39.2, 46.0, 31.5, 4.3],
"reals": [1, 1, 1, 1, 0, 2, 1, 2, 0, 1],
})
event_table = prepare_event_table(data)
print(event_table)Native Polars output:
shape: (10, 16)
βββββββββ¬ββββββββββ¬ββββββββββ¬ββββββββββ¬ββββ¬βββββββββββββββββ¬βββββββββββββββββ¬βββββββββββββββββ¬βββββββββββββββββ
β times β count_0 β count_1 β count_2 β β¦ β trainsition_pr β state_occupanc β state_occupanc β state_occupanc β
β --- β --- β --- β --- β β obabilities_to β y_probability_ β y_probability_ β y_probability_ β
β f64 β u32 β u32 β u32 β β _2_at_times β 0_at_times β 1_at_times β 2_at_times β
β β β β β β --- β --- β --- β --- β
β β β β β β f64 β f64 β f64 β f64 β
βββββββββͺββββββββββͺββββββββββͺββββββββββͺββββͺβββββββββββββββββͺβββββββββββββββββͺβββββββββββββββββͺβββββββββββββββββ‘
β 4.3 β 0 β 1 β 0 β β¦ β 0.0 β 0.9 β 0.1 β 0.0 β
β 9.7 β 0 β 1 β 0 β β¦ β 0.0 β 0.8 β 0.2 β 0.0 β
β 14.2 β 0 β 0 β 1 β β¦ β 0.1 β 0.7 β 0.2 β 0.1 β
β 18.6 β 0 β 1 β 0 β β¦ β 0.0 β 0.6 β 0.3 β 0.1 β
β 24.1 β 0 β 1 β 0 β β¦ β 0.0 β 0.5 β 0.4 β 0.1 β
β 31.5 β 1 β 0 β 0 β β¦ β 0.0 β 0.5 β 0.4 β 0.1 β
β 34.8 β 1 β 0 β 0 β β¦ β 0.0 β 0.5 β 0.4 β 0.1 β
β 39.2 β 0 β 1 β 0 β β¦ β 0.0 β 0.333333 β 0.566667 β 0.1 β
β 46.0 β 0 β 0 β 1 β β¦ β 0.166667 β 0.166667 β 0.566667 β 0.266667 β
β 49.9 β 0 β 1 β 0 β β¦ β 0.0 β 0.0 β 0.733333 β 0.266667 β
βββββββββ΄ββββββββββ΄ββββββββββ΄ββββββββββ΄ββββ΄βββββββββββββββββ΄βββββββββββββββββ΄βββββββββββββββββ΄βββββββββββββββββPrediction for specific time horizons
estimates = predict_aj_estimates(
event_table,
pl.Series([10.0, 20.0, 30.0, 40.0, 50.0]),
)
print(estimates)Native Polars output:
shape: (5, 5)
βββββββββ¬βββββββββββββββββββββββββ¬βββββββββββββββββββββββββ¬βββββββββββββββββββββββββ¬ββββββββββββββββββββββ
β times β state_occupancy_probβ¦ β state_occupancy_probβ¦ β state_occupancy_probβ¦ β estimate_origin β
β --- β --- β --- β --- β --- β
β f64 β f64 β f64 β f64 β enum β
βββββββββͺβββββββββββββββββββββββββͺβββββββββββββββββββββββββͺβββββββββββββββββββββββββͺββββββββββββββββββββββ‘
β 10.0 β 0.8 β 0.2 β 0.0 β fixed_time_horizons β
β 20.0 β 0.6 β 0.3 β 0.1 β fixed_time_horizons β
β 30.0 β 0.5 β 0.4 β 0.1 β fixed_time_horizons β
β 40.0 β 0.333333 β 0.566667 β 0.1 β fixed_time_horizons β
β 50.0 β 0.0 β 0.733333 β 0.266667 β fixed_time_horizons β
βββββββββ΄βββββββββββββββββββββββββ΄βββββββββββββββββββββββββ΄βββββββββββββββββββββββββ΄ββββββββββββββββββββββlifelines fits and returns one cumulative-incidence table for each event type.
import pandas as pd
from lifelines import AalenJohansenFitter
data = pd.DataFrame({
"duration": [24.1, 9.7, 49.9, 18.6, 34.8, 14.2, 39.2, 46.0, 31.5, 4.3],
"event": [1, 1, 1, 1, 0, 2, 1, 2, 0, 1],
})
primary = AalenJohansenFitter(calculate_variance=False)
primary.fit(data["duration"], data["event"], event_of_interest=1)
print("Event table for event type 1")
print(primary.cumulative_density_)
competing = AalenJohansenFitter(calculate_variance=False)
competing.fit(data["duration"], data["event"], event_of_interest=2)
print("Event table for event type 2")
print(competing.cumulative_density_)Native pandas output:
Event table for event type 1
CIF_1
event_at
0.0 0.000000
4.3 0.100000
9.7 0.200000
14.2 0.200000
18.6 0.300000
24.1 0.400000
31.5 0.400000
34.8 0.400000
39.2 0.566667
46.0 0.566667
49.9 0.733333
Event table for event type 2
CIF_2
event_at
0.0 0.000000
4.3 0.000000
9.7 0.000000
14.2 0.100000
18.6 0.100000
24.1 0.100000
31.5 0.100000
34.8 0.100000
39.2 0.100000
46.0 0.266667
49.9 0.266667Prediction for specific time horizons
print(primary.predict([10.0, 20.0, 30.0, 40.0, 50.0]))
print(competing.predict([10.0, 20.0, 30.0, 40.0, 50.0]))Native pandas output:
10.0 0.200000
20.0 0.300000
30.0 0.400000
40.0 0.566667
50.0 0.733333
Name: CIF_1, dtype: float64
10.0 0.000000
20.0 0.100000
30.0 0.100000
40.0 0.100000
50.0 0.266667
Name: CIF_2, dtype: float64tidycmprsk returns outcome-specific rows in a tibble.
library(dplyr)
library(tidycmprsk)
data <- tibble(
times = c(24.1, 9.7, 49.9, 18.6, 34.8, 14.2, 39.2, 46.0, 31.5, 4.3),
reals = factor(
c(1, 1, 1, 1, 0, 2, 1, 2, 0, 1),
levels = c(0, 1, 2),
labels = c("censored", "primary_event", "competing_event")
)
)
fit <- cuminc(Surv(times, reals) ~ 1, data)
tidy(fit) |>
select(outcome, time, estimate) |>
print(n = 20)Native tibble output:
# A tibble: 20 Γ 3
outcome time estimate
<chr> <dbl> <dbl>
1 primary_event 4.3 0.100
2 primary_event 9.7 0.200
3 primary_event 14.2 0.200
4 primary_event 18.6 0.300
5 primary_event 24.1 0.400
6 primary_event 31.5 0.400
7 primary_event 34.8 0.400
8 primary_event 39.2 0.567
9 primary_event 46.0 0.567
10 primary_event 49.9 0.733
11 competing_event 4.3 0
12 competing_event 9.7 0
13 competing_event 14.2 0.100
14 competing_event 18.6 0.100
15 competing_event 24.1 0.100
16 competing_event 31.5 0.100
17 competing_event 34.8 0.100
18 competing_event 39.2 0.100
19 competing_event 46.0 0.267
20 competing_event 49.9 0.267