diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 9e8db0b..8cfc871 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-11-25T10:55:50","documenter_version":"1.8.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-12-02T14:28:27","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/benches/index.html b/dev/benches/index.html index 8508e27..2166f94 100644 --- a/dev/benches/index.html +++ b/dev/benches/index.html @@ -51,4 +51,4 @@ # using PrettyTables # pretty_table(rez, backend = Val(:markdown)) -rez
5×3 DataFrame
RowAlgorithmunstratifiedstratified
StringFloat64Float64
1Pohar Perme5.744234.97847
2EdererI2.755761.65064
3EdererII12.501210.5824
4Hakulinen8.516235.77065
5Graffeo's LRT9.0047112.5814

Benchmarking across time

The following charts provide a glimpse of NetSurvival.jl's performance along time, also ran on github CI:

+rez
5×3 DataFrame
RowAlgorithmunstratifiedstratified
StringFloat64Float64
1Pohar Perme6.03315.23631
2EdererI2.776221.66969
3EdererII12.949110.3885
4Hakulinen8.677965.68397
5Graffeo's LRT9.1963813.3658

Benchmarking across time

The following charts provide a glimpse of NetSurvival.jl's performance along time, also ran on github CI:

diff --git a/dev/example/9aecfb2a.svg b/dev/example/0020f193.svg similarity index 93% rename from dev/example/9aecfb2a.svg rename to dev/example/0020f193.svg index 827c957..657684c 100644 --- a/dev/example/9aecfb2a.svg +++ b/dev/example/0020f193.svg @@ -1,43 +1,43 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/example/17b99c97.svg b/dev/example/17b99c97.svg deleted file mode 100644 index 5078fc5..0000000 --- a/dev/example/17b99c97.svg +++ /dev/null @@ -1,78 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/example/362ecc30.svg b/dev/example/362ecc30.svg new file mode 100644 index 0000000..58da4fa --- /dev/null +++ b/dev/example/362ecc30.svg @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/example/3e1ca2fa.svg b/dev/example/3e1ca2fa.svg deleted file mode 100644 index 940dc01..0000000 --- a/dev/example/3e1ca2fa.svg +++ /dev/null @@ -1,90 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/example/76fdcb1c.svg b/dev/example/8bb03385.svg similarity index 77% rename from dev/example/76fdcb1c.svg rename to dev/example/8bb03385.svg index 0e06acf..1238179 100644 --- a/dev/example/76fdcb1c.svg +++ b/dev/example/8bb03385.svg @@ -1,247 +1,247 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/example/c486d20b.svg b/dev/example/c486d20b.svg new file mode 100644 index 0000000..34c5e34 --- /dev/null +++ b/dev/example/c486d20b.svg @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/example/e26f4cf4.svg b/dev/example/e26f4cf4.svg new file mode 100644 index 0000000..0e61623 --- /dev/null +++ b/dev/example/e26f4cf4.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/example/index.html b/dev/example/index.html index 4eac5c8..a0dd9b2 100644 --- a/dev/example/index.html +++ b/dev/example/index.html @@ -1,22 +1,15 @@ -Case study · NetSurvival.jl
GitHub

Case study

Introduction and datasets

We will illustrate with an example using the dataset colrec, which comprises $5971$ patients diagnosed with colon or rectal cancer between $1994$ and $2000$. This dataset is sourced from the Slovenia cancer registry. Given the high probability that the patients are Slovenian, we will be using the Slovenian mortality table slopop as reference for the populational rates. Subsequently, we can apply various non-parametric estimators for net survival analysis.

N.B.

Mortality tables may vary in structure, with options such as the addition or removal of specific covariates. To confirm that the mortality table is in the correct format, please refer to RateTables.jl's documentation, or directly extract it from there.

Cohort details

The patients in the study are diagnosed between January 1st 1994 and December 31st 2000. Before we move on to the survival probabilities, it is important to be aware of how your data is distributed and of what it comprises. 

using NetSurvival, RateTables, DataFrames
-
+Case study · NetSurvival.jl
GitHub
Case study
Introduction and datasets
We will illustrate with an example using the dataset colrec, which comprises $5971$ patients diagnosed with colon or rectal cancer between $1994$ and $2000$. This dataset is sourced from the Slovenia cancer registry. Given the high probability that the patients are Slovenian, we will be using the Slovenian mortality table slopop as reference for the population mortality rates. Subsequently, we can apply various non-parametric estimators for net survival analysis.
N.B.
Mortality tables may vary in structure, with options such as the addition or removal of specific covariates. To confirm that the mortality table is in the correct format, please refer to RateTables.jl's documentation, or directly extract it from there.
Cohort details
The patients in the study are diagnosed between January 1st $1994$ and December 31st $2000$. Before we move on to the survival probabilities, it is important to be aware of how your data is distributed and of what it comprises. 
using NetSurvival, RateTables, DataFrames
 first(colrec,10)
10×7 DataFrame
Rowageyearsextimestatusstagesite
Float64Float64SymbolFloat64BoolInt64String7
123004.0728528.0male16.0false1rectum
212082.0729260.0female504.0false3rectum
324277.0728583.0male22.0false3colon
429256.0729843.0female3998.0false1colon
530416.0728869.0female9.0false99colon
615156.0729686.0male88.0false2colon
723953.0728713.0female7769.0false1rectum
820775.0728867.0male7366.0false1rectum
917555.0730141.0female6480.0false1rectum
1020329.0730137.0female6449.0false1rectum
Let's explore how our data is distributed first, starting with the age variable.
println(minimum(colrec.age./365.241), " ; ",maximum(colrec.age./365.241))
12.482169307388821 ; 96.7169622249419
The study can be considered diverse in terms of age seeing as the patients are between 12 and 96 years old, approximately. 
using Plots
-
-plot(
-    histogram(colrec.age./365.241, label="Age"),
-    histogram(colrec.time./365.241, label="Follow-up time")
-)
Example block output
The graph above show us that although it has a wide range of patients within all age groups, it is mostly centered around older adults and elderly, with the majority of the patients being between 60 and 80 years old. 
Looking at the second graph that details the distribution of the follow-up times, we notice that the values quickly drop. Unfortunately, this is a common theme in cancer studies. 
Let's take a look at the sex variable now: 
combine(groupby(colrec, :sex), nrow)
2×2 DataFrame
Rowsexnrow
SymbolInt64
1male3289
2female2682
This dataframe shows the number of male and female patients. There isn't too big of a difference between the two. We can say this study includes both gender relatively equally, thus, reducing bias. 
With these two observations, it is also worth noting that colorectal cancer is most common with men and people older than 50.
In total, we note that we have $5971$ patients. By taking a look at the status variable, we can determine the deaths and censorship:
sum(colrec.status)
4979
Out of the $5971$ patients, $4979$ have died. This is a very high mortality rate, and again, unfortunately common in cancer studies.
(nrow(colrec) - sum(colrec.status)) / nrow(colrec)
0.16613632557360575
In other terms, the censorship rate is of 16.6%, meaning the event, in this case death, was not observed for only 16.6% of the patients. This is a low censorship rate and thus the quality of the signal will be pretty good. 
Mortality table
We will be using the mortality table slopop taken from the RateTables.jl package as the study is done on Slovenian patients. 
slopop
RateTable(:sex,)
By examining slopop, we notice it contains information regarding age and year, as expected for mortality tables. Additionally, it incorporates the covariate sex, which has two possible entries (:male or :female). The ratetable is then three dimensional, with the covariate sex added. For example, the daily hazard rate for a woman turning $45$ on the January 1st $2006$ can be accessed through the following command:
daily_hazard(slopop, 45*365.241, 2006*365.241; sex=:female)
2.02680951796113e-6
Making the survival probability easily calculated with:
exp(-(daily_hazard(slopop, 45*365.241, 2006*365.241; sex=:female))*365)
0.9992604880997519
Overall and expected survival
For this part, we will be using the Survival.jl package to apply the Kaplan Meier estimator for the overall survival.
using Survival
-
+plot(histogram(colrec.age./365.241, label="Age"),
+    histogram(colrec.time./365.241, label="Follow-up time"))
Example block output
The graph above show us that although the dataset has a wide range of patients within all age groups, it is mostly centered around older adults and elderly, with the majority of the patients being between $60$ and $80$ years old. Looking at the second graph that details the distribution of the follow-up times. We notice there that the values quickly drop. Unfortunately, this is a common theme in cancer studies. 
Let's take a look at the sex variable now, by looking at the number of male and female patients:
combine(groupby(colrec, :sex), nrow)
2×2 DataFrame
Rowsexnrow
SymbolInt64
1male3289
2female2682
There isn't too big of a difference between the two. We can say this study includes both gender relatively equally, thus, reducing bias. With these two observations, it is also worth noting that colorectal cancer is most common with men and people older than $50$.
In total, we note that we have $5971$ patients. By taking a look at the status variable, we can determine the deaths and censorship:
sum(colrec.status)
4979
Out of the $5971$ patients, $4979$ have died. This is a very high mortality rate, and again, unfortunately common in cancer studies.
(nrow(colrec) - sum(colrec.status)) / nrow(colrec)
0.16613632557360575
In other terms, the censorship rate is of 16.6%, meaning the event, in this case death, was not observed for only 16.6% of the patients. This is a low censorship rate and thus the quality of the signal will be pretty good. 
Mortality table
We will be using the mortality table slopop taken from the RateTables.jl package as the study is done on Slovenian patients. 
slopop
RateTable(:sex,)
The show method of the RateTable class shows the additional covariate sex that the rate table has on top of the (mandatory) age and year variables. The sex variable has two madalities, :male and :female. The ratetable is then three dimensional. For example, the daily hazard rate for a woman turning $45$ on the January 1st $2006$ can be accessed through the following command:
λ  = daily_hazard(slopop, 45*365.241, 2006*365.241; sex=:female)
2.02680951796113e-6
Making the daily survival probability easily calculated with:
exp(-λ)
0.999997973192536
Overall and expected survival
For this part, we will be using the Survival.jl package to apply the Kaplan Meier estimator for the overall survival.
using Survival
 km = fit(KaplanMeier, colrec.time./365.241, colrec.status)
-plot(km.events.time, km.survival, label=false, title = "Kaplan-Meier Estimator for the Overall Survival")
Example block output
The graph above indicates a significant dip in survival probability within the first $5$ years, and even more at $10$ years. This period in the study is then crucial for the analysis. 
Estimated net survival
In this part, we are interested in the first $5$ years of the study. We will thus limit the follow-up time to $5$ years, meaning we will censor all individuals with a follow-up time that is higher than this. Then, we will apply the different net survival methods.
colrec.time5 .= 0.0
-colrec.status5 .= Bool(true)
-for i in 1:nrow(colrec)
-    colrec.time5[i] = min(colrec.time[i], round(5*365.241))
-    if colrec.time[i] > 5*365.241
-        colrec.status5[i] = false
+plot(km.events.time, km.survival, label=false, title = "Kaplan-Meier Estimator for the Overall Survival")
Example block output
The graph above indicates a significant dip in survival probability within the first $5$ years, and even more at $10$ years. This period in the study is then crucial for the analysis. 
Estimated net survival
We will restrict ourselves to the first $5$ years of the study. For that, let us re-censor the dataset as follows: 
for i in 1:nrow(colrec)
+    if colrec.time[i] > 1826 # five years
+        colrec.status[i] = false
+        colrec.time[i] = 1826
     end
-end
Now that we have defined our own time and status variables according to the observations made, we can apply the different non parametric methods for relative survival.
e1 = fit(EdererI, @formula(Surv(time5,status5)~1), colrec, slopop)
EdererI(t ∈ (1.0, 1826.0)) with summary stats:
+end
We can now apply the different non-parametric methods to compute the relative survival.
e1 = fit(EdererI, @formula(Surv(time,status)~1), colrec, slopop)
EdererI(t ∈ (1.0, 1826.0)) with summary stats:
  1826×5 DataFrame
   Row │ Sₑ        ∂Λₑ           σₑ           lower_95_CI  upper_95_CI 
       │ Float64   Float64       Float64      Float64      Float64     
@@ -30,15 +23,15 @@
     7 │ 0.981734   0.00158443   0.00180287      0.978271     0.98521
     8 │ 0.976991   0.00483135   0.00202379      0.973124     0.980874
   ⋮   │    ⋮           ⋮             ⋮            ⋮            ⋮
- 1820 │ 0.456494  -0.000133259  0.017137        0.441416     0.472087
- 1821 │ 0.456555  -0.000133262  0.017137        0.441475     0.47215
- 1822 │ 0.456616  -0.000133307  0.017137        0.441534     0.472213
- 1823 │ 0.456677  -0.000133276  0.017137        0.441593     0.472276
- 1824 │ 0.456527   0.000328197  0.0171432       0.441442     0.472127
- 1825 │ 0.456377   0.000328384  0.0171494       0.441292     0.471977
- 1826 │ 0.456016   0.000790464  0.0171618       0.440932     0.471616
+ 1820 │ 0.456957  -0.000133259  0.0171307       0.441869     0.47256
+ 1821 │ 0.457018  -0.000133262  0.0171307       0.441928     0.472623
+ 1822 │ 0.457079  -0.000133307  0.0171307       0.441987     0.472686
+ 1823 │ 0.457139  -0.000133276  0.0171307       0.442045     0.472749
+ 1824 │ 0.456989   0.000328197  0.017137        0.441895     0.472599
+ 1825 │ 0.456839   0.000328384  0.0171432       0.441745     0.47245
+ 1826 │ 0.456478   0.000790464  0.0171556       0.441385     0.472088
                                                      1811 rows omitted
With the EdererI method, after $1826$ days have passed, we can say that the survival rate at this mark is around $0.456$, in the hypothetical world where patients can only die of cancer.
crude_e1 = CrudeMortality(e1)
-println(crude_e1.Mₒ[1826], " , ", crude_e1.Mₑ[1826], " , ", crude_e1.Mₚ[1826])
0.6368404579125496 , 0.5158124650476976 , 0.121027992864852
Out of the 0.63 patients that have died, according to the EdererI method, 0.51 died because of colorectal cancer and 0.12 died of other causes.
e2 = fit(EdererII, @formula(Surv(time5,status5)~1), colrec, slopop)
EdererII(t ∈ (1.0, 1826.0)) with summary stats:
+println(crude_e1.Mₒ[1826], " , ", crude_e1.Mₑ[1826], " , ", crude_e1.Mₚ[1826])
0.636476800600479 , 0.515341453104102 , 0.121135347496377
Out of the 0.63 patients that have died, according to the EdererI method, 0.51 died because of colorectal cancer and 0.12 died of other causes.
e2 = fit(EdererII, @formula(Surv(time,status)~1), colrec, slopop)
EdererII(t ∈ (1.0, 1826.0)) with summary stats:
  1826×5 DataFrame
   Row │ Sₑ        ∂Λₑ           σₑ           lower_95_CI  upper_95_CI 
       │ Float64   Float64       Float64      Float64      Float64     
@@ -52,15 +45,15 @@
     7 │ 0.981728   0.00158595   0.00180287      0.978265     0.985203
     8 │ 0.976983   0.00483311   0.00202379      0.973116     0.980866
   ⋮   │    ⋮           ⋮             ⋮            ⋮            ⋮
- 1820 │ 0.441142  -0.00011516   0.017137        0.426571     0.45621
- 1821 │ 0.441193  -0.000115156  0.017137        0.42662      0.456263
- 1822 │ 0.441243  -0.000115176  0.017137        0.426669     0.456316
- 1823 │ 0.441294  -0.000115213  0.017137        0.426718     0.456368
- 1824 │ 0.441142   0.00034624   0.0171432       0.426565     0.456216
- 1825 │ 0.440989   0.000346509  0.0171494       0.426412     0.456063
- 1826 │ 0.440632   0.000808582  0.0171618       0.426057     0.455706
-                                                     1811 rows omitted
Similarily, the EdererII method, also known as the conditional method, shows that at the $5$ year mark, the survival probability is of $0.44$ in this hypothetical world.
crude_e2 = CrudeMortality(e2)
-println(crude_e2.Mₒ[1826], " , ", crude_e2.Mₑ[1826], " , ", crude_e2.Mₚ[1826])
0.6368404579125496 , 0.5334518856851325 , 0.10338857222741712
Here, out of the 0.63 patients that have dued, 0.53 are due to colorectal cancer and 0.1 due to other causes.
pp = fit(PoharPerme, @formula(Surv(time5,status5)~1), colrec, slopop)
PoharPerme(t ∈ (1.0, 1826.0)) with summary stats:
+ 1820 │ 0.441589  -0.00011516   0.0171307       0.427008     0.456667
+ 1821 │ 0.44164   -0.000115156  0.0171307       0.427058     0.45672
+ 1822 │ 0.441691  -0.000115176  0.0171307       0.427107     0.456772
+ 1823 │ 0.441741  -0.000115213  0.0171307       0.427156     0.456825
+ 1824 │ 0.441588   0.00034624   0.017137        0.427003     0.456672
+ 1825 │ 0.441435   0.000346509  0.0171432       0.42685      0.45652
+ 1826 │ 0.441079   0.000808582  0.0171556       0.426494     0.456162
+                                                     1811 rows omitted
Similarly, the EdererII method, also known as the conditional method, shows that at the $5$ year mark, the survival probability is of $0.44$ in this hypothetical world.
crude_e2 = CrudeMortality(e2)
+println(crude_e2.Mₒ[1826], " , ", crude_e2.Mₑ[1826], " , ", crude_e2.Mₚ[1826])
0.636476800600479 , 0.5329969626451269 , 0.10347983795535214
Here, out of the 0.63 patients that have died, 0.53 are due to colorectal cancer and 0.1 due to other causes.
pp = fit(PoharPerme, @formula(Surv(time,status)~1), colrec, slopop)
PoharPerme(t ∈ (1.0, 1826.0)) with summary stats:
  1826×5 DataFrame
   Row │ Sₑ        ∂Λₑ           σₑ           lower_95_CI  upper_95_CI 
       │ Float64   Float64       Float64      Float64      Float64     
