Get started with polarstate

This is a reproducible example that shows how to use polarstate, and compare the outputs to {lifelines} package.

First: Import polars as a dependency:

import polars as pl

Competing Risks Example

from polarstate import prepare_event_table

times_and_reals = 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(times_and_reals)

print(event_table)

shape: (10, 15)
┌───────┬─────────┬─────────┬─────────┬───┬──────────────┬─────────────┬─────────────┬─────────────┐
│ times ┆ count_0 ┆ count_1 ┆ count_2 ┆ … ┆ trainsition_ ┆ trainsition ┆ state_occup ┆ state_occup │
│ ---   ┆ ---     ┆ ---     ┆ ---     ┆   ┆ probabilitie ┆ _probabilit ┆ ancy_probab ┆ ancy_probab │
│ f64   ┆ i64     ┆ i64     ┆ i64     ┆   ┆ s_to_1…      ┆ ies_to_2…   ┆ ility_1_…   ┆ ility_2_…   │
│       ┆         ┆         ┆         ┆   ┆ ---          ┆ ---         ┆ ---         ┆ ---         │
│       ┆         ┆         ┆         ┆   ┆ f64          ┆ f64         ┆ f64         ┆ f64         │
╞═══════╪═════════╪═════════╪═════════╪═══╪══════════════╪═════════════╪═════════════╪═════════════╡
│ 4.3   ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.1          ┆ 0.0         ┆ 0.1         ┆ 0.0         │
│ 9.7   ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.1          ┆ 0.0         ┆ 0.2         ┆ 0.0         │
│ 14.2  ┆ 0       ┆ 0       ┆ 1       ┆ … ┆ 0.0          ┆ 0.1         ┆ 0.2         ┆ 0.1         │
│ 18.6  ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.1          ┆ 0.0         ┆ 0.3         ┆ 0.1         │
│ 24.1  ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.1          ┆ 0.0         ┆ 0.4         ┆ 0.1         │
│ 31.5  ┆ 1       ┆ 0       ┆ 0       ┆ … ┆ 0.0          ┆ 0.0         ┆ 0.4         ┆ 0.1         │
│ 34.8  ┆ 1       ┆ 0       ┆ 0       ┆ … ┆ 0.0          ┆ 0.0         ┆ 0.4         ┆ 0.1         │
│ 39.2  ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.166667     ┆ 0.0         ┆ 0.566667    ┆ 0.1         │
│ 46.0  ┆ 0       ┆ 0       ┆ 1       ┆ … ┆ 0.0          ┆ 0.166667    ┆ 0.566667    ┆ 0.266667    │
│ 49.9  ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.166667     ┆ 0.0         ┆ 0.733333    ┆ 0.266667    │
└───────┴─────────┴─────────┴─────────┴───┴──────────────┴─────────────┴─────────────┴─────────────┘

Prediction for Specific Time-Horizons

from polarstate import predict_aj_estimates

predict_aj_estimates(event_table, pl.Series([10.0, 20.0, 30.0, 40.0, 50.0]))

shape: (5, 5)

times	state_occupancy_probability_0	state_occupancy_probability_1	state_occupancy_probability_2	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	-5.5511e-17	0.733333	0.266667	"fixed_time_horizons"

import pandas as pd
from lifelines import AalenJohansenFitter

df = 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]
})

ajf_priamry = AalenJohansenFitter()
ajf_priamry.fit(df["duration"], df["event"], event_of_interest=1)

print("Event Table for event type 1")
print(ajf_priamry.cumulative_density_)

ajf_competing = AalenJohansenFitter()
ajf_competing.fit(df["duration"], df["event"], event_of_interest=2)

print("Event Table for event type 2")
print(ajf_competing.cumulative_density_)

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.266667

Prediction for Specific Time-Horizons

print("Cumulative incidence for event type 1 at times 10.0, 20.0, 30.0, 40.0, 50.0:")
print(ajf_priamry.predict([10.0, 20.0, 30.0, 40.0, 50.0]))

print("Cumulative incidence for event type 2 at times 10.0, 20.0, 30.0, 40.0, 50.0:")
print(ajf_competing.predict([10.0, 20.0, 30.0, 40.0, 50.0]))

Cumulative incidence for event type 1 at times 10.0, 20.0, 30.0, 40.0, 50.0:
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
Cumulative incidence for event type 2 at times 10.0, 20.0, 30.0, 40.0, 50.0:
10.0    0.000000
20.0    0.100000
30.0    0.100000
40.0    0.100000
50.0    0.266667
Name: CIF_2, dtype: float64

library(tidycmprsk)
cuminc(Surv(ttdeath, death_cr) ~ 1, trial)