@@ -74,24 +67,28 @@
     7 │ 0.981722   0.00158763   0.0018035       0.978258     0.985199
     8 │ 0.976973   0.00483737   0.00202472      0.973104     0.980858
   ⋮   │    ⋮           ⋮             ⋮            ⋮            ⋮
- 1820 │ 0.44143   -0.000148858  0.0178917       0.426219     0.457184
- 1821 │ 0.441496  -0.000148885  0.0178917       0.426282     0.457252
- 1822 │ 0.441561  -0.000148937  0.0178917       0.426345     0.45732
- 1823 │ 0.441627  -0.000149005  0.0178917       0.426409     0.457389
- 1824 │ 0.441403   0.00050879   0.0179038       0.426182     0.457167
- 1825 │ 0.441281   0.000276025  0.0179089       0.42606      0.457045
- 1826 │ 0.440922   0.000813485  0.0179221       0.425702     0.456685
-                                                     1811 rows omitted
We conclude for the Poher-Perme method, that in a world where cancer patients could only die due to cancer, only 41% of these patients would still be alive $5$ year after their diagnosis.
crude_pp = CrudeMortality(pp)
-println(crude_pp.Mₒ[1826], " , ", crude_pp.Mₑ[1826], " , ", crude_pp.Mₚ[1826])
0.6458893614663985 , 0.5319347797905097 , 0.11395458167588875
Finally, for this estimator, we have that of the 0.64 patients that have died, 0.53 is due to colorectal cancer while 0.11 is due to other causes.
We will plot the Pohar Perme method only.
conf_int = confint(pp; level = 0.05)
-lower_bounds = [lower[1] for lower in conf_int]
-upper_bounds = [upper[2] for upper in conf_int]
+ 1820 │ 0.44187   -0.000148858  0.017886        0.426649     0.457635
+ 1821 │ 0.441936  -0.000148885  0.017886        0.426712     0.457703
+ 1822 │ 0.442002  -0.000148937  0.017886        0.426776     0.457772
+ 1823 │ 0.442068  -0.000149005  0.017886        0.426839     0.45784
+ 1824 │ 0.441843   0.00050879   0.0178981       0.426612     0.457618
+ 1825 │ 0.441721   0.000276025  0.0179032       0.42649      0.457496
+ 1826 │ 0.441362   0.000813485  0.0179164       0.426132     0.457136
+                                                     1811 rows omitted
We conclude for the Pohar Perme method, that in a world where cancer patients could only die due to cancer, only 41% of these patients would still be alive $5$ year after their diagnosis. The Pohar Perme estimator is the best estimator of the excess hazard under the standard hypotheses. 
crude_pp = CrudeMortality(pp)
+println(crude_pp.Mₒ[1826], " , ", crude_pp.Mₑ[1826], " , ", crude_pp.Mₚ[1826])
0.6455403791714492 , 0.5314853737936457 , 0.1140550053778036
Finally, for this estimator, we have that of the 0.64 patients that have died, 0.53 is due to colorectal cancer while 0.11 is due to other causes.
We will plot the Pohar Perme method only.
function mkribbon(pp)
+    S = pp.Sₑ
+    ci = confint(pp; level = 0.05)
+    l,u = getindex.(ci, 1), getindex.(ci,2)
+    rb = (S - l, u - S)
+    return rb
+end
 
-p1 = plot(pp.grid, pp.Sₑ, ribbon=(pp.Sₑ - lower_bounds, upper_bounds - pp.Sₑ), xlab = "Time (days)", ylab = "Net survival", label = false)
+p1 = plot(pp.grid, pp.Sₑ, ribbon=mkribbon(pp), xlab = "Time (days)", ylab = "Net survival", label = false)
 
 p2 = plot(pp.grid, crude_pp.Mₑ, label = "Excess Mortality Rate")
 p2 = plot!(pp.grid, crude_pp.Mₚ, label = "Population Mortality Rate")
 
-plot(p1,p2)
Example block output
Looking at the graph, and the rapid dip it takes, it is evident that the first $5$ years are crucial and that the survival probability is highly affected in these years. Additionnally, the crude mortality graph allows us to see how much of this curve is due to the colorectacl cancer studied versus other undefined causes. It is clear that the large majority is due to the cancer.
Net survival with respect to covariates
We are now interested in comparing the different groups of patients defined by various covariates. 
pp_sex = fit(PoharPerme, @formula(Surv(time5,status5)~sex), colrec, slopop)
+plot(p1,p2)
Example block output
Looking at the graph, and the rapid dip it takes, it is evident that the first $5$ years are crucial and that the survival probability is highly affected in these years. Additionally, the crude mortality graph allows us to see how much of this curve is due to the colorectal cancer studied versus other undefined causes. It is clear that the large majority is due to the cancer.
Net survival with respect to covariates
We are now interested in comparing the different groups of patients defined by various covariates. 
pp_sex = fit(PoharPerme, @formula(Surv(time,status)~sex), colrec, slopop)
 pp_males = pp_sex[pp_sex.sex .== :male,:estimator][1]
 pp_females = pp_sex[pp_sex.sex .== :female,:estimator][1]
PoharPerme(t ∈ (1.0, 1826.0)) with summary stats:
  1826×5 DataFrame
@@ -107,20 +104,20 @@
     7 │ 0.97871    0.00142212   0.00289239     0.973177     0.984274
     8 │ 0.97321    0.00561896   0.00324789     0.967035     0.979425
   ⋮   │    ⋮           ⋮            ⋮            ⋮            ⋮
- 1820 │ 0.450196  -0.000129968  0.025504       0.428245     0.473272
- 1821 │ 0.450254  -0.000129989  0.025504       0.428301     0.473333
- 1822 │ 0.450313  -0.000130047  0.025504       0.428356     0.473395
- 1823 │ 0.450371  -0.000130101  0.025504       0.428412     0.473456
- 1824 │ 0.45043   -0.000130144  0.025504       0.428468     0.473518
- 1825 │ 0.450071   0.000797147  0.0255209      0.428112     0.473156
- 1826 │ 0.449584   0.00108151   0.0255496      0.427625     0.472671
-                                                    1811 rows omitted
When comparing at time $1826$, we notice that the survival probability is slightly inferior for men than for women ($0.433 < 0.449$). It is also more probable for the women to die from other causes than the men seeing as $0.0255 > 0.025$. Still, the differences are minimal. Let's confirm this with the Grafféo log-rank test:
test_sex = fit(GraffeoTest, @formula(Surv(time5,status5)~sex), colrec, slopop)
Grafféo's log-rank-type-test
+ 1820 │ 0.450638  -0.000129968  0.0254941      0.428674     0.473727
+ 1821 │ 0.450697  -0.000129989  0.0254941      0.42873      0.473789
+ 1822 │ 0.450755  -0.000130047  0.0254941      0.428786     0.473851
+ 1823 │ 0.450814  -0.000130101  0.0254941      0.428841     0.473912
+ 1824 │ 0.450873  -0.000130144  0.0254941      0.428897     0.473974
+ 1825 │ 0.450513   0.000797147  0.025511       0.428541     0.473612
+ 1826 │ 0.450026   0.00108151   0.0255398      0.428054     0.473126
+                                                    1811 rows omitted
When comparing at time $1826$, we notice that the survival probability is slightly inferior for men than for women ($0.433 < 0.449$). It is also more probable for the women to die from other causes than the men seeing as $0.0255 > 0.025$. Still, the differences are minimal. Let's confirm this with the Grafféo log-rank test:
test_sex = fit(GraffeoTest, @formula(Surv(time,status)~sex), colrec, slopop)
Grafféo's log-rank-type-test
 1×3 DataFrame
  Row │ test_statistic  degrees_of_freedom  p_value  
      │ Float64         Int64               Float64  
 ─────┼──────────────────────────────────────────────
-   1 │        0.54326                   1  0.461085
The p-value is indeed above $0.05$. We cannot reject the null hypothesis $H_0$ and thus we dismiss the differences between the two sexes.
As for the age, we will define two different groups: individuals aged 65 and above and those who are not.
colrec.age65 .= ifelse.(colrec.age .>= 65*365.241, :old, :young)
-pp_age65 = fit(PoharPerme, @formula(Surv(time5,status5)~age65), colrec, slopop)
+   1 │       0.533424                   1  0.465171
The p-value is indeed above $0.05$. We cannot reject the null hypothesis $H_0$ and thus we dismiss the differences between the two sexes.
As for the age, we will define two different groups: individuals aged 65 and above and those who are not.
colrec.age65 .= ifelse.(colrec.age .>= 65*365.241, :old, :young)
+pp_age65 = fit(PoharPerme, @formula(Surv(time,status)~age65), colrec, slopop)
 pp_young = pp_age65[pp_age65.age65 .== :young, :estimator][1]
 pp_old = pp_age65[pp_age65.age65 .== :old, :estimator][1]
PoharPerme(t ∈ (1.0, 1826.0)) with summary stats:
  1826×5 DataFrame
@@ -136,59 +133,34 @@
     7 │ 0.974672   0.00180987   0.00273909     0.969454     0.979919
     8 │ 0.967646   0.00720934   0.0030984      0.961787     0.97354
   ⋮   │    ⋮           ⋮            ⋮            ⋮            ⋮
- 1820 │ 0.40224   -0.000235167  0.0269335      0.381557     0.424044
- 1821 │ 0.402335  -0.0002352    0.0269335      0.381647     0.424144
- 1822 │ 0.402429  -0.000235279  0.0269335      0.381737     0.424244
- 1823 │ 0.402524  -0.000235379  0.0269335      0.381826     0.424344
- 1824 │ 0.402139   0.000955393  0.0269598      0.381442     0.42396
- 1825 │ 0.401925   0.000534307  0.0269708      0.38123      0.423743
- 1826 │ 0.401615   0.000770012  0.0269896      0.380922     0.423432
-                                                    1811 rows omitted
Here, the difference between the two is much more important. In the first group, the individuals are aged under 65 and at $5$ years time, they have a $50.1$% chance of survival. On the other hand, the individuals aged 65 and up have a $40.1$% chance of survival. 
It is also worth noting that their chances of dying from other causes is higher than the younger group, given their age. 
When applying the Grafféo test, we get the results below:
test_age65 = fit(GraffeoTest, @formula(Surv(time5,status5)~age65), colrec, slopop)
Grafféo's log-rank-type-test
+ 1820 │ 0.402479  -0.000235167  0.0269303      0.381786     0.424294
+ 1821 │ 0.402574  -0.0002352    0.0269303      0.381876     0.424393
+ 1822 │ 0.402669  -0.000235279  0.0269303      0.381966     0.424493
+ 1823 │ 0.402763  -0.000235379  0.0269303      0.382056     0.424593
+ 1824 │ 0.402379   0.000955393  0.0269566      0.381671     0.424209
+ 1825 │ 0.402164   0.000534307  0.0269676      0.381459     0.423992
+ 1826 │ 0.401854   0.000770012  0.0269863      0.381151     0.423681
+                                                    1811 rows omitted
Here, the difference between the two is much more important. In the first group, the individuals are aged under 65 and at $5$ years time, they have a $50.1$% chance of survival. On the other hand, the individuals aged 65 and up have a $40.1$% chance of survival. 
It is also worth noting that their chances of dying from other causes is higher than the younger group, given their age. 
When applying the Grafféo test, we get the results below:
test_age65 = fit(GraffeoTest, @formula(Surv(time,status)~age65), colrec, slopop)
Grafféo's log-rank-type-test
 1×3 DataFrame
  Row │ test_statistic  degrees_of_freedom  p_value     
      │ Float64         Int64               Float64     
 ─────┼─────────────────────────────────────────────────
-   1 │        76.8388                   1  1.85495e-18
The p-value is well under $0.05$, meaning we reject the $H_0$ hypothesis and must admit there are differences between the individuals aged 65 and above and the others.
When plotting both we get:
conf_int_men = confint(pp_males; level = 0.05)
-lower_bounds_men = [lower[1] for lower in conf_int_men]
-upper_bounds_men = [upper[2] for upper in conf_int_men]
-
-conf_int_women = confint(pp_females; level = 0.05)
-lower_bounds_women = [lower[1] for lower in conf_int_women]
-upper_bounds_women = [upper[2] for upper in conf_int_women]
-
-conf_int_under65 = confint(pp_young; level = 0.05)
-lower_bounds_under65 = [lower[1] for lower in conf_int_under65]
-upper_bounds_under65 = [upper[2] for upper in conf_int_under65]
-
-conf_int_65 = confint(pp_old; level = 0.05)
-lower_bounds_65 = [lower[1] for lower in conf_int_65]
-upper_bounds_65 = [upper[2] for upper in conf_int_65]
-
-plot1 = plot(pp_males.grid, pp_males.Sₑ, ribbon=(pp_males.Sₑ - lower_bounds_men, upper_bounds_men - pp_males.Sₑ), xlab = "Time (days)", ylab = "Net survival", label = "men")
-plot1 = plot!(pp_females.grid, pp_females.Sₑ, ribbon=(pp_females.Sₑ - lower_bounds_women, upper_bounds_women - pp_females.Sₑ), xlab = "Time (days)", ylab = "Net survival", label = "women")
+   1 │        77.4381                   1  1.36951e-18
The p-value is well under $0.05$, meaning we reject the $H_0$ hypothesis and must admit there are differences between the individuals aged 65 and above and the others.
When plotting both we get:
plot1 = plot(pp_males.grid, pp_males.Sₑ, ribbon=mkribbon(pp_males), xlab = "Time (days)", ylab = "Net survival", label = "men")
+plot1 = plot!(plot1, pp_females.grid, pp_females.Sₑ, ribbon=mkribbon(pp_females), label = "women")
 
-plot2 = plot(pp_young.grid, pp_young.Sₑ, ribbon=(pp_young.Sₑ - lower_bounds_under65, upper_bounds_under65 - pp_young.Sₑ), xlab = "Time (days)", ylab = "Net survival", label = "Under 65")
-plot2 = plot!(pp_old.grid, pp_old.Sₑ, ribbon=(pp_old.Sₑ - lower_bounds_65, upper_bounds_65 - pp_old.Sₑ), xlab = "Time (days)", ylab = "Net survival", label = "65 and up")
+plot2 = plot(pp_young.grid, pp_young.Sₑ, ribbon=mkribbon(pp_young), xlab = "Time (days)", ylab = "Net survival", label = "Under 65")
+plot2 = plot!(pp_old.grid, pp_old.Sₑ, ribbon=mkribbon(pp_old), label = "65 and up")
 
-plot(plot1, plot2, layout = (1, 2))
Example block output
Visually, it is almost immediately understood that there are no worthy differences between the two sexes whereas the age65 variable seems to play a big role.
Estimated sample size and life expectancy
Given that the age group plays a significant role in the study, we will now estimate the sample size by yearly intervals in order to better compare the age groups.
elt, ess = nessie(@formula(Surv(time,status)~age65), colrec, slopop)
-elt
2×2 DataFrame
Rowage65expected_life_time
SymbolFloat64
1young24.7882
2old10.2949
The expected life time for the younger patients is significatively higher than for older patients (24.78 years > 10.29 years).
hcat(ess[:,3]...)
23×2 Matrix{Float64}:
+plot(plot1, plot2, layout = (1, 2))
Example block output
Visually, it is almost immediately understood that there are no worthy differences between the two sexes whereas the age65 variable seems to play a big role.
The same kind of graph can be made on the stage: 
pp3 = fit(PoharPerme, @formula(Surv(time,status)~stage), colrec, slopop)
+plot3 = plot(xlab = "Time (days)", ylab = "Net survival",title="Net survival per cancer stages")
+for i in 1:nrow(pp3)
+    e = pp3[i,:estimator]
+    plot!(plot3, e.grid, e.Sₑ, ribbon=mkribbon(e), label=pp3[i,:stage])
+end
+plot(plot3)
Example block output
Estimated sample size and life expectancy
Given that the age group plays a significant role in the study, we will now estimate the sample size by yearly intervals in order to better compare the age groups.
elt, ess = nessie(@formula(Surv(time,status)~age65), colrec, slopop)
+elt
2×2 DataFrame
Rowage65expected_life_time
SymbolFloat64
1young24.7882
2old10.2949
The expected life time for the younger patients is significatively higher than for older patients (24.78 years > 10.29 years).
hcat(ess[:,3]...)
5×2 Matrix{Float64}:
  2352.0   3619.0
  2323.53  3389.94
  2293.83  3168.85
  2262.95  2954.21
- 2230.76  2746.24
- 2197.39  2546.26
- 2162.79  2352.33
- 2127.0   2162.2
- 2090.07  1974.6
- 2052.34  1792.37
-    ⋮     
- 1846.31  1047.55
- 1801.39   933.738
- 1754.86   827.664
- 1706.36   727.389
- 1655.44   633.358
- 1601.95   546.265
- 1545.51   464.941
- 1484.84   388.018
- 1419.49   316.588
Finally, the table above represents yearly expected sample sizes for both age groups under 65 and above, with the second column representing the latter. We can see that the sample size decreases for the older patients in a much more dramatic way than for the younger ages.
Unsurprisingly, we can thus conclude that age plays an important role in the study.

+ 2230.76  2746.24

Finally, the table above represents yearly expected sample sizes for both age groups under 65 and above, with the second column representing the latter. We can see that the sample size decreases for the older patients in a much more dramatic way than for the younger ages.

Unsurprisingly, we can thus conclude that age plays an important role in the study.

diff --git a/dev/getting_started/index.html b/dev/getting_started/index.html index a0f6ff7..1aa14ff 100644 --- a/dev/getting_started/index.html +++ b/dev/getting_started/index.html @@ -1,2 +1,2 @@ -Getting Started · NetSurvival.jl
GitHub

Getting Started

In many population-based studies, the specific cause of death is unidentified, unreliable or even unavailable. Relative survival analysis addresses this scenario, previously unexplored in general survival analysis. Different methods were created with the aim to construct a consistant and reliable estimator for this purpose.

Net survival settings

Note first that, for any positive random variable $X$, we will use extensively the following functions (that each fully characterize the distribution of $X$):

Function symbol and definitionFunction Name
$S_X(t) = \mathbb P(X > t)$Survival function
$\Lambda_X(t) = -\ln S_X(t)$Cumulative Hazard function
$\lambda_X(t) = \partial \Lambda_X(t)$Instantaneous hazard function

Consider a study that consists of censored survival times from a specific cause. Such a study consists of several random objects:

Random variableNameIs it observed ?
$E$"Excess" lifetime✘
$P$"Population" lifetime✘
$O = E \wedge P$"Overall" lifetime✘
$C$"Censoring" time✘
$\mathbf D$Vector of covariates✔
$T = O \wedge C$Event time✔
$\Delta = \mathbf{1}\{T \leq C\}$Event status✔

It is important to note that we do not observe a a potential indicator $\mathbf{1}\{E \geq P\}$. This is one of the key differences between net survival and standard survival. The standard Net survival analysis solves this problem by assuming that the underlying times $E$ and $P$ are independent from each other.

One central hypothesis in net survival (on top of non-informative censoring) is that $P$ and $E$ are independent. This independence can be written in terms of hazard rates $\lambda_O(t) = \lambda_P(t) + \lambda_E(t)$, or in terms of survival functions $S_O(t) = S_P(t)S_E(t)$.

The population hazard for each individual $\lambda_{P_i}$ is usally drawn from a reference life table, and may depend on covariates $\mathbf D_i$ such as age and date, sex, country, race, etc... See the RateTables.jl package for more details on the potential covariates. On the other hand, the excess mortality is assumed to be i.i.d. between individuals and not to depend on covariates at all. Thus, we mostly omit these covariates from our notations.

Available Estimators

The estimation of net survival is usually discussed in terms of the estimation of the cumulative excess hazard $\Lambda_E(t)$ and/or the instantaneous hazard $\lambda_E = \partial\Lambda_E$. To describe the estimators, we use the following counting processes notations, similar to standard survival analysis(see e.g. [6] or [7]).

  • The uncensored event indicatrix $\partial N_i(t)$ for individual $i$ at time $t$
  • The total number of uncensored events process $\partial N(t) = \sum_i \partial N_i(t)$ at time $t$
  • The at-risk indicatrix $Y_i(t)$, for whether an individual is still at risk
  • The total number at risk process $Y(t) = \sum_i Y_i(t)$ at time $t$

With these definitions and assumptions in mind, we will now present the four different methods implemented in this package, commonly used in literature, to estimate the excess hazard function $\partial\Lambda_E(t)$ and its variance. Recall that we estimate the variance as $\sigma_E^2(t) = \int_{0}^t \partial\sigma_E^2(s)$.

NameReferenceProposed (partial) Excess Hazard $\partial\hat{\Lambda}_E(s)$Proposed (partial) Variance $\partial\hat{\sigma}_E^2(s)$
Ederer I[1]$\frac{\sum_i N_i(s)}{\sum_i Y_i(s)} - \frac{\sum_i S_{P_i}(s)\partial\Lambda_{P_i}(s)}{\sum_i S_{P_i}(s)}$$\frac{\sum_i N_i(s)}{\left(\sum_i Y_i(s)\right)^2}$
Ederer II[2]$\frac{\sum_i N_i(s)}{\sum_i Y_i(s)} - \frac{\sum_i Y_i(s)\partial\Lambda_{P_i}(s)}{\sum_i Y_i(s)}$$\frac{\sum_i N_i(s)}{\left(\sum_i Y_i(s)\right)^2}$
Hakulinen[3]$\frac{\sum_i N_i(s)}{\sum_i Y_i(s)} - \frac{\sum_i \frac{Y_i(s)}{ S_{P_i}(s)}\partial\Lambda_{P_i}(s)}{\sum_i \frac{Y_i(s)}{ S_{P_i}(s)}}$$\frac{\sum_i N_i(s)}{\left(\sum_i Y_i(s)\right)^2}$
Pohar Perme[4]$\frac{\sum_i \frac{\partial N_i(s)}{S_{P_i}(s)} - \sum_i \frac{Y_i(s)}{S_{P_i}(s)}\partial\Lambda_{P_i}(s)}{\sum_i \frac{Y_i(s)}{S_{P_i}(s)}}$$\frac{\sum_{i=1}^n \frac{\partial N_i(s)}{S^2_{P_i}}}{\left(\sum_i \frac{Y_i(s)}{S_{p_i}(s)}\right)^2}$

where, in the variances, it is understood that when no more individuals are at risk $0/0$ gives $0$.

The Pohar Perme estimator [4] is the newest addition to relative survival analysis between the four methods, particularly designed to handle situations where covariates may change over time. It is trusted from the field (see e.g. [8] and [9]) that only this estimator should really be used, the other ones being included mostly for historical reasons and comparisons.

NetSurvival.PoharPerme — Type
PoharPerme

This method estimates net survival probabilities by applying the following estimation:

\[\partial\hat{\Lambda}_E(t) = \frac{\sum_i \frac{dN_i(u)}{S_{P_i}(u)} - \sum_i \frac{Y_i(u)}{S_{P_i}(u)}d\Lambda_{P_i}(u)}{\sum_i \frac{Y_i(u)}{S_{P_i}(u)}}\]

To fit the Pohar Perme to your data based on a certain rate table, apply the example below to your code : 

fit(PoharPerme, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)

References:

  • [4] Perme, Maja Pohar and Stare, Janez and Estève, Jacques (2012). On Estimation in Relative Survival.
source
NetSurvival.EdererI — Type
EdererI

This method estimates net survival probabilities by applying the following estimation:

\[\partial\hat{\Lambda}_E(t) = \frac{\sum_i N_i(u)}{\sum_i Y_i(u)} - \frac{\sum_i S_{P_i}(u)d\Lambda_{P_i}(u)}{\sum_i S_{P_i}(u)}\]

To call this function: 

fit(EdererI, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)
source
NetSurvival.EdererII — Type
EdererII

This method estimates net survival probabilities by applying the following estimation:

\[\partial\hat{\Lambda}_E(t) = \frac{\sum_i N_i(u)}{\sum_i Y_i(u)} - \frac{\sum_i Y_i(u)d\Lambda_{P_i}(u)}{\sum_i Y_i(u)}\]

To call this function: 

fit(EdererII, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)
source
NetSurvival.Hakulinen — Type
Hakulinen

To call this function: 

fit(Hakulinen, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)
source

Crude Mortality

The crude mortality rate is the global mortality rate (from all causes of death) for a population. This measure does not use the cause of death information which, as we previously mentioned, can be unreliable and incomplete. It is given as:

\[M_E(t) = \int_0^t S_O(u-) \lambda_E(u)du\]

There exists a few estimators of this quantity, the most known one being the Cronin-Feuer estimator [10], given by:

\[\hat{M}_E(t) = \int_0^t \hat{S}_O(u-) \hat{\lambda}_E(u)du\]

where, traditionally, Kaplan-Meier is used to estimate the overall survival function $\hat{S}_O$ and Ederer II is used for the excess hazard rate $\hat{\lambda}_E$ and the population hazard rate $\hat{\lambda}_P$.

Our implementation is a bit more permissive, as any net survival estimators can be used for $\hat{\lambda}_E$. Of course, the default is still the original Ederer II which provides the original Cronin-Feuer estimator.

NetSurvival.CrudeMortality — Type
CrudeMortality

This method calculates the crude mortality and presents both the excess mortality and population mortality rates.

The default Cronin-Feuer estimator can be fitted to data with the following interface: 

fit(CrudeMortality, args...)

where the args are passed to fit(EdererII,args...) to compute the excess hazard.

A more direct syntax can be used, specifying directly the estimator for the excess hazard:

CrudeMortality(fit(EdererII,args...))

To use another excess hazard estimator, simply replace EdererII with the method of your choice (PoharPerme, EdererI, Hakulinen).

References:

  • [10] Cronin, Kathleen A and Feuer, Eric J (2000). Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival.
source

Grafféo Log-Rank Test

The Grafféo Log-Rank Test [5] was constructed as a complement to the Pohar Perme estimator, aiming to compare the net survival functions provided by the latter. The test is designed to compare these functions across multiple groups, including stratified covariables, and to ultimately determine, with the given p-value, which covariables are impactful to the study.

The null $(H_0)$ hypothesis tests the following assumption:

\[\forall t \in [0,T], \; \; \Lambda_{E,g_1}(t) = \Lambda_{E,g_2}(t) = ... = \Lambda_{E,g_k}(t),\]

where $G = \{g_1,...,g_k\}$ is a partition of $1,...,n$ consisting of disjoint groups of individuals that we wish to compare to each other. For all group $g \in G$, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by:

  • $\partial N_{E,g}(s) = \sum_{i \in g} \frac{\partial N_i(s)}{S_{P_i}(s)} - \frac{Y_i(s)}{S_{P_i}(s)}\partial\Lambda_{P_i}(s)$
  • $Y_{E,g}(s) = \sum_{i \in g} \frac{Y_i(s)}{S_{P_i}(s)}$
  • $R_{g}(s) = \frac{Y_{E,g}(s)}{\sum_{g\in G} Y_{E,g}(s)}$

Then, define the vector $\mathbf Z = \left(Z_{g_r}: \; r \in 1,...,k-1 \right)$ with entries:

\[Z_g(T) = N_{E,g}(s) - \int_{0}^T Y_{E,g}(s) \partial\hat{\Lambda}_E(s)\]

The test statistic is then given by:

\[U(T) = \mathbf Z(T)'\hat{\Sigma}_Z^{-1} \mathbf Z(T)\]

where the entries of the $\hat{\Sigma}_Z$ matrix are given by:

\[\sigma_{g,h}(T) = \int_0^T \sum_{\ell \in G} \left(\delta_{g,\ell} - R_g(t) \right)\left(\delta_{h,\ell} - R_h(t)\right) \left(\sum_{i\in\ell} \frac{\partial N_i(s)}{S^2_{P_i}}\right)\]

Under $H_0$, the statistic $U(T)$ is asymptotically $\chi^2(k-1)$-distributed. We thus reject the $H_0$ hypothesis when the p-value obtained is under $0.05$, admitting the notable difference between the groups.

NetSurvival.GraffeoTest — Type
GraffeoTest

The Grafféo test is a log-rank type test and is typically used in net survival analysis to determine the impact of certain covariates in the study.

The null $(H_0)$ hypothesis tests the following assumption:

\[\forall t \in [0,T], \Lambda_{E,g_1}(t) = \Lambda_{E,g_2}(t) = ... = \Lambda_{E,g_k}(t)\]

where $G = {g_1,...,g_k}$ is a partition of $1,...,n$ consisting of disjoint groups of individuals that we wish to compare to each other. For all group $g \in G$, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by:

  • $\partial N_{E,g}(s) = \sum_{i \in g} \frac{\partial N_i(s)}{S_{P_i}(s)} - \frac{Y_i(s)}{S_{P_i}(s)}\partial\Lambda_{P_i}(s)$
  • $Y_{E,g}(s) = \sum_{i \in g} \frac{Y_i(s)}{S_{P_i}(s)}$
  • $R_{g}(s) = \frac{Y_{E,g}(s)}{\sum_{g\in G} Y_{E,g}(s)}$

Then, define the vector $\mathbf Z = \left(Z_{g_r}: r \in 1,...,k-1 \right)$ with entries:

\[Z_g(T) = N_{E,g}(s) - \int_{0}^T Y_{E,g}(s) \partial\hat{\Lambda}_E(s)\]

The test statistic is then given by:

\[U(T) = \mathbf Z(T)'\hat{\Sigma}_Z^{-1} \mathbf Z(T)\]

where the entries of the $\hat{\Sigma}_Z$ matrix are given by:

\[\sigma_{g,h}(T) = \int_0^T \sum_{\ell \in G} \left(\delta_{g,\ell} - R_g(t) \right)\left(\delta_{h,\ell} - R_h(t)\right) \left(\sum_{i\in\ell} \frac{\partial N_i(s)}{S^2_{P_i}}\right)\]

Under $H_0$, the statistic $U(T)$ is asymptotically $\chi^2(k-1)$-distributed.

To apply the test to your data based on a certain rate table, apply the example below to your code : 

fit(GraffeoTest, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)

If you wish to stratify a covariate:

fit(GraffeoTest, @formula(Surv(time,status)~covariable1 + Strata(covariable2)), data, ratetable)

The produced test statistics is supposed to follow a chi squared distribution under $(H_0)$. You can fetch the results using the .stat, .df and .pval fields of the returned object.

References:

  • [5] Grafféo, Nathalie and Castell, Fabienne and Belot, Aurélien and Giorgi, Roch (2016). A Log-Rank-Type Test to Compare Net Survival Distributions.
source

Nessie

The Nessie function estimates the sample size by yearly intervals as well as averages an estimated lifespan left for a given group.

This function is highly dependant on the Life function taken from the RateTables.jl package which you can find documented here.

The sample size is thus taken by the following formula:

\[ESS(t) = \sum_i^N S_{P_i}(t)\]

While the estimated lifepsan is directly taken from the expectation function.

NetSurvival.nessie — Function
nessie

To call this function, use the formula below: 

nessie(@formula(Surv(time,status)~covariate), data, ratetable)
source

References

[1]
F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).
[2]
F. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).
[3]
T. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).
[4]
M. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).
[5]
N. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).
[6]
T. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).
[7]
P. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).
[8]
M. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).
[9]
H. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).
[10]
K. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).
+Getting Started · NetSurvival.jl
GitHub

Getting Started

In many population-based studies, the specific cause of death is unidentified, unreliable or even unavailable. Relative survival analysis addresses this scenario, previously unexplored in general survival analysis. Different methods were created with the aim to construct a consistant and reliable estimator for this purpose.

Net survival settings

Note first that, for any positive random variable $X$, we will use extensively the following functions (that each fully characterize the distribution of $X$):

Function symbol and definitionFunction Name
$S_X(t) = \mathbb P(X > t)$Survival function
$\Lambda_X(t) = -\ln S_X(t)$Cumulative Hazard function
$\lambda_X(t) = \partial \Lambda_X(t)$Instantaneous hazard function

Consider a study that consists of censored survival times from a specific cause. Such a study consists of several random objects:

Random variableNameIs it observed ?
$E$"Excess" lifetime✘
$P$"Population" lifetime✘
$O = E \wedge P$"Overall" lifetime✘
$C$"Censoring" time✘
$\mathbf D$Vector of covariates✔
$T = O \wedge C$Event time✔
$\Delta = \mathbf{1}\{T \leq C\}$Event status✔

It is important to note that we do not observe a a potential indicator $\mathbf{1}\{E \geq P\}$. This is one of the key differences between net survival and standard survival. The standard Net survival analysis solves this problem by assuming that the underlying times $E$ and $P$ are independent from each other.

One central hypothesis in net survival (on top of non-informative censoring) is that $P$ and $E$ are independent. This independence can be written in terms of hazard rates $\lambda_O(t) = \lambda_P(t) + \lambda_E(t)$, or in terms of survival functions $S_O(t) = S_P(t)S_E(t)$.