times_and_reals <- data.frame(
  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"))
)

(cuminc(Surv(times, reals) ~ 1, times_and_reals)) |>
    tidy() |>
    print()

No Compating Risks Example

from lifelines import AalenJohansenFitter
from lifelines.datasets import load_waltons
from polarstate import prepare_event_table


T, E = load_waltons()['T'], load_waltons()['E']

times_and_reals = pl.DataFrame({
        "times": T,
        "reals": E
    })

result = prepare_event_table(times_and_reals)

print(result)

shape: (32, 15)
┌───────┬─────────┬─────────┬─────────┬───┬──────────────┬─────────────┬─────────────┬─────────────┐
│ times ┆ count_0 ┆ count_1 ┆ count_2 ┆ … ┆ trainsition_ ┆ trainsition ┆ state_occup ┆ state_occup │
│ ---   ┆ ---     ┆ ---     ┆ ---     ┆   ┆ probabilitie ┆ _probabilit ┆ ancy_probab ┆ ancy_probab │
│ f64   ┆ i64     ┆ i64     ┆ i64     ┆   ┆ s_to_1…      ┆ ies_to_2…   ┆ ility_1_…   ┆ ility_2_…   │
│       ┆         ┆         ┆         ┆   ┆ ---          ┆ ---         ┆ ---         ┆ ---         │
│       ┆         ┆         ┆         ┆   ┆ f64          ┆ f64         ┆ f64         ┆ f64         │
╞═══════╪═════════╪═════════╪═════════╪═══╪══════════════╪═════════════╪═════════════╪═════════════╡
│ 6.0   ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.006135     ┆ 0.0         ┆ 0.006135    ┆ 0.0         │
│ 7.0   ┆ 1       ┆ 1       ┆ 0       ┆ … ┆ 0.006135     ┆ 0.0         ┆ 0.01227     ┆ 0.0         │
│ 9.0   ┆ 0       ┆ 3       ┆ 0       ┆ … ┆ 0.01852      ┆ 0.0         ┆ 0.03079     ┆ 0.0         │
│ 13.0  ┆ 0       ┆ 3       ┆ 0       ┆ … ┆ 0.01852      ┆ 0.0         ┆ 0.04931     ┆ 0.0         │
│ 15.0  ┆ 0       ┆ 2       ┆ 0       ┆ … ┆ 0.012347     ┆ 0.0         ┆ 0.061656    ┆ 0.0         │
│ …     ┆ …       ┆ …       ┆ …       ┆ … ┆ …            ┆ …           ┆ …           ┆ …           │
│ 63.0  ┆ 0       ┆ 9       ┆ 0       ┆ … ┆ 0.06023      ┆ 0.0         ┆ 0.81931     ┆ 0.0         │
│ 66.0  ┆ 0       ┆ 3       ┆ 0       ┆ … ┆ 0.020077     ┆ 0.0         ┆ 0.839386    ┆ 0.0         │
│ 68.0  ┆ 1       ┆ 9       ┆ 0       ┆ … ┆ 0.06023      ┆ 0.0         ┆ 0.899616    ┆ 0.0         │
│ 69.0  ┆ 1       ┆ 12      ┆ 0       ┆ … ┆ 0.086043     ┆ 0.0         ┆ 0.985659    ┆ 0.0         │
│ 75.0  ┆ 0       ┆ 1       ┆ 0       ┆ … ┆ 0.014341     ┆ 0.0         ┆ 1.0         ┆ 0.0         │
└───────┴─────────┴─────────┴─────────┴───┴──────────────┴─────────────┴─────────────┴─────────────┘

from lifelines import AalenJohansenFitter
from lifelines.datasets import load_waltons
T, E = load_waltons()['T'], load_waltons()['E']
ajf = AalenJohansenFitter(calculate_variance=True)
ajf.fit(T, E, event_of_interest=1)
ajf.cumulative_density_

	CIF_1
event_at
0.0	0.000000
6.0	0.006135
7.0	0.012270
9.0	0.030790
13.0	0.049310
15.0	0.061656
17.0	0.067830
19.0	0.086350
22.0	0.111043
26.0	0.141910
29.0	0.172776
32.0	0.178949
33.0	0.197469
36.0	0.209816
38.0	0.222163
41.0	0.265376
43.0	0.271549
45.0	0.327109
47.0	0.333339
48.0	0.383183
51.0	0.401875
53.0	0.445488
54.0	0.457949
56.0	0.570097
58.0	0.595019
60.0	0.688476
61.0	0.745695
62.0	0.759079
63.0	0.819310
66.0	0.839386
68.0	0.899616
69.0	0.985659
75.0	1.000000