The population hazard for each individual $\lambda_{P_i}$ is usally drawn from a reference life table, and may depend on covariates $\mathbf D_i$ such as age and date, sex, country, race, etc... See the RateTables.jl package for more details on the potential covariates. On the other hand, the excess mortality is assumed to be i.i.d. between individuals and not to depend on covariates at all. Thus, we mostly omit these covariates from our notations.

Available Estimators

The estimation of net survival is usually discussed in terms of the estimation of the cumulative excess hazard $\Lambda_E(t)$ and/or the instantaneous hazard $\lambda_E = \partial\Lambda_E$. To describe the estimators, we use the following counting processes notations, similar to standard survival analysis(see e.g. [6] or [7]).

  • The uncensored event indicatrix $\partial N_i(t)$ for individual $i$ at time $t$
  • The total number of uncensored events process $\partial N(t) = \sum_i \partial N_i(t)$ at time $t$
  • The at-risk indicatrix $Y_i(t)$, for whether an individual is still at risk
  • The total number at risk process $Y(t) = \sum_i Y_i(t)$ at time $t$

With these definitions and assumptions in mind, we will now present the four different methods implemented in this package, commonly used in literature, to estimate the excess hazard function $\partial\Lambda_E(t)$ and its variance. Recall that we estimate the variance as $\sigma_E^2(t) = \int_{0}^t \partial\sigma_E^2(s)$.

NameReferenceProposed (partial) Excess Hazard $\partial\hat{\Lambda}_E(s)$Proposed (partial) Variance $\partial\hat{\sigma}_E^2(s)$
Ederer I[1]$\frac{\sum_i N_i(s)}{\sum_i Y_i(s)} - \frac{\sum_i S_{P_i}(s)\partial\Lambda_{P_i}(s)}{\sum_i S_{P_i}(s)}$$\frac{\sum_i N_i(s)}{\left(\sum_i Y_i(s)\right)^2}$
Ederer II[2]$\frac{\sum_i N_i(s)}{\sum_i Y_i(s)} - \frac{\sum_i Y_i(s)\partial\Lambda_{P_i}(s)}{\sum_i Y_i(s)}$$\frac{\sum_i N_i(s)}{\left(\sum_i Y_i(s)\right)^2}$
Hakulinen[3]$\frac{\sum_i N_i(s)}{\sum_i Y_i(s)} - \frac{\sum_i \frac{Y_i(s)}{ S_{P_i}(s)}\partial\Lambda_{P_i}(s)}{\sum_i \frac{Y_i(s)}{ S_{P_i}(s)}}$$\frac{\sum_i N_i(s)}{\left(\sum_i Y_i(s)\right)^2}$
Pohar Perme[4]$\frac{\sum_i \frac{\partial N_i(s)}{S_{P_i}(s)} - \sum_i \frac{Y_i(s)}{S_{P_i}(s)}\partial\Lambda_{P_i}(s)}{\sum_i \frac{Y_i(s)}{S_{P_i}(s)}}$$\frac{\sum_{i=1}^n \frac{\partial N_i(s)}{S^2_{P_i}}}{\left(\sum_i \frac{Y_i(s)}{S_{p_i}(s)}\right)^2}$

where, in the variances, it is understood that when no more individuals are at risk $0/0$ gives $0$.

The Pohar Perme estimator [4] is the newest addition to relative survival analysis between the four methods, particularly designed to handle situations where covariates may change over time. It is trusted from the field (see e.g. [8] and [9]) that only this estimator should really be used, the other ones being included mostly for historical reasons and comparisons.

NetSurvival.PoharPerme — Type
PoharPerme

This method estimates net survival probabilities by applying the following estimation:

\[\partial\hat{\Lambda}_E(t) = \frac{\sum_i \frac{dN_i(u)}{S_{P_i}(u)} - \sum_i \frac{Y_i(u)}{S_{P_i}(u)}d\Lambda_{P_i}(u)}{\sum_i \frac{Y_i(u)}{S_{P_i}(u)}}\]

To fit the Pohar Perme to your data based on a certain rate table, apply the example below to your code : 

fit(PoharPerme, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)

References:

  • [4] Perme, Maja Pohar and Stare, Janez and Estève, Jacques (2012). On Estimation in Relative Survival.
source
NetSurvival.EdererI — Type
EdererI

This method estimates net survival probabilities by applying the following estimation:

\[\partial\hat{\Lambda}_E(t) = \frac{\sum_i N_i(u)}{\sum_i Y_i(u)} - \frac{\sum_i S_{P_i}(u)d\Lambda_{P_i}(u)}{\sum_i S_{P_i}(u)}\]

To call this function: 

fit(EdererI, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)
source
NetSurvival.EdererII — Type
EdererII

This method estimates net survival probabilities by applying the following estimation:

\[\partial\hat{\Lambda}_E(t) = \frac{\sum_i N_i(u)}{\sum_i Y_i(u)} - \frac{\sum_i Y_i(u)d\Lambda_{P_i}(u)}{\sum_i Y_i(u)}\]

To call this function: 

fit(EdererII, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)
source
NetSurvival.Hakulinen — Type
Hakulinen

To call this function: 

fit(Hakulinen, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)
source

Crude Mortality

The crude mortality rate is the global mortality rate (from all causes of death) for a population. This measure does not use the cause of death information which, as we previously mentioned, can be unreliable and incomplete. It is given as:

\[M_E(t) = \int_0^t S_O(u-) \lambda_E(u)du\]

There exists a few estimators of this quantity, the most known one being the Cronin-Feuer estimator [10], given by:

\[\hat{M}_E(t) = \int_0^t \hat{S}_O(u-) \hat{\lambda}_E(u)du\]

where, traditionally, Kaplan-Meier is used to estimate the overall survival function $\hat{S}_O$ and Ederer II is used for the excess hazard rate $\hat{\lambda}_E$ and the population hazard rate $\hat{\lambda}_P$.

Our implementation is a bit more permissive, as any net survival estimators can be used for $\hat{\lambda}_E$. Of course, the default is still the original Ederer II which provides the original Cronin-Feuer estimator.

NetSurvival.CrudeMortality — Type
CrudeMortality

This method calculates the crude mortality and presents both the excess mortality and population mortality rates.

The default Cronin-Feuer estimator can be fitted to data with the following interface: 

fit(CrudeMortality, args...)

where the args are passed to fit(EdererII,args...) to compute the excess hazard.

A more direct syntax can be used, specifying directly the estimator for the excess hazard:

CrudeMortality(fit(EdererII,args...))

To use another excess hazard estimator, simply replace EdererII with the method of your choice (PoharPerme, EdererI, Hakulinen).

References:

  • [10] Cronin, Kathleen A and Feuer, Eric J (2000). Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival.
source

Grafféo Log-Rank Test

The Grafféo Log-Rank Test [5] was constructed as a complement to the Pohar Perme estimator, aiming to compare the net survival functions provided by the latter. The test is designed to compare these functions across multiple groups, including stratified covariables, and to ultimately determine, with the given p-value, which covariables are impactful to the study.

The null $(H_0)$ hypothesis tests the following assumption:

\[\forall t \in [0,T], \; \; \Lambda_{E,g_1}(t) = \Lambda_{E,g_2}(t) = ... = \Lambda_{E,g_k}(t),\]

where $G = \{g_1,...,g_k\}$ is a partition of $1,...,n$ consisting of disjoint groups of individuals that we wish to compare to each other. For all group $g \in G$, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by:

  • $\partial N_{E,g}(s) = \sum_{i \in g} \frac{\partial N_i(s)}{S_{P_i}(s)} - \frac{Y_i(s)}{S_{P_i}(s)}\partial\Lambda_{P_i}(s)$
  • $Y_{E,g}(s) = \sum_{i \in g} \frac{Y_i(s)}{S_{P_i}(s)}$
  • $R_{g}(s) = \frac{Y_{E,g}(s)}{\sum_{g\in G} Y_{E,g}(s)}$

Then, define the vector $\mathbf Z = \left(Z_{g_r}: \; r \in 1,...,k-1 \right)$ with entries:

\[Z_g(T) = N_{E,g}(s) - \int_{0}^T Y_{E,g}(s) \partial\hat{\Lambda}_E(s)\]

The test statistic is then given by:

\[U(T) = \mathbf Z(T)'\hat{\Sigma}_Z^{-1} \mathbf Z(T)\]

where the entries of the $\hat{\Sigma}_Z$ matrix are given by:

\[\sigma_{g,h}(T) = \int_0^T \sum_{\ell \in G} \left(\delta_{g,\ell} - R_g(t) \right)\left(\delta_{h,\ell} - R_h(t)\right) \left(\sum_{i\in\ell} \frac{\partial N_i(s)}{S^2_{P_i}}\right)\]

Under $H_0$, the statistic $U(T)$ is asymptotically $\chi^2(k-1)$-distributed. We thus reject the $H_0$ hypothesis when the p-value obtained is under $0.05$, admitting the notable difference between the groups.

NetSurvival.GraffeoTest — Type
GraffeoTest

The Grafféo test is a log-rank type test and is typically used in net survival analysis to determine the impact of certain covariates in the study.

The null $(H_0)$ hypothesis tests the following assumption:

\[\forall t \in [0,T], \Lambda_{E,g_1}(t) = \Lambda_{E,g_2}(t) = ... = \Lambda_{E,g_k}(t)\]

where $G = {g_1,...,g_k}$ is a partition of $1,...,n$ consisting of disjoint groups of individuals that we wish to compare to each other. For all group $g \in G$, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by:

  • $\partial N_{E,g}(s) = \sum_{i \in g} \frac{\partial N_i(s)}{S_{P_i}(s)} - \frac{Y_i(s)}{S_{P_i}(s)}\partial\Lambda_{P_i}(s)$
  • $Y_{E,g}(s) = \sum_{i \in g} \frac{Y_i(s)}{S_{P_i}(s)}$
  • $R_{g}(s) = \frac{Y_{E,g}(s)}{\sum_{g\in G} Y_{E,g}(s)}$

Then, define the vector $\mathbf Z = \left(Z_{g_r}: r \in 1,...,k-1 \right)$ with entries:

\[Z_g(T) = N_{E,g}(s) - \int_{0}^T Y_{E,g}(s) \partial\hat{\Lambda}_E(s)\]

The test statistic is then given by:

\[U(T) = \mathbf Z(T)'\hat{\Sigma}_Z^{-1} \mathbf Z(T)\]

where the entries of the $\hat{\Sigma}_Z$ matrix are given by:

\[\sigma_{g,h}(T) = \int_0^T \sum_{\ell \in G} \left(\delta_{g,\ell} - R_g(t) \right)\left(\delta_{h,\ell} - R_h(t)\right) \left(\sum_{i\in\ell} \frac{\partial N_i(s)}{S^2_{P_i}}\right)\]

Under $H_0$, the statistic $U(T)$ is asymptotically $\chi^2(k-1)$-distributed.

To apply the test to your data based on a certain rate table, apply the example below to your code : 

fit(GraffeoTest, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)

If you wish to stratify a covariate:

fit(GraffeoTest, @formula(Surv(time,status)~covariable1 + Strata(covariable2)), data, ratetable)

The produced test statistics is supposed to follow a chi squared distribution under $(H_0)$. You can fetch the results using the .stat, .df and .pval fields of the returned object.

References:

  • [5] Grafféo, Nathalie and Castell, Fabienne and Belot, Aurélien and Giorgi, Roch (2016). A Log-Rank-Type Test to Compare Net Survival Distributions.
source

Nessie

The Nessie function estimates the sample size by yearly intervals as well as averages an estimated lifespan left for a given group.

This function is highly dependant on the Life function taken from the RateTables.jl package which you can find documented here.

The sample size is thus taken by the following formula:

\[ESS(t) = \sum_i^N S_{P_i}(t)\]

While the estimated lifepsan is directly taken from the expectation function.

NetSurvival.nessie — Function
nessie

To call this function, use the formula below: 

nessie(@formula(Surv(time,status)~covariate), data, ratetable)
source

References

[1]
F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).
[2]
F. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).
[3]
T. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).
[4]
M. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).
[5]
N. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).
[6]
T. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).
[7]
P. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).
[8]
M. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).
[9]
H. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).
[10]
K. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).
diff --git a/dev/index.html b/dev/index.html index 06b1072..8ca3b8c 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ -Home · NetSurvival.jl
GitHub

Introduction

The NetSurvival.jl package provides the necessary tools to perform estimations and analysis in the Net Survival field. This specialized branch of Survival Analysis focuses on estimating the probability of survival from a specific event of interest, for example a given cancer, without considering other causes of death. This is especially relevant in the (unfortunately quite common) case where the cause of death indicatrix is either unavailable or untrustworthy. Consequently, the so-called missing indicatrix issue forbids the use of standard competitive risks survival analysis methods on these datasets. For that, a few standard estimators were established in the last 50 years, backed by a wide literature.

By integrating observed data from the target population with historical population mortality data (usually sourced from national census datasets), Net Survival allows the extraction of the specific mortality hazard associated with the particular disease, even under the missing indicatrix issue. The concept of relative survival analysis dates back several decades to the seminal article by Ederer, Axtell, and Cutler in 1961 [1] and the one by Ederer and Heise in 1959 [2].

For years, the Hakulinen estimator (1977) [3] and the Ederer I and II estimators were widely regarded as the gold standard for non-parametric survival curve estimation. However, the introduction of the Pohar-Perme, Stare, and Estève estimator in 2012 [4] resolved several issues inherent in previous estimators, providing a reliable and consistent non-parametric estimator for net survival analysis.

Features

Standard tools nowadays are composed of R packages, with underlying C and C++ routines, that are hard to read, maintain, and use. This package is an attempt to bring standard relative survival analysis modeling routines to Julia, while providing an interface that is close to the relsurv standard, albeit significantly faster and easier to maintain in the future. Our hope is that the junction with classical modeling API in Julia will allow later extensions of the existing modeling methods, with a simple interface for the practitioners.

Some key features in NetSurvival.jl are:

  • A panel of different non-parametric net survival estimators (Ederer I [1], Ederer II [2], Hakulinen [3], Pohar Perme [4]) with an interface compliant with Julia's standards.
  • Grafféo's log-rank test [5] to compare net survival curves accross groups, including stratified testing.
  • Crude mortality, Expected Sample Size, and other useful metrics in net survival field.
  • A 'Nessie' function that outputs the estimated sample size by yearly intervals and the average lifespan expectancy left for a given group.
  • A compact, readable and efficient codebase (up to 1000x less LOC than relsurv for the same functionalities), ensuring long-term maintenability.
  • Significant performance improvements (up to 50x) compared to the R package relsurv.

Installation

The package is not yet available on Julia's general registry, and thus can be installed through the following command:

using Pkg
-Pkg.add("https://github.com/JuliaSurv/NetSurvival.jl.git")

See the rest of this documentation to have a glimpse of the functionalities!

References

[1]
F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).
[2]
F. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).
[3]
T. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).
[4]
M. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).
[5]
N. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).
+Home · NetSurvival.jl
GitHub

Introduction

The NetSurvival.jl package provides the necessary tools to perform estimations and analysis in the Net Survival field. This specialized branch of Survival Analysis focuses on estimating the probability of survival from a specific event of interest, for example a given cancer, without considering other causes of death. This is especially relevant in the (unfortunately quite common) case where the cause of death indicatrix is either unavailable or untrustworthy. Consequently, the so-called missing indicatrix issue forbids the use of standard competitive risks survival analysis methods on these datasets. For that, a few standard estimators were established in the last 50 years, backed by a wide literature.

By integrating observed data from the target population with historical population mortality data (usually sourced from national census datasets), Net Survival allows the extraction of the specific mortality hazard associated with the particular disease, even under the missing indicatrix issue. The concept of relative survival analysis dates back several decades to the seminal article by Ederer, Axtell, and Cutler in 1961 [1] and the one by Ederer and Heise in 1959 [2].

For years, the Hakulinen estimator (1977) [3] and the Ederer I and II estimators were widely regarded as the gold standard for non-parametric survival curve estimation. However, the introduction of the Pohar Perme, Stare, and Estève estimator in 2012 [4] resolved several issues inherent in previous estimators, providing a reliable and consistent non-parametric estimator for net survival analysis.

Features

Standard tools nowadays are composed of R packages, with underlying C and C++ routines, that are hard to read, maintain, and use. This package is an attempt to bring standard relative survival analysis modeling routines to Julia, while providing an interface that is close to the relsurv standard, albeit significantly faster and easier to maintain in the future. Our hope is that the junction with classical modeling API in Julia will allow later extensions of the existing modeling methods, with a simple interface for the practitioners.

Some key features in NetSurvival.jl are:

  • A panel of different non-parametric net survival estimators (Ederer I [1], Ederer II [2], Hakulinen [3], Pohar Perme [4]) with an interface compliant with Julia's standards.
  • Grafféo's log-rank test [5] to compare net survival curves accross groups, including stratified testing.
  • Crude mortality, Expected Sample Size, and other useful metrics in net survival field.
  • A 'Nessie' function that outputs the estimated sample size by yearly intervals and the average lifespan expectancy left for a given group.
  • A compact, readable and efficient codebase (up to 1000x less LOC than relsurv for the same functionalities), ensuring long-term maintenability.
  • Significant performance improvements (up to 50x) compared to the R package relsurv.

Installation

The package is not yet available on Julia's general registry, and thus can be installed through the following command:

using Pkg
+Pkg.add("https://github.com/JuliaSurv/NetSurvival.jl.git")

See the rest of this documentation to have a glimpse of the functionalities!

References

[1]
F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).
[2]
F. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).
[3]
T. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).
[4]
M. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).
[5]
N. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).
diff --git a/dev/references/index.html b/dev/references/index.html index 6217dd0..52482d9 100644 --- a/dev/references/index.html +++ b/dev/references/index.html @@ -1,2 +1,2 @@ -References · NetSurvival.jl
GitHub

References

[1]
F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).
[2]
F. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).
[3]
T. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).
[4]
M. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).
[5]
N. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).
[6]
T. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).
[7]
P. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).
[8]
M. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).
[9]
H. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).
[10]
K. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).
[11]
T. Hakulinen and K. H. Abeywickrama. A computer program package for relative survival analysis. Computer programs in biomedicine 19, 197–207 (1985).
[12]
T. Hakulinen and L. Tenkanen. Regression analysis of relative survival rates. Journal of the Royal Statistical Society Series C: Applied Statistics 36, 309–317 (1987).
+References · NetSurvival.jl
GitHub

References

[1]
F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).
[2]
F. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).
[3]
T. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).
[4]
M. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).
[5]
N. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).
[6]
T. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).
[7]
P. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).
[8]
M. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).
[9]
H. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).
[10]
K. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).
[11]
T. Hakulinen and K. H. Abeywickrama. A computer program package for relative survival analysis. Computer programs in biomedicine 19, 197–207 (1985).
[12]
T. Hakulinen and L. Tenkanen. Regression analysis of relative survival rates. Journal of the Royal Statistical Society Series C: Applied Statistics 36, 309–317 (1987).
diff --git a/dev/search_index.js b/dev/search_index.js index 3b4b21f..d346948 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"references/#References","page":"References","title":"References","text":"","category":"section"},{"location":"references/","page":"References","title":"References","text":"F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).\n\n\n\nF. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).\n\n\n\nT. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).\n\n\n\nM. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).\n\n\n\nN. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).\n\n\n\nT. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).\n\n\n\nP. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).\n\n\n\nM. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).\n\n\n\nH. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).\n\n\n\nK. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).\n\n\n\n","category":"page"},{"location":"references/","page":"References","title":"References","text":"T. Hakulinen and K. H. Abeywickrama. A computer program package for relative survival analysis. Computer programs in biomedicine 19, 197–207 (1985).\n\n\n\nT. Hakulinen and L. Tenkanen. Regression analysis of relative survival rates. Journal of the Royal Statistical Society Series C: Applied Statistics 36, 309–317 (1987).\n\n\n\n","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"CurrentModule = NetSurvival","category":"page"},{"location":"benches/#Benchmarking-results","page":"Benchmarking results","title":"Benchmarking results","text":"","category":"section"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"This page provides benchmark results of several standard net survival routines implemented in this package. Note that the runtime also depends on other packages, in particular on RateTables.jl. ","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"note: Take numbers displayed here carefully\nAll the following benchmarks are run on github action contunous integration platform, which is a very slow computing engine. Thus, these numbers may not represent your local performance. A locally ran version of these benchmarks in availiable on the github readme, and the below code blocks can be used to check performance on your own hardware.","category":"page"},{"location":"benches/#Benchmarks-w.r.t.-relsurv","page":"Benchmarking results","title":"Benchmarks w.r.t. relsurv","text":"","category":"section"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"This first set of benchmark compares standard functionalities with their implementation in relsurv. Below numbers gives runtime mulitpliers w.r.t. R::relsurv, computed on github action CI.","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"using RateTables, NetSurvival, RCall, DataFrames\n\nfunction test_surv(r_method,::Type{E}, stratified) where E\n if stratified\n jl = @timed fit(E, @formula(Surv(time,status)~sex), colrec, slopop)\n @rput r_method\n r = @timed R\"\"\"\n rez = relsurv::rs.surv(survival::Surv(time, stat) ~ sex, rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop, method = r_method, add.times=1:8149)\n \"\"\"\n else\n jl = @timed fit(E, @formula(Surv(time,status)~1), colrec, slopop)\n @rput r_method\n r = @timed R\"\"\"\n rez = relsurv::rs.surv(survival::Surv(time, stat) ~ 1, rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop, method = r_method, add.times=1:8149)\n \"\"\"\n end\n return r.time / jl.time\nend\n\nfunction test_graffeo(stratified)\n if stratified\n jl = @timed fit(GraffeoTest, @formula(Surv(time,status)~stage+Strata(sex)), colrec, slopop)\n r = @timed R\"\"\"\n rez = relsurv::rs.diff(survival::Surv(time, stat) ~ stage + survival::strata(sex), rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop)\n \"\"\"\n else \n jl = @timed fit(GraffeoTest, @formula(Surv(time,status)~stage), colrec, slopop)\n r = @timed R\"\"\"\n rez = relsurv::rs.diff(survival::Surv(time, stat) ~ stage, rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop)\n \"\"\"\n end\n return r.time / jl.time \nend\ntest_all(stratified) = [\n test_surv(\"pohar-perme\", PoharPerme, stratified),\n test_surv(\"ederer1\", EdererI, stratified),\n test_surv(\"ederer2\", EdererII, stratified),\n test_surv(\"hakulinen\", Hakulinen, stratified),\n test_graffeo(stratified),\n]\ntest_all() = DataFrame(\n Algorithm = [\"Pohar Perme\", \"EdererI\", \"EdererII\", \"Hakulinen\", \"Graffeo's LRT\"], \n unstratified = test_all(false), \n stratified = test_all(true)\n)\nrez = test_all()\nrez = test_all() # discard first run.\n\n# note: to obtain the pretty printing from the readme, you need to install PrettyTables.jl and do : \n# using PrettyTables\n# pretty_table(rez, backend = Val(:markdown))\n\nrez","category":"page"},{"location":"benches/#Benchmarking-across-time","page":"Benchmarking results","title":"Benchmarking across time","text":"","category":"section"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"The following charts provide a glimpse of NetSurvival.jl's performance along time, also ran on github CI: ","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"CurrentModule = NetSurvival","category":"page"},{"location":"getting_started/#Getting-Started","page":"Getting Started","title":"Getting Started","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"In many population-based studies, the specific cause of death is unidentified, unreliable or even unavailable. Relative survival analysis addresses this scenario, previously unexplored in general survival analysis. Different methods were created with the aim to construct a consistant and reliable estimator for this purpose.","category":"page"},{"location":"getting_started/#Net-survival-settings","page":"Getting Started","title":"Net survival settings","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Note first that, for any positive random variable X, we will use extensively the following functions (that each fully characterize the distribution of X): ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Function symbol and definition Function Name\nS_X(t) = mathbb P(X t) Survival function\nLambda_X(t) = -ln S_X(t) Cumulative Hazard function\nlambda_X(t) = partial Lambda_X(t) Instantaneous hazard function","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Consider a study that consists of censored survival times from a specific cause. Such a study consists of several random objects: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Random variable Name Is it observed ?\nE \"Excess\" lifetime ✘\nP \"Population\" lifetime ✘\nO = E wedge P \"Overall\" lifetime ✘\nC \"Censoring\" time ✘\nmathbf D Vector of covariates ✔\nT = O wedge C Event time ✔\nDelta = mathbf1T leq C Event status ✔","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"It is important to note that we do not observe a a potential indicator mathbf1E geq P. This is one of the key differences between net survival and standard survival. The standard Net survival analysis solves this problem by assuming that the underlying times E and P are independent from each other.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"One central hypothesis in net survival (on top of non-informative censoring) is that P and E are independent. This independence can be written in terms of hazard rates lambda_O(t) = lambda_P(t) + lambda_E(t), or in terms of survival functions S_O(t) = S_P(t)S_E(t).","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The population hazard for each individual lambda_P_i is usally drawn from a reference life table, and may depend on covariates mathbf D_i such as age and date, sex, country, race, etc... See the RateTables.jl package for more details on the potential covariates. On the other hand, the excess mortality is assumed to be i.i.d. between individuals and not to depend on covariates at all. Thus, we mostly omit these covariates from our notations.","category":"page"},{"location":"getting_started/#Available-Estimators","page":"Getting Started","title":"Available Estimators","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The estimation of net survival is usually discussed in terms of the estimation of the cumulative excess hazard Lambda_E(t) and/or the instantaneous hazard lambda_E = partialLambda_E. To describe the estimators, we use the following counting processes notations, similar to standard survival analysis(see e.g. [6] or [7]). ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The uncensored event indicatrix partial N_i(t) for individual i at time t \nThe total number of uncensored events process partial N(t) = sum_i partial N_i(t) at time t\nThe at-risk indicatrix Y_i(t), for whether an individual is still at risk \nThe total number at risk process Y(t) = sum_i Y_i(t) at time t","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"With these definitions and assumptions in mind, we will now present the four different methods implemented in this package, commonly used in literature, to estimate the excess hazard function partialLambda_E(t) and its variance. Recall that we estimate the variance as sigma_E^2(t) = int_0^t partialsigma_E^2(s). ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Name Reference Proposed (partial) Excess Hazard partialhatLambda_E(s) Proposed (partial) Variance partialhatsigma_E^2(s)\nEderer I [1] fracsum_i N_i(s)sum_i Y_i(s) - fracsum_i S_P_i(s)partialLambda_P_i(s)sum_i S_P_i(s) fracsum_i N_i(s)left(sum_i Y_i(s)right)^2\nEderer II [2] fracsum_i N_i(s)sum_i Y_i(s) - fracsum_i Y_i(s)partialLambda_P_i(s)sum_i Y_i(s) fracsum_i N_i(s)left(sum_i Y_i(s)right)^2\nHakulinen [3] fracsum_i N_i(s)sum_i Y_i(s) - fracsum_i fracY_i(s) S_P_i(s)partialLambda_P_i(s)sum_i fracY_i(s) S_P_i(s) fracsum_i N_i(s)left(sum_i Y_i(s)right)^2\nPohar Perme [4] fracsum_i fracpartial N_i(s)S_P_i(s) - sum_i fracY_i(s)S_P_i(s)partialLambda_P_i(s)sum_i fracY_i(s)S_P_i(s) fracsum_i=1^n fracpartial N_i(s)S^2_P_ileft(sum_i fracY_i(s)S_p_i(s)right)^2","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where, in the variances, it is understood that when no more individuals are at risk 00 gives 0. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The Pohar Perme estimator [4] is the newest addition to relative survival analysis between the four methods, particularly designed to handle situations where covariates may change over time. It is trusted from the field (see e.g. [8] and [9]) that only this estimator should really be used, the other ones being included mostly for historical reasons and comparisons. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"PoharPerme\nEdererI\nEdererII\nHakulinen","category":"page"},{"location":"getting_started/#NetSurvival.PoharPerme","page":"Getting Started","title":"NetSurvival.PoharPerme","text":"PoharPerme\n\nThis method estimates net survival probabilities by applying the following estimation:\n\npartialhatLambda_E(t) = fracsum_i fracdN_i(u)S_P_i(u) - sum_i fracY_i(u)S_P_i(u)dLambda_P_i(u)sum_i fracY_i(u)S_P_i(u)\n\nTo fit the Pohar Perme to your data based on a certain rate table, apply the example below to your code : \n\nfit(PoharPerme, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\nReferences: \n\n[4] Perme, Maja Pohar and Stare, Janez and Estève, Jacques (2012). On Estimation in Relative Survival.\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#NetSurvival.EdererI","page":"Getting Started","title":"NetSurvival.EdererI","text":"EdererI\n\nThis method estimates net survival probabilities by applying the following estimation:\n\npartialhatLambda_E(t) = fracsum_i N_i(u)sum_i Y_i(u) - fracsum_i S_P_i(u)dLambda_P_i(u)sum_i S_P_i(u)\n\nTo call this function: \n\nfit(EdererI, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#NetSurvival.EdererII","page":"Getting Started","title":"NetSurvival.EdererII","text":"EdererII\n\nThis method estimates net survival probabilities by applying the following estimation:\n\npartialhatLambda_E(t) = fracsum_i N_i(u)sum_i Y_i(u) - fracsum_i Y_i(u)dLambda_P_i(u)sum_i Y_i(u)\n\nTo call this function: \n\nfit(EdererII, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#NetSurvival.Hakulinen","page":"Getting Started","title":"NetSurvival.Hakulinen","text":"Hakulinen\n\nTo call this function: \n\nfit(Hakulinen, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#Crude-Mortality","page":"Getting Started","title":"Crude Mortality","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The crude mortality rate is the global mortality rate (from all causes of death) for a population. This measure does not use the cause of death information which, as we previously mentioned, can be unreliable and incomplete. It is given as:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"M_E(t) = int_0^t S_O(u-) lambda_E(u)du","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"There exists a few estimators of this quantity, the most known one being the Cronin-Feuer estimator [10], given by:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"hatM_E(t) = int_0^t hatS_O(u-) hatlambda_E(u)du","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where, traditionally, Kaplan-Meier is used to estimate the overall survival function hatS_O and Ederer II is used for the excess hazard rate hatlambda_E and the population hazard rate hatlambda_P.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Our implementation is a bit more permissive, as any net survival estimators can be used for hatlambda_E. Of course, the default is still the original Ederer II which provides the original Cronin-Feuer estimator. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"CrudeMortality","category":"page"},{"location":"getting_started/#NetSurvival.CrudeMortality","page":"Getting Started","title":"NetSurvival.CrudeMortality","text":"CrudeMortality\n\nThis method calculates the crude mortality and presents both the excess mortality and population mortality rates.\n\nThe default Cronin-Feuer estimator can be fitted to data with the following interface: \n\nfit(CrudeMortality, args...)\n\nwhere the args are passed to fit(EdererII,args...) to compute the excess hazard.\n\nA more direct syntax can be used, specifying directly the estimator for the excess hazard:\n\nCrudeMortality(fit(EdererII,args...))\n\nTo use another excess hazard estimator, simply replace EdererII with the method of your choice (PoharPerme, EdererI, Hakulinen).\n\nReferences: \n\n[10] Cronin, Kathleen A and Feuer, Eric J (2000). Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival.\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#Grafféo-Log-Rank-Test","page":"Getting Started","title":"Grafféo Log-Rank Test","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The Grafféo Log-Rank Test [5] was constructed as a complement to the Pohar Perme estimator, aiming to compare the net survival functions provided by the latter. The test is designed to compare these functions across multiple groups, including stratified covariables, and to ultimately determine, with the given p-value, which covariables are impactful to the study.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The null (H_0) hypothesis tests the following assumption:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"forall t in 0T Lambda_Eg_1(t) = Lambda_Eg_2(t) = = Lambda_Eg_k(t)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where G = g_1g_k is a partition of 1n consisting of disjoint groups of individuals that we wish to compare to each other. For all group g in G, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"partial N_Eg(s) = sum_i in g fracpartial N_i(s)S_P_i(s) - fracY_i(s)S_P_i(s)partialLambda_P_i(s)\nY_Eg(s) = sum_i in g fracY_i(s)S_P_i(s)\nR_g(s) = fracY_Eg(s)sum_gin G Y_Eg(s)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Then, define the vector mathbf Z = left(Z_g_r r in 1k-1 right) with entries: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Z_g(T) = N_Eg(s) - int_0^T Y_Eg(s) partialhatLambda_E(s)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The test statistic is then given by:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"U(T) = mathbf Z(T)hatSigma_Z^-1 mathbf Z(T)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where the entries of the hatSigma_Z matrix are given by: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"sigma_gh(T) = int_0^T sum_ell in G left(delta_gell - R_g(t) right)left(delta_hell - R_h(t)right) left(sum_iinell fracpartial N_i(s)S^2_P_iright)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Under H_0, the statistic U(T) is asymptotically chi^2(k-1)-distributed. We thus reject the H_0 hypothesis when the p-value obtained is under 005, admitting the notable difference between the groups. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"GraffeoTest","category":"page"},{"location":"getting_started/#NetSurvival.GraffeoTest","page":"Getting Started","title":"NetSurvival.GraffeoTest","text":"GraffeoTest\n\nThe Grafféo test is a log-rank type test and is typically used in net survival analysis to determine the impact of certain covariates in the study.\n\nThe null (H_0) hypothesis tests the following assumption:\n\nforall t in 0T Lambda_Eg_1(t) = Lambda_Eg_2(t) = = Lambda_Eg_k(t)\n\nwhere G = g_1g_k is a partition of 1n consisting of disjoint groups of individuals that we wish to compare to each other. For all group g in G, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by: \n\npartial N_Eg(s) = sum_i in g fracpartial N_i(s)S_P_i(s) - fracY_i(s)S_P_i(s)partialLambda_P_i(s)\nY_Eg(s) = sum_i in g fracY_i(s)S_P_i(s)\nR_g(s) = fracY_Eg(s)sum_gin G Y_Eg(s)\n\nThen, define the vector mathbf Z = left(Z_g_r r in 1k-1 right) with entries: \n\nZ_g(T) = N_Eg(s) - int_0^T Y_Eg(s) partialhatLambda_E(s)\n\nThe test statistic is then given by:\n\nU(T) = mathbf Z(T)hatSigma_Z^-1 mathbf Z(T)\n\nwhere the entries of the hatSigma_Z matrix are given by: \n\nsigma_gh(T) = int_0^T sum_ell in G left(delta_gell - R_g(t) right)left(delta_hell - R_h(t)right) left(sum_iinell fracpartial N_i(s)S^2_P_iright)\n\nUnder H_0, the statistic U(T) is asymptotically chi^2(k-1)-distributed.\n\nTo apply the test to your data based on a certain rate table, apply the example below to your code : \n\nfit(GraffeoTest, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\nIf you wish to stratify a covariate:\n\nfit(GraffeoTest, @formula(Surv(time,status)~covariable1 + Strata(covariable2)), data, ratetable)\n\nThe produced test statistics is supposed to follow a chi squared distribution under (H_0). You can fetch the results using the .stat, .df and .pval fields of the returned object. \n\nReferences: \n\n[5] Grafféo, Nathalie and Castell, Fabienne and Belot, Aurélien and Giorgi, Roch (2016). A Log-Rank-Type Test to Compare Net Survival Distributions. \n\n\n\n\n\n","category":"type"},{"location":"getting_started/#Nessie","page":"Getting Started","title":"Nessie","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The Nessie function estimates the sample size by yearly intervals as well as averages an estimated lifespan left for a given group. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"This function is highly dependant on the Life function taken from the RateTables.jl package which you can find documented here.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The sample size is thus taken by the following formula:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"ESS(t) = sum_i^N S_P_i(t)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"While the estimated lifepsan is directly taken from the expectation function. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"nessie","category":"page"},{"location":"getting_started/#NetSurvival.nessie","page":"Getting Started","title":"NetSurvival.nessie","text":"nessie\n\nTo call this function, use the formula below: \n\nnessie(@formula(Surv(time,status)~covariate), data, ratetable)\n\n\n\n\n\n","category":"function"},{"location":"getting_started/#References","page":"Getting Started","title":"References","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).\n\n\n\nF. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).\n\n\n\nT. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).\n\n\n\nM. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).\n\n\n\nN. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).\n\n\n\nT. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).\n\n\n\nP. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).\n\n\n\nM. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).\n\n\n\nH. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).\n\n\n\nK. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).\n\n\n\n","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = NetSurvival","category":"page"},{"location":"#Introduction","page":"Home","title":"Introduction","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"The NetSurvival.jl package provides the necessary tools to perform estimations and analysis in the Net Survival field. This specialized branch of Survival Analysis focuses on estimating the probability of survival from a specific event of interest, for example a given cancer, without considering other causes of death. This is especially relevant in the (unfortunately quite common) case where the cause of death indicatrix is either unavailable or untrustworthy. Consequently, the so-called missing indicatrix issue forbids the use of standard competitive risks survival analysis methods on these datasets. For that, a few standard estimators were established in the last 50 years, backed by a wide literature.","category":"page"},{"location":"","page":"Home","title":"Home","text":"By integrating observed data from the target population with historical population mortality data (usually sourced from national census datasets), Net Survival allows the extraction of the specific mortality hazard associated with the particular disease, even under the missing indicatrix issue. The concept of relative survival analysis dates back several decades to the seminal article by Ederer, Axtell, and Cutler in 1961 [1] and the one by Ederer and Heise in 1959 [2].","category":"page"},{"location":"","page":"Home","title":"Home","text":"For years, the Hakulinen estimator (1977) [3] and the Ederer I and II estimators were widely regarded as the gold standard for non-parametric survival curve estimation. However, the introduction of the Pohar-Perme, Stare, and Estève estimator in 2012 [4] resolved several issues inherent in previous estimators, providing a reliable and consistent non-parametric estimator for net survival analysis.","category":"page"},{"location":"#Features","page":"Home","title":"Features","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Standard tools nowadays are composed of R packages, with underlying C and C++ routines, that are hard to read, maintain, and use. This package is an attempt to bring standard relative survival analysis modeling routines to Julia, while providing an interface that is close to the relsurv standard, albeit significantly faster and easier to maintain in the future. Our hope is that the junction with classical modeling API in Julia will allow later extensions of the existing modeling methods, with a simple interface for the practitioners.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Some key features in NetSurvival.jl are:","category":"page"},{"location":"","page":"Home","title":"Home","text":"A panel of different non-parametric net survival estimators (Ederer I [1], Ederer II [2], Hakulinen [3], Pohar Perme [4]) with an interface compliant with Julia's standards. \nGrafféo's log-rank test [5] to compare net survival curves accross groups, including stratified testing.\nCrude mortality, Expected Sample Size, and other useful metrics in net survival field.\nA 'Nessie' function that outputs the estimated sample size by yearly intervals and the average lifespan expectancy left for a given group. \nA compact, readable and efficient codebase (up to 1000x less LOC than relsurv for the same functionalities), ensuring long-term maintenability.\nSignificant performance improvements (up to 50x) compared to the R package relsurv.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"The package is not yet available on Julia's general registry, and thus can be installed through the following command:","category":"page"},{"location":"","page":"Home","title":"Home","text":"using Pkg\nPkg.add(\"https://github.com/JuliaSurv/NetSurvival.jl.git\")","category":"page"},{"location":"","page":"Home","title":"Home","text":"See the rest of this documentation to have a glimpse of the functionalities!","category":"page"},{"location":"#References","page":"Home","title":"References","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).\n\n\n\nF. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).\n\n\n\nT. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).\n\n\n\nM. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).\n\n\n\nN. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).\n\n\n\n","category":"page"},{"location":"example/#Case-study","page":"Case study","title":"Case study","text":"","category":"section"},{"location":"example/#Introduction-and-datasets","page":"Case study","title":"Introduction and datasets","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"We will illustrate with an example using the dataset colrec, which comprises 5971 patients diagnosed with colon or rectal cancer between 1994 and 2000. This dataset is sourced from the Slovenia cancer registry. Given the high probability that the patients are Slovenian, we will be using the Slovenian mortality table slopop as reference for the populational rates. Subsequently, we can apply various non-parametric estimators for net survival analysis.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"note: N.B.\nMortality tables may vary in structure, with options such as the addition or removal of specific covariates. To confirm that the mortality table is in the correct format, please refer to RateTables.jl's documentation, or directly extract it from there.","category":"page"},{"location":"example/#Cohort-details","page":"Case study","title":"Cohort details","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"The patients in the study are diagnosed between January 1st 1994 and December 31st 2000. Before we move on to the survival probabilities, it is important to be aware of how your data is distributed and of what it comprises. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"using NetSurvival, RateTables, DataFrames\n\nfirst(colrec,10)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Let's explore how our data is distributed first, starting with the age variable.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"println(minimum(colrec.age./365.241), \" ; \",maximum(colrec.age./365.241))","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The study can be considered diverse in terms of age seeing as the patients are between 12 and 96 years old, approximately. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"using Plots \n\nplot(\n histogram(colrec.age./365.241, label=\"Age\"),\n histogram(colrec.time./365.241, label=\"Follow-up time\")\n)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The graph above show us that although it has a wide range of patients within all age groups, it is mostly centered around older adults and elderly, with the majority of the patients being between 60 and 80 years old. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Looking at the second graph that details the distribution of the follow-up times, we notice that the values quickly drop. Unfortunately, this is a common theme in cancer studies. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Let's take a look at the sex variable now: ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"combine(groupby(colrec, :sex), nrow)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"This dataframe shows the number of male and female patients. There isn't too big of a difference between the two. We can say this study includes both gender relatively equally, thus, reducing bias. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"With these two observations, it is also worth noting that colorectal cancer is most common with men and people older than 50.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"In total, we note that we have 5971 patients. By taking a look at the status variable, we can determine the deaths and censorship:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"sum(colrec.status)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Out of the 5971 patients, 4979 have died. This is a very high mortality rate, and again, unfortunately common in cancer studies.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"(nrow(colrec) - sum(colrec.status)) / nrow(colrec)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"In other terms, the censorship rate is of 16.6%, meaning the event, in this case death, was not observed for only 16.6% of the patients. This is a low censorship rate and thus the quality of the signal will be pretty good. ","category":"page"},{"location":"example/#Mortality-table","page":"Case study","title":"Mortality table","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"We will be using the mortality table slopop taken from the RateTables.jl package as the study is done on Slovenian patients. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"slopop","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"By examining slopop, we notice it contains information regarding age and year, as expected for mortality tables. Additionally, it incorporates the covariate sex, which has two possible entries (:male or :female). The ratetable is then three dimensional, with the covariate sex added. For example, the daily hazard rate for a woman turning 45 on the January 1st 2006 can be accessed through the following command:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"daily_hazard(slopop, 45*365.241, 2006*365.241; sex=:female)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Making the survival probability easily calculated with:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"exp(-(daily_hazard(slopop, 45*365.241, 2006*365.241; sex=:female))*365)","category":"page"},{"location":"example/#Overall-and-expected-survival","page":"Case study","title":"Overall and expected survival","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"For this part, we will be using the Survival.jl package to apply the Kaplan Meier estimator for the overall survival.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"using Survival \n\nkm = fit(KaplanMeier, colrec.time./365.241, colrec.status)\nplot(km.events.time, km.survival, label=false, title = \"Kaplan-Meier Estimator for the Overall Survival\")","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The graph above indicates a significant dip in survival probability within the first 5 years, and even more at 10 years. This period in the study is then crucial for the analysis. ","category":"page"},{"location":"example/#Estimated-net-survival","page":"Case study","title":"Estimated net survival","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"In this part, we are interested in the first 5 years of the study. We will thus limit the follow-up time to 5 years, meaning we will censor all individuals with a follow-up time that is higher than this. Then, we will apply the different net survival methods.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"colrec.time5 .= 0.0\ncolrec.status5 .= Bool(true)\nfor i in 1:nrow(colrec)\n colrec.time5[i] = min(colrec.time[i], round(5*365.241))\n if colrec.time[i] > 5*365.241\n colrec.status5[i] = false\n end\nend","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Now that we have defined our own time and status variables according to the observations made, we can apply the different non parametric methods for relative survival.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"e1 = fit(EdererI, @formula(Surv(time5,status5)~1), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"With the EdererI method, after 1826 days have passed, we can say that the survival rate at this mark is around 0456, in the hypothetical world where patients can only die of cancer.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"crude_e1 = CrudeMortality(e1)\nprintln(crude_e1.Mₒ[1826], \" , \", crude_e1.Mₑ[1826], \" , \", crude_e1.Mₚ[1826])","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Out of the 0.63 patients that have died, according to the EdererI method, 0.51 died because of colorectal cancer and 0.12 died of other causes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"e2 = fit(EdererII, @formula(Surv(time5,status5)~1), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Similarily, the EdererII method, also known as the conditional method, shows that at the 5 year mark, the survival probability is of 044 in this hypothetical world.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"crude_e2 = CrudeMortality(e2)\nprintln(crude_e2.Mₒ[1826], \" , \", crude_e2.Mₑ[1826], \" , \", crude_e2.Mₚ[1826])","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Here, out of the 0.63 patients that have dued, 0.53 are due to colorectal cancer and 0.1 due to other causes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"pp = fit(PoharPerme, @formula(Surv(time5,status5)~1), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"We conclude for the Poher-Perme method, that in a world where cancer patients could only die due to cancer, only 41% of these patients would still be alive 5 year after their diagnosis.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"crude_pp = CrudeMortality(pp)\nprintln(crude_pp.Mₒ[1826], \" , \", crude_pp.Mₑ[1826], \" , \", crude_pp.Mₚ[1826])","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Finally, for this estimator, we have that of the 0.64 patients that have died, 0.53 is due to colorectal cancer while 0.11 is due to other causes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"We will plot the Pohar Perme method only.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"conf_int = confint(pp; level = 0.05)\nlower_bounds = [lower[1] for lower in conf_int]\nupper_bounds = [upper[2] for upper in conf_int] \n\np1 = plot(pp.grid, pp.Sₑ, ribbon=(pp.Sₑ - lower_bounds, upper_bounds - pp.Sₑ), xlab = \"Time (days)\", ylab = \"Net survival\", label = false)\n\np2 = plot(pp.grid, crude_pp.Mₑ, label = \"Excess Mortality Rate\")\np2 = plot!(pp.grid, crude_pp.Mₚ, label = \"Population Mortality Rate\")\n\nplot(p1,p2)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Looking at the graph, and the rapid dip it takes, it is evident that the first 5 years are crucial and that the survival probability is highly affected in these years. Additionnally, the crude mortality graph allows us to see how much of this curve is due to the colorectacl cancer studied versus other undefined causes. It is clear that the large majority is due to the cancer.","category":"page"},{"location":"example/#Net-survival-with-respect-to-covariates","page":"Case study","title":"Net survival with respect to covariates","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"We are now interested in comparing the different groups of patients defined by various covariates. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"pp_sex = fit(PoharPerme, @formula(Surv(time5,status5)~sex), colrec, slopop)\npp_males = pp_sex[pp_sex.sex .== :male,:estimator][1]\npp_females = pp_sex[pp_sex.sex .== :female,:estimator][1]","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"When comparing at time 1826, we notice that the survival probability is slightly inferior for men than for women (0433 0449). It is also more probable for the women to die from other causes than the men seeing as 00255 0025. Still, the differences are minimal. Let's confirm this with the Grafféo log-rank test:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"test_sex = fit(GraffeoTest, @formula(Surv(time5,status5)~sex), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The p-value is indeed above 005. We cannot reject the null hypothesis H_0 and thus we dismiss the differences between the two sexes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"As for the age, we will define two different groups: individuals aged 65 and above and those who are not.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"colrec.age65 .= ifelse.(colrec.age .>= 65*365.241, :old, :young)\npp_age65 = fit(PoharPerme, @formula(Surv(time5,status5)~age65), colrec, slopop)\npp_young = pp_age65[pp_age65.age65 .== :young, :estimator][1]\npp_old = pp_age65[pp_age65.age65 .== :old, :estimator][1]","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Here, the difference between the two is much more important. In the first group, the individuals are aged under 65 and at 5 years time, they have a 501% chance of survival. On the other hand, the individuals aged 65 and up have a 401% chance of survival. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"It is also worth noting that their chances of dying from other causes is higher than the younger group, given their age. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"When applying the Grafféo test, we get the results below:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"test_age65 = fit(GraffeoTest, @formula(Surv(time5,status5)~age65), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The p-value is well under 005, meaning we reject the H_0 hypothesis and must admit there are differences between the individuals aged 65 and above and the others.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"When plotting both we get:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"conf_int_men = confint(pp_males; level = 0.05)\nlower_bounds_men = [lower[1] for lower in conf_int_men]\nupper_bounds_men = [upper[2] for upper in conf_int_men] \n\nconf_int_women = confint(pp_females; level = 0.05)\nlower_bounds_women = [lower[1] for lower in conf_int_women]\nupper_bounds_women = [upper[2] for upper in conf_int_women]\n\nconf_int_under65 = confint(pp_young; level = 0.05)\nlower_bounds_under65 = [lower[1] for lower in conf_int_under65]\nupper_bounds_under65 = [upper[2] for upper in conf_int_under65] \n\nconf_int_65 = confint(pp_old; level = 0.05)\nlower_bounds_65 = [lower[1] for lower in conf_int_65]\nupper_bounds_65 = [upper[2] for upper in conf_int_65] \n\nplot1 = plot(pp_males.grid, pp_males.Sₑ, ribbon=(pp_males.Sₑ - lower_bounds_men, upper_bounds_men - pp_males.Sₑ), xlab = \"Time (days)\", ylab = \"Net survival\", label = \"men\")\nplot1 = plot!(pp_females.grid, pp_females.Sₑ, ribbon=(pp_females.Sₑ - lower_bounds_women, upper_bounds_women - pp_females.Sₑ), xlab = \"Time (days)\", ylab = \"Net survival\", label = \"women\")\n\nplot2 = plot(pp_young.grid, pp_young.Sₑ, ribbon=(pp_young.Sₑ - lower_bounds_under65, upper_bounds_under65 - pp_young.Sₑ), xlab = \"Time (days)\", ylab = \"Net survival\", label = \"Under 65\")\nplot2 = plot!(pp_old.grid, pp_old.Sₑ, ribbon=(pp_old.Sₑ - lower_bounds_65, upper_bounds_65 - pp_old.Sₑ), xlab = \"Time (days)\", ylab = \"Net survival\", label = \"65 and up\")\n\nplot(plot1, plot2, layout = (1, 2))","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Visually, it is almost immediately understood that there are no worthy differences between the two sexes whereas the age65 variable seems to play a big role.","category":"page"},{"location":"example/#Estimated-sample-size-and-life-expectancy","page":"Case study","title":"Estimated sample size and life expectancy","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"Given that the age group plays a significant role in the study, we will now estimate the sample size by yearly intervals in order to better compare the age groups.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"elt, ess = nessie(@formula(Surv(time,status)~age65), colrec, slopop)\nelt","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The expected life time for the younger patients is significatively higher than for older patients (24.78 years > 10.29 years).","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"hcat(ess[:,3]...)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Finally, the table above represents yearly expected sample sizes for both age groups under 65 and above, with the second column representing the latter. We can see that the sample size decreases for the older patients in a much more dramatic way than for the younger ages.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Unsurprisingly, we can thus conclude that age plays an important role in the study.","category":"page"}] +[{"location":"references/#References","page":"References","title":"References","text":"","category":"section"},{"location":"references/","page":"References","title":"References","text":"F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).\n\n\n\nF. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).\n\n\n\nT. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).\n\n\n\nM. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).\n\n\n\nN. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).\n\n\n\nT. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).\n\n\n\nP. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).\n\n\n\nM. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).\n\n\n\nH. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).\n\n\n\nK. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).\n\n\n\n","category":"page"},{"location":"references/","page":"References","title":"References","text":"T. Hakulinen and K. H. Abeywickrama. A computer program package for relative survival analysis. Computer programs in biomedicine 19, 197–207 (1985).\n\n\n\nT. Hakulinen and L. Tenkanen. Regression analysis of relative survival rates. Journal of the Royal Statistical Society Series C: Applied Statistics 36, 309–317 (1987).\n\n\n\n","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"CurrentModule = NetSurvival","category":"page"},{"location":"benches/#Benchmarking-results","page":"Benchmarking results","title":"Benchmarking results","text":"","category":"section"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"This page provides benchmark results of several standard net survival routines implemented in this package. Note that the runtime also depends on other packages, in particular on RateTables.jl. ","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"note: Take numbers displayed here carefully\nAll the following benchmarks are run on github action contunous integration platform, which is a very slow computing engine. Thus, these numbers may not represent your local performance. A locally ran version of these benchmarks in availiable on the github readme, and the below code blocks can be used to check performance on your own hardware.","category":"page"},{"location":"benches/#Benchmarks-w.r.t.-relsurv","page":"Benchmarking results","title":"Benchmarks w.r.t. relsurv","text":"","category":"section"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"This first set of benchmark compares standard functionalities with their implementation in relsurv. Below numbers gives runtime mulitpliers w.r.t. R::relsurv, computed on github action CI.","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"using RateTables, NetSurvival, RCall, DataFrames\n\nfunction test_surv(r_method,::Type{E}, stratified) where E\n if stratified\n jl = @timed fit(E, @formula(Surv(time,status)~sex), colrec, slopop)\n @rput r_method\n r = @timed R\"\"\"\n rez = relsurv::rs.surv(survival::Surv(time, stat) ~ sex, rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop, method = r_method, add.times=1:8149)\n \"\"\"\n else\n jl = @timed fit(E, @formula(Surv(time,status)~1), colrec, slopop)\n @rput r_method\n r = @timed R\"\"\"\n rez = relsurv::rs.surv(survival::Surv(time, stat) ~ 1, rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop, method = r_method, add.times=1:8149)\n \"\"\"\n end\n return r.time / jl.time\nend\n\nfunction test_graffeo(stratified)\n if stratified\n jl = @timed fit(GraffeoTest, @formula(Surv(time,status)~stage+Strata(sex)), colrec, slopop)\n r = @timed R\"\"\"\n rez = relsurv::rs.diff(survival::Surv(time, stat) ~ stage + survival::strata(sex), rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop)\n \"\"\"\n else \n jl = @timed fit(GraffeoTest, @formula(Surv(time,status)~stage), colrec, slopop)\n r = @timed R\"\"\"\n rez = relsurv::rs.diff(survival::Surv(time, stat) ~ stage, rmap=list(age = age, sex = sex, year = diag), data = relsurv::colrec, ratetable = relsurv::slopop)\n \"\"\"\n end\n return r.time / jl.time \nend\ntest_all(stratified) = [\n test_surv(\"pohar-perme\", PoharPerme, stratified),\n test_surv(\"ederer1\", EdererI, stratified),\n test_surv(\"ederer2\", EdererII, stratified),\n test_surv(\"hakulinen\", Hakulinen, stratified),\n test_graffeo(stratified),\n]\ntest_all() = DataFrame(\n Algorithm = [\"Pohar Perme\", \"EdererI\", \"EdererII\", \"Hakulinen\", \"Graffeo's LRT\"], \n unstratified = test_all(false), \n stratified = test_all(true)\n)\nrez = test_all()\nrez = test_all() # discard first run.\n\n# note: to obtain the pretty printing from the readme, you need to install PrettyTables.jl and do : \n# using PrettyTables\n# pretty_table(rez, backend = Val(:markdown))\n\nrez","category":"page"},{"location":"benches/#Benchmarking-across-time","page":"Benchmarking results","title":"Benchmarking across time","text":"","category":"section"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"The following charts provide a glimpse of NetSurvival.jl's performance along time, also ran on github CI: ","category":"page"},{"location":"benches/","page":"Benchmarking results","title":"Benchmarking results","text":"","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"CurrentModule = NetSurvival","category":"page"},{"location":"getting_started/#Getting-Started","page":"Getting Started","title":"Getting Started","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"In many population-based studies, the specific cause of death is unidentified, unreliable or even unavailable. Relative survival analysis addresses this scenario, previously unexplored in general survival analysis. Different methods were created with the aim to construct a consistant and reliable estimator for this purpose.","category":"page"},{"location":"getting_started/#Net-survival-settings","page":"Getting Started","title":"Net survival settings","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Note first that, for any positive random variable X, we will use extensively the following functions (that each fully characterize the distribution of X): ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Function symbol and definition Function Name\nS_X(t) = mathbb P(X t) Survival function\nLambda_X(t) = -ln S_X(t) Cumulative Hazard function\nlambda_X(t) = partial Lambda_X(t) Instantaneous hazard function","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Consider a study that consists of censored survival times from a specific cause. Such a study consists of several random objects: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Random variable Name Is it observed ?\nE \"Excess\" lifetime ✘\nP \"Population\" lifetime ✘\nO = E wedge P \"Overall\" lifetime ✘\nC \"Censoring\" time ✘\nmathbf D Vector of covariates ✔\nT = O wedge C Event time ✔\nDelta = mathbf1T leq C Event status ✔","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"It is important to note that we do not observe a a potential indicator mathbf1E geq P. This is one of the key differences between net survival and standard survival. The standard Net survival analysis solves this problem by assuming that the underlying times E and P are independent from each other.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"One central hypothesis in net survival (on top of non-informative censoring) is that P and E are independent. This independence can be written in terms of hazard rates lambda_O(t) = lambda_P(t) + lambda_E(t), or in terms of survival functions S_O(t) = S_P(t)S_E(t).","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The population hazard for each individual lambda_P_i is usally drawn from a reference life table, and may depend on covariates mathbf D_i such as age and date, sex, country, race, etc... See the RateTables.jl package for more details on the potential covariates. On the other hand, the excess mortality is assumed to be i.i.d. between individuals and not to depend on covariates at all. Thus, we mostly omit these covariates from our notations.","category":"page"},{"location":"getting_started/#Available-Estimators","page":"Getting Started","title":"Available Estimators","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The estimation of net survival is usually discussed in terms of the estimation of the cumulative excess hazard Lambda_E(t) and/or the instantaneous hazard lambda_E = partialLambda_E. To describe the estimators, we use the following counting processes notations, similar to standard survival analysis(see e.g. [6] or [7]). ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The uncensored event indicatrix partial N_i(t) for individual i at time t \nThe total number of uncensored events process partial N(t) = sum_i partial N_i(t) at time t\nThe at-risk indicatrix Y_i(t), for whether an individual is still at risk \nThe total number at risk process Y(t) = sum_i Y_i(t) at time t","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"With these definitions and assumptions in mind, we will now present the four different methods implemented in this package, commonly used in literature, to estimate the excess hazard function partialLambda_E(t) and its variance. Recall that we estimate the variance as sigma_E^2(t) = int_0^t partialsigma_E^2(s). ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Name Reference Proposed (partial) Excess Hazard partialhatLambda_E(s) Proposed (partial) Variance partialhatsigma_E^2(s)\nEderer I [1] fracsum_i N_i(s)sum_i Y_i(s) - fracsum_i S_P_i(s)partialLambda_P_i(s)sum_i S_P_i(s) fracsum_i N_i(s)left(sum_i Y_i(s)right)^2\nEderer II [2] fracsum_i N_i(s)sum_i Y_i(s) - fracsum_i Y_i(s)partialLambda_P_i(s)sum_i Y_i(s) fracsum_i N_i(s)left(sum_i Y_i(s)right)^2\nHakulinen [3] fracsum_i N_i(s)sum_i Y_i(s) - fracsum_i fracY_i(s) S_P_i(s)partialLambda_P_i(s)sum_i fracY_i(s) S_P_i(s) fracsum_i N_i(s)left(sum_i Y_i(s)right)^2\nPohar Perme [4] fracsum_i fracpartial N_i(s)S_P_i(s) - sum_i fracY_i(s)S_P_i(s)partialLambda_P_i(s)sum_i fracY_i(s)S_P_i(s) fracsum_i=1^n fracpartial N_i(s)S^2_P_ileft(sum_i fracY_i(s)S_p_i(s)right)^2","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where, in the variances, it is understood that when no more individuals are at risk 00 gives 0. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The Pohar Perme estimator [4] is the newest addition to relative survival analysis between the four methods, particularly designed to handle situations where covariates may change over time. It is trusted from the field (see e.g. [8] and [9]) that only this estimator should really be used, the other ones being included mostly for historical reasons and comparisons. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"PoharPerme\nEdererI\nEdererII\nHakulinen","category":"page"},{"location":"getting_started/#NetSurvival.PoharPerme","page":"Getting Started","title":"NetSurvival.PoharPerme","text":"PoharPerme\n\nThis method estimates net survival probabilities by applying the following estimation:\n\npartialhatLambda_E(t) = fracsum_i fracdN_i(u)S_P_i(u) - sum_i fracY_i(u)S_P_i(u)dLambda_P_i(u)sum_i fracY_i(u)S_P_i(u)\n\nTo fit the Pohar Perme to your data based on a certain rate table, apply the example below to your code : \n\nfit(PoharPerme, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\nReferences: \n\n[4] Perme, Maja Pohar and Stare, Janez and Estève, Jacques (2012). On Estimation in Relative Survival.\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#NetSurvival.EdererI","page":"Getting Started","title":"NetSurvival.EdererI","text":"EdererI\n\nThis method estimates net survival probabilities by applying the following estimation:\n\npartialhatLambda_E(t) = fracsum_i N_i(u)sum_i Y_i(u) - fracsum_i S_P_i(u)dLambda_P_i(u)sum_i S_P_i(u)\n\nTo call this function: \n\nfit(EdererI, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#NetSurvival.EdererII","page":"Getting Started","title":"NetSurvival.EdererII","text":"EdererII\n\nThis method estimates net survival probabilities by applying the following estimation:\n\npartialhatLambda_E(t) = fracsum_i N_i(u)sum_i Y_i(u) - fracsum_i Y_i(u)dLambda_P_i(u)sum_i Y_i(u)\n\nTo call this function: \n\nfit(EdererII, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#NetSurvival.Hakulinen","page":"Getting Started","title":"NetSurvival.Hakulinen","text":"Hakulinen\n\nTo call this function: \n\nfit(Hakulinen, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#Crude-Mortality","page":"Getting Started","title":"Crude Mortality","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The crude mortality rate is the global mortality rate (from all causes of death) for a population. This measure does not use the cause of death information which, as we previously mentioned, can be unreliable and incomplete. It is given as:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"M_E(t) = int_0^t S_O(u-) lambda_E(u)du","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"There exists a few estimators of this quantity, the most known one being the Cronin-Feuer estimator [10], given by:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"hatM_E(t) = int_0^t hatS_O(u-) hatlambda_E(u)du","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where, traditionally, Kaplan-Meier is used to estimate the overall survival function hatS_O and Ederer II is used for the excess hazard rate hatlambda_E and the population hazard rate hatlambda_P.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Our implementation is a bit more permissive, as any net survival estimators can be used for hatlambda_E. Of course, the default is still the original Ederer II which provides the original Cronin-Feuer estimator. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"CrudeMortality","category":"page"},{"location":"getting_started/#NetSurvival.CrudeMortality","page":"Getting Started","title":"NetSurvival.CrudeMortality","text":"CrudeMortality\n\nThis method calculates the crude mortality and presents both the excess mortality and population mortality rates.\n\nThe default Cronin-Feuer estimator can be fitted to data with the following interface: \n\nfit(CrudeMortality, args...)\n\nwhere the args are passed to fit(EdererII,args...) to compute the excess hazard.\n\nA more direct syntax can be used, specifying directly the estimator for the excess hazard:\n\nCrudeMortality(fit(EdererII,args...))\n\nTo use another excess hazard estimator, simply replace EdererII with the method of your choice (PoharPerme, EdererI, Hakulinen).\n\nReferences: \n\n[10] Cronin, Kathleen A and Feuer, Eric J (2000). Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival.\n\n\n\n\n\n","category":"type"},{"location":"getting_started/#Grafféo-Log-Rank-Test","page":"Getting Started","title":"Grafféo Log-Rank Test","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The Grafféo Log-Rank Test [5] was constructed as a complement to the Pohar Perme estimator, aiming to compare the net survival functions provided by the latter. The test is designed to compare these functions across multiple groups, including stratified covariables, and to ultimately determine, with the given p-value, which covariables are impactful to the study.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The null (H_0) hypothesis tests the following assumption:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"forall t in 0T Lambda_Eg_1(t) = Lambda_Eg_2(t) = = Lambda_Eg_k(t)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where G = g_1g_k is a partition of 1n consisting of disjoint groups of individuals that we wish to compare to each other. For all group g in G, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"partial N_Eg(s) = sum_i in g fracpartial N_i(s)S_P_i(s) - fracY_i(s)S_P_i(s)partialLambda_P_i(s)\nY_Eg(s) = sum_i in g fracY_i(s)S_P_i(s)\nR_g(s) = fracY_Eg(s)sum_gin G Y_Eg(s)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Then, define the vector mathbf Z = left(Z_g_r r in 1k-1 right) with entries: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Z_g(T) = N_Eg(s) - int_0^T Y_Eg(s) partialhatLambda_E(s)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The test statistic is then given by:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"U(T) = mathbf Z(T)hatSigma_Z^-1 mathbf Z(T)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"where the entries of the hatSigma_Z matrix are given by: ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"sigma_gh(T) = int_0^T sum_ell in G left(delta_gell - R_g(t) right)left(delta_hell - R_h(t)right) left(sum_iinell fracpartial N_i(s)S^2_P_iright)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"Under H_0, the statistic U(T) is asymptotically chi^2(k-1)-distributed. We thus reject the H_0 hypothesis when the p-value obtained is under 005, admitting the notable difference between the groups. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"GraffeoTest","category":"page"},{"location":"getting_started/#NetSurvival.GraffeoTest","page":"Getting Started","title":"NetSurvival.GraffeoTest","text":"GraffeoTest\n\nThe Grafféo test is a log-rank type test and is typically used in net survival analysis to determine the impact of certain covariates in the study.\n\nThe null (H_0) hypothesis tests the following assumption:\n\nforall t in 0T Lambda_Eg_1(t) = Lambda_Eg_2(t) = = Lambda_Eg_k(t)\n\nwhere G = g_1g_k is a partition of 1n consisting of disjoint groups of individuals that we wish to compare to each other. For all group g in G, let's denote the numerator and denominator of the Pohar Perme (partial) excess hazard estimators, restricted to individuals in the group, by: \n\npartial N_Eg(s) = sum_i in g fracpartial N_i(s)S_P_i(s) - fracY_i(s)S_P_i(s)partialLambda_P_i(s)\nY_Eg(s) = sum_i in g fracY_i(s)S_P_i(s)\nR_g(s) = fracY_Eg(s)sum_gin G Y_Eg(s)\n\nThen, define the vector mathbf Z = left(Z_g_r r in 1k-1 right) with entries: \n\nZ_g(T) = N_Eg(s) - int_0^T Y_Eg(s) partialhatLambda_E(s)\n\nThe test statistic is then given by:\n\nU(T) = mathbf Z(T)hatSigma_Z^-1 mathbf Z(T)\n\nwhere the entries of the hatSigma_Z matrix are given by: \n\nsigma_gh(T) = int_0^T sum_ell in G left(delta_gell - R_g(t) right)left(delta_hell - R_h(t)right) left(sum_iinell fracpartial N_i(s)S^2_P_iright)\n\nUnder H_0, the statistic U(T) is asymptotically chi^2(k-1)-distributed.\n\nTo apply the test to your data based on a certain rate table, apply the example below to your code : \n\nfit(GraffeoTest, @formula(Surv(time,status)~covariable1 + covariable2), data, ratetable)\n\nIf you wish to stratify a covariate:\n\nfit(GraffeoTest, @formula(Surv(time,status)~covariable1 + Strata(covariable2)), data, ratetable)\n\nThe produced test statistics is supposed to follow a chi squared distribution under (H_0). You can fetch the results using the .stat, .df and .pval fields of the returned object. \n\nReferences: \n\n[5] Grafféo, Nathalie and Castell, Fabienne and Belot, Aurélien and Giorgi, Roch (2016). A Log-Rank-Type Test to Compare Net Survival Distributions. \n\n\n\n\n\n","category":"type"},{"location":"getting_started/#Nessie","page":"Getting Started","title":"Nessie","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The Nessie function estimates the sample size by yearly intervals as well as averages an estimated lifespan left for a given group. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"This function is highly dependant on the Life function taken from the RateTables.jl package which you can find documented here.","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"The sample size is thus taken by the following formula:","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"ESS(t) = sum_i^N S_P_i(t)","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"While the estimated lifepsan is directly taken from the expectation function. ","category":"page"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"nessie","category":"page"},{"location":"getting_started/#NetSurvival.nessie","page":"Getting Started","title":"NetSurvival.nessie","text":"nessie\n\nTo call this function, use the formula below: \n\nnessie(@formula(Surv(time,status)~covariate), data, ratetable)\n\n\n\n\n\n","category":"function"},{"location":"getting_started/#References","page":"Getting Started","title":"References","text":"","category":"section"},{"location":"getting_started/","page":"Getting Started","title":"Getting Started","text":"F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).\n\n\n\nF. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).\n\n\n\nT. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).\n\n\n\nM. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).\n\n\n\nN. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).\n\n\n\nT. R. Fleming and D. P. Harrington. Counting Processes and Survival Analysis. Vol. 625 (John Wiley & Sons, 2013).\n\n\n\nP. K. Andersen, O. Borgan, R. D. Gill and N. Keiding. Statistical Models Based on Counting Processes. Springer Series in Statistics (Springer US, New York, NY, 1993).\n\n\n\nM. P. Perme and K. Pavlic. Nonparametric relative survival analysis with the R package relsurv. Journal of Statistical Software 87, 1–27 (2018).\n\n\n\nH. Charvat and A. Belot. Mexhaz: An R package for fitting flexible hazard-based regression models for overall and excess mortality with a random effect. Journal of Statistical Software 98, 1–36 (2021).\n\n\n\nK. A. Cronin and E. J. Feuer. Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival. Statistics in medicine 19, 1729–1740 (2000).\n\n\n\n","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = NetSurvival","category":"page"},{"location":"#Introduction","page":"Home","title":"Introduction","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"The NetSurvival.jl package provides the necessary tools to perform estimations and analysis in the Net Survival field. This specialized branch of Survival Analysis focuses on estimating the probability of survival from a specific event of interest, for example a given cancer, without considering other causes of death. This is especially relevant in the (unfortunately quite common) case where the cause of death indicatrix is either unavailable or untrustworthy. Consequently, the so-called missing indicatrix issue forbids the use of standard competitive risks survival analysis methods on these datasets. For that, a few standard estimators were established in the last 50 years, backed by a wide literature.","category":"page"},{"location":"","page":"Home","title":"Home","text":"By integrating observed data from the target population with historical population mortality data (usually sourced from national census datasets), Net Survival allows the extraction of the specific mortality hazard associated with the particular disease, even under the missing indicatrix issue. The concept of relative survival analysis dates back several decades to the seminal article by Ederer, Axtell, and Cutler in 1961 [1] and the one by Ederer and Heise in 1959 [2].","category":"page"},{"location":"","page":"Home","title":"Home","text":"For years, the Hakulinen estimator (1977) [3] and the Ederer I and II estimators were widely regarded as the gold standard for non-parametric survival curve estimation. However, the introduction of the Pohar Perme, Stare, and Estève estimator in 2012 [4] resolved several issues inherent in previous estimators, providing a reliable and consistent non-parametric estimator for net survival analysis.","category":"page"},{"location":"#Features","page":"Home","title":"Features","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Standard tools nowadays are composed of R packages, with underlying C and C++ routines, that are hard to read, maintain, and use. This package is an attempt to bring standard relative survival analysis modeling routines to Julia, while providing an interface that is close to the relsurv standard, albeit significantly faster and easier to maintain in the future. Our hope is that the junction with classical modeling API in Julia will allow later extensions of the existing modeling methods, with a simple interface for the practitioners.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Some key features in NetSurvival.jl are:","category":"page"},{"location":"","page":"Home","title":"Home","text":"A panel of different non-parametric net survival estimators (Ederer I [1], Ederer II [2], Hakulinen [3], Pohar Perme [4]) with an interface compliant with Julia's standards. \nGrafféo's log-rank test [5] to compare net survival curves accross groups, including stratified testing.\nCrude mortality, Expected Sample Size, and other useful metrics in net survival field.\nA 'Nessie' function that outputs the estimated sample size by yearly intervals and the average lifespan expectancy left for a given group. \nA compact, readable and efficient codebase (up to 1000x less LOC than relsurv for the same functionalities), ensuring long-term maintenability.\nSignificant performance improvements (up to 50x) compared to the R package relsurv.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"The package is not yet available on Julia's general registry, and thus can be installed through the following command:","category":"page"},{"location":"","page":"Home","title":"Home","text":"using Pkg\nPkg.add(\"https://github.com/JuliaSurv/NetSurvival.jl.git\")","category":"page"},{"location":"","page":"Home","title":"Home","text":"See the rest of this documentation to have a glimpse of the functionalities!","category":"page"},{"location":"#References","page":"Home","title":"References","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"F. Ederer. The relative survival rate: a statistical methodology. Natl. Cancer Inst. Monogr. 6, 101–121 (1961).\n\n\n\nF. Ederer and H. Heise. The effect of eliminating deaths from cancer in general survival rates, methodological notes 11. End Result Evaluation Section, National Cancer Institute (1959).\n\n\n\nT. Hakulinen. On long-term relative survival rates. Journal of Chronic Diseases 30, 431–443 (1977).\n\n\n\nM. P. Perme, J. Stare and J. Estève. On Estimation in Relative Survival. Biometrics 68, 113–120 (2011).\n\n\n\nN. Grafféo, F. Castell, A. Belot and R. Giorgi. A Log-Rank-Type Test to Compare Net Survival Distributions. Biometrics 72, 760–769 (2016).\n\n\n\n","category":"page"},{"location":"example/#Case-study","page":"Case study","title":"Case study","text":"","category":"section"},{"location":"example/#Introduction-and-datasets","page":"Case study","title":"Introduction and datasets","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"We will illustrate with an example using the dataset colrec, which comprises 5971 patients diagnosed with colon or rectal cancer between 1994 and 2000. This dataset is sourced from the Slovenia cancer registry. Given the high probability that the patients are Slovenian, we will be using the Slovenian mortality table slopop as reference for the population mortality rates. Subsequently, we can apply various non-parametric estimators for net survival analysis.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"note: N.B.\nMortality tables may vary in structure, with options such as the addition or removal of specific covariates. To confirm that the mortality table is in the correct format, please refer to RateTables.jl's documentation, or directly extract it from there.","category":"page"},{"location":"example/#Cohort-details","page":"Case study","title":"Cohort details","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"The patients in the study are diagnosed between January 1st 1994 and December 31st 2000. Before we move on to the survival probabilities, it is important to be aware of how your data is distributed and of what it comprises. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"using NetSurvival, RateTables, DataFrames\nfirst(colrec,10)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Let's explore how our data is distributed first, starting with the age variable.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"println(minimum(colrec.age./365.241), \" ; \",maximum(colrec.age./365.241))","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The study can be considered diverse in terms of age seeing as the patients are between 12 and 96 years old, approximately. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"using Plots \nplot(histogram(colrec.age./365.241, label=\"Age\"),\n histogram(colrec.time./365.241, label=\"Follow-up time\"))","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The graph above show us that although the dataset has a wide range of patients within all age groups, it is mostly centered around older adults and elderly, with the majority of the patients being between 60 and 80 years old. Looking at the second graph that details the distribution of the follow-up times. We notice there that the values quickly drop. Unfortunately, this is a common theme in cancer studies. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Let's take a look at the sex variable now, by looking at the number of male and female patients:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"combine(groupby(colrec, :sex), nrow)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"There isn't too big of a difference between the two. We can say this study includes both gender relatively equally, thus, reducing bias. With these two observations, it is also worth noting that colorectal cancer is most common with men and people older than 50.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"In total, we note that we have 5971 patients. By taking a look at the status variable, we can determine the deaths and censorship:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"sum(colrec.status)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Out of the 5971 patients, 4979 have died. This is a very high mortality rate, and again, unfortunately common in cancer studies.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"(nrow(colrec) - sum(colrec.status)) / nrow(colrec)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"In other terms, the censorship rate is of 16.6%, meaning the event, in this case death, was not observed for only 16.6% of the patients. This is a low censorship rate and thus the quality of the signal will be pretty good. ","category":"page"},{"location":"example/#Mortality-table","page":"Case study","title":"Mortality table","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"We will be using the mortality table slopop taken from the RateTables.jl package as the study is done on Slovenian patients. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"slopop","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The show method of the RateTable class shows the additional covariate sex that the rate table has on top of the (mandatory) age and year variables. The sex variable has two madalities, :male and :female. The ratetable is then three dimensional. For example, the daily hazard rate for a woman turning 45 on the January 1st 2006 can be accessed through the following command:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"λ = daily_hazard(slopop, 45*365.241, 2006*365.241; sex=:female)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Making the daily survival probability easily calculated with:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"exp(-λ)","category":"page"},{"location":"example/#Overall-and-expected-survival","page":"Case study","title":"Overall and expected survival","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"For this part, we will be using the Survival.jl package to apply the Kaplan Meier estimator for the overall survival.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"using Survival \nkm = fit(KaplanMeier, colrec.time./365.241, colrec.status)\nplot(km.events.time, km.survival, label=false, title = \"Kaplan-Meier Estimator for the Overall Survival\")","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The graph above indicates a significant dip in survival probability within the first 5 years, and even more at 10 years. This period in the study is then crucial for the analysis. ","category":"page"},{"location":"example/#Estimated-net-survival","page":"Case study","title":"Estimated net survival","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"We will restrict ourselves to the first 5 years of the study. For that, let us re-censor the dataset as follows: ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"for i in 1:nrow(colrec)\n if colrec.time[i] > 1826 # five years\n colrec.status[i] = false\n colrec.time[i] = 1826\n end\nend","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"We can now apply the different non-parametric methods to compute the relative survival.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"e1 = fit(EdererI, @formula(Surv(time,status)~1), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"With the EdererI method, after 1826 days have passed, we can say that the survival rate at this mark is around 0456, in the hypothetical world where patients can only die of cancer.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"crude_e1 = CrudeMortality(e1)\nprintln(crude_e1.Mₒ[1826], \" , \", crude_e1.Mₑ[1826], \" , \", crude_e1.Mₚ[1826])","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Out of the 0.63 patients that have died, according to the EdererI method, 0.51 died because of colorectal cancer and 0.12 died of other causes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"e2 = fit(EdererII, @formula(Surv(time,status)~1), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Similarly, the EdererII method, also known as the conditional method, shows that at the 5 year mark, the survival probability is of 044 in this hypothetical world.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"crude_e2 = CrudeMortality(e2)\nprintln(crude_e2.Mₒ[1826], \" , \", crude_e2.Mₑ[1826], \" , \", crude_e2.Mₚ[1826])","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Here, out of the 0.63 patients that have died, 0.53 are due to colorectal cancer and 0.1 due to other causes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"pp = fit(PoharPerme, @formula(Surv(time,status)~1), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"We conclude for the Pohar Perme method, that in a world where cancer patients could only die due to cancer, only 41% of these patients would still be alive 5 year after their diagnosis. The Pohar Perme estimator is the best estimator of the excess hazard under the standard hypotheses. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"crude_pp = CrudeMortality(pp)\nprintln(crude_pp.Mₒ[1826], \" , \", crude_pp.Mₑ[1826], \" , \", crude_pp.Mₚ[1826])","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Finally, for this estimator, we have that of the 0.64 patients that have died, 0.53 is due to colorectal cancer while 0.11 is due to other causes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"We will plot the Pohar Perme method only.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"function mkribbon(pp)\n S = pp.Sₑ\n ci = confint(pp; level = 0.05)\n l,u = getindex.(ci, 1), getindex.(ci,2)\n rb = (S - l, u - S)\n return rb\nend\n\np1 = plot(pp.grid, pp.Sₑ, ribbon=mkribbon(pp), xlab = \"Time (days)\", ylab = \"Net survival\", label = false)\n\np2 = plot(pp.grid, crude_pp.Mₑ, label = \"Excess Mortality Rate\")\np2 = plot!(pp.grid, crude_pp.Mₚ, label = \"Population Mortality Rate\")\n\nplot(p1,p2)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Looking at the graph, and the rapid dip it takes, it is evident that the first 5 years are crucial and that the survival probability is highly affected in these years. Additionally, the crude mortality graph allows us to see how much of this curve is due to the colorectal cancer studied versus other undefined causes. It is clear that the large majority is due to the cancer.","category":"page"},{"location":"example/#Net-survival-with-respect-to-covariates","page":"Case study","title":"Net survival with respect to covariates","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"We are now interested in comparing the different groups of patients defined by various covariates. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"pp_sex = fit(PoharPerme, @formula(Surv(time,status)~sex), colrec, slopop)\npp_males = pp_sex[pp_sex.sex .== :male,:estimator][1]\npp_females = pp_sex[pp_sex.sex .== :female,:estimator][1]","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"When comparing at time 1826, we notice that the survival probability is slightly inferior for men than for women (0433 0449). It is also more probable for the women to die from other causes than the men seeing as 00255 0025. Still, the differences are minimal. Let's confirm this with the Grafféo log-rank test:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"test_sex = fit(GraffeoTest, @formula(Surv(time,status)~sex), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The p-value is indeed above 005. We cannot reject the null hypothesis H_0 and thus we dismiss the differences between the two sexes.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"As for the age, we will define two different groups: individuals aged 65 and above and those who are not.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"colrec.age65 .= ifelse.(colrec.age .>= 65*365.241, :old, :young)\npp_age65 = fit(PoharPerme, @formula(Surv(time,status)~age65), colrec, slopop)\npp_young = pp_age65[pp_age65.age65 .== :young, :estimator][1]\npp_old = pp_age65[pp_age65.age65 .== :old, :estimator][1]","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Here, the difference between the two is much more important. In the first group, the individuals are aged under 65 and at 5 years time, they have a 501% chance of survival. On the other hand, the individuals aged 65 and up have a 401% chance of survival. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"It is also worth noting that their chances of dying from other causes is higher than the younger group, given their age. ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"When applying the Grafféo test, we get the results below:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"test_age65 = fit(GraffeoTest, @formula(Surv(time,status)~age65), colrec, slopop)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The p-value is well under 005, meaning we reject the H_0 hypothesis and must admit there are differences between the individuals aged 65 and above and the others.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"When plotting both we get:","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"plot1 = plot(pp_males.grid, pp_males.Sₑ, ribbon=mkribbon(pp_males), xlab = \"Time (days)\", ylab = \"Net survival\", label = \"men\")\nplot1 = plot!(plot1, pp_females.grid, pp_females.Sₑ, ribbon=mkribbon(pp_females), label = \"women\")\n\nplot2 = plot(pp_young.grid, pp_young.Sₑ, ribbon=mkribbon(pp_young), xlab = \"Time (days)\", ylab = \"Net survival\", label = \"Under 65\")\nplot2 = plot!(pp_old.grid, pp_old.Sₑ, ribbon=mkribbon(pp_old), label = \"65 and up\")\n\nplot(plot1, plot2, layout = (1, 2))","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Visually, it is almost immediately understood that there are no worthy differences between the two sexes whereas the age65 variable seems to play a big role.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The same kind of graph can be made on the stage: ","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"pp3 = fit(PoharPerme, @formula(Surv(time,status)~stage), colrec, slopop)\nplot3 = plot(xlab = \"Time (days)\", ylab = \"Net survival\",title=\"Net survival per cancer stages\")\nfor i in 1:nrow(pp3)\n e = pp3[i,:estimator]\n plot!(plot3, e.grid, e.Sₑ, ribbon=mkribbon(e), label=pp3[i,:stage])\nend\nplot(plot3)","category":"page"},{"location":"example/#Estimated-sample-size-and-life-expectancy","page":"Case study","title":"Estimated sample size and life expectancy","text":"","category":"section"},{"location":"example/","page":"Case study","title":"Case study","text":"Given that the age group plays a significant role in the study, we will now estimate the sample size by yearly intervals in order to better compare the age groups.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"elt, ess = nessie(@formula(Surv(time,status)~age65), colrec, slopop)\nelt","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"The expected life time for the younger patients is significatively higher than for older patients (24.78 years > 10.29 years).","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"hcat(ess[:,3]...)","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Finally, the table above represents yearly expected sample sizes for both age groups under 65 and above, with the second column representing the latter. We can see that the sample size decreases for the older patients in a much more dramatic way than for the younger ages.","category":"page"},{"location":"example/","page":"Case study","title":"Case study","text":"Unsurprisingly, we can thus conclude that age plays an important role in the study.","category":"page"}] }