Survival analysis of
Skin Cutaneous Melanoma[SKCM]

A case study of The Cancer Genome Atlas，TCGA

作者

Kili

发布于

2024-08-14

摘要

癌症基因组图谱计划(The Cancer Genome,Atlas,TCGA)是美国的国家肿瘤研究所（National Cancer Institute，NCI）与国家人类基因研究所（National Human Genome Research Institute，NHGRI）共同开展的一个广泛而综合的合作计划。其目标是获取、刻画并分析人类癌症中大规模、多种变异的分子特征，且为癌症研究迅速地提供数据。选取发病于皮肤的皮肤黑色素瘤样本[SKCM],即Skin Cutaneous Melanoma作为生存分析对象.分析不同样本特征下的死亡率差别.

代码
```{r}
library(survival)
library(survminer)
library(autoReg)
library(readr)
library(survMisc)
library(easystats)
library(tidyverse)
library(gtsummary)
library(rms)
library(conflicted)
library(ggplot2)
library(scales)
conflicts_prefer(dplyr::filter)
conflicts_prefer(dplyr::summarize)
conflicts_prefer(dplyr::select)
skcm <- read_delim("data/skcm_new.csv",
  delim = "\t", escape_double = FALSE,
  col_types = cols(form_completion_date = col_date(format = "%Y-%m-%d"), birth_days_to = col_number()),
  trim_ws = TRUE
)
```
代码
```{r}
sum_of_na <- function(x) {
  sum(is.na(x))
}
skcm %>% summarise(
  across(everything(), sum_of_na)
)
```
# A tibble: 1 × 58
  bcr_patient_barcode bcr_patient_uuid form_completion_date  days status
                <int>            <int>                <int> <int>  <int>
1                   0                0                    0     8      0
# ℹ 53 more variables: birth_days_to <int>, height_cm_at_diagnosis <int>,
#   weight_kg_at_diagnosis <int>, primary_at_dx_count <int>,
#   breslow_thickness_at_diagnosis <int>, primary_melanoma_mitotic_rate <int>,
#   initial_pathologic_dx_year <int>, age_at_diagnosis <int>,
#   submitted_tumor_dx_days_to <int>, new_tumor_event_melanoma_count <int>,
#   days_to_initial_pathologic_diagnosis <int>, prospective_collection <int>,
#   retrospective_collection <int>, gender <int>, race <int>, …

1 数据背景

Melanoma is a cancer in a type of skin cells called melanocytes. Melanocyes are the cells that produce melanin, which colors the skin. When exposed to sun, these cells make more melanin, causing the skin to darken or tan. Melanoma can occur anywhere on the body and risk factors include fair complexion, family history of melanoma, and being exposed to natural or artificial sunlight over long periods of time. Melanoma is most often discovered because it has metastasized, or spread, to another organ, such as the lymph nodes. In many cases, the primary skin melanoma site is never found. Because of this challenge, TCGA is studying primarily metastatic cases (in contrast to other cancers selected for study, where metastatic cases are excluded). For 2011, it was estimated that there were 70,230 new cases of melanoma and 8,790 deaths from the disease

2 数据清洗与数据预览

先预览数据:

代码
```{r}
#| tbl-cap: 数据预览
skcm <- skcm |>
  select(!bcr_patient_barcode:bcr_patient_uuid) |>
  select(!submitted_tumor_dx_days_to)
skcm <- skcm |> dplyr::select(days, status, height_cm_at_diagnosis, weight_kg_at_diagnosis, breslow_thickness_at_diagnosis, primary_melanoma_mitotic_rate, age_at_diagnosis, prospective_collection, gender, tumor_status, clark_level_at_diagnosis, primary_melanoma_tumor_ulceration, ajcc_tumor_pathologic_pt, ajcc_nodes_pathologic_pn, submitted_tumor_site.1, radiation_treatment_adjuvant, new_tumor_event_dx_indicator)
skcm <- skcm %>%
  mutate(across(function(col) n_distinct(col) < 10, factor)) %>%
  as_tibble()

gtsummary::tbl_summary(skcm, by = "tumor_status") %>%
  as_flex_table() |>
  flextable::set_caption("Summary statistics.by~tumor_status") %>%
  flextable::theme_zebra()
```

显然这是右删失数据,status=1代表死亡,=2代表右删失

观察到有缺失值,进行删除与补全

代码
```{r}
#| output: false
skcm <- skcm |>
  dplyr::filter(days > 0) |>
  filter(!(is.na(radiation_treatment_adjuvant) | is.na(new_tumor_event_dx_indicator) | is.na(submitted_tumor_site.1) | is.na(submitted_tumor_site.1) | is.na(ajcc_nodes_pathologic_pn))) |>
  filter(!is.na(tumor_status))
skcm <- skcm |> filter(!is.na(ajcc_tumor_pathologic_pt))
skcm$height_cm_at_diagnosis[is.na(skcm$height_cm_at_diagnosis)] <- mean(na.omit(skcm$height_cm_at_diagnosis))
skcm$weight_kg_at_diagnosis[is.na(skcm$weight_kg_at_diagnosis)] <- mean(na.omit(skcm$weight_kg_at_diagnosis))
skcm$breslow_thickness_at_diagnosis[is.na(skcm$breslow_thickness_at_diagnosis)] <- mean(na.omit(skcm$breslow_thickness_at_diagnosis))
library(mice)
tempData <- mice(skcm, m = 5, maxit = 50, seed = 500)
completedData <- complete(tempData, action = 5)
skcm <- complete(tempData, action = 5)
```

将高度除以重量创造新特征BMI:

代码
```{r}
skcm <- skcm |>
  mutate(BMI = height_cm_at_diagnosis / weight_kg_at_diagnosis) |>
  select(-height_cm_at_diagnosis, -weight_kg_at_diagnosis)
```
代码
```{r}
head(skcm)
gaze(skcm)
```
  days status breslow_thickness_at_diagnosis primary_melanoma_mitotic_rate
1  387      0                           13.0                             4
2   14      0                            9.0                            13
3  282      1                           12.0                            10
4   12      0                            8.0                             9
5   17      0                            5.0                            10
6  504      1                            0.4                             3
  age_at_diagnosis prospective_collection gender tumor_status
1               46                      2      1            1
2               74                      2      1            1
3               56                      2      2            1
4               71                      2      1            1
5               80                      2      2            1
6               79                      1      1            2
  clark_level_at_diagnosis primary_melanoma_tumor_ulceration
1                        3                                 2
2                        4                                 2
3                        3                                 2
4                        3                                 1
5                        3                                 2
6                        2                                 1
  ajcc_tumor_pathologic_pt ajcc_nodes_pathologic_pn submitted_tumor_site.1
1                       13                        1                      2
2                       13                        1                      2
3                       13                        6                      2
4                       12                        1                      2
5                       13                        1                      2
6                        4                        9                      4
  radiation_treatment_adjuvant new_tumor_event_dx_indicator      BMI
1                            1                            2 2.758621
2                            1                            1 2.285714
3                            1                            1 2.243590
4                            1                            1 2.910714
5                            1                            1 2.507246
6                            1                            2 2.745763
———————————————————————————————————————————————————————————————
               name                  levels         stats      
———————————————————————————————————————————————————————————————
days                               Mean ± SD    1741.0 ± 1883.8 
status                             0                218 (61.4%) 
                                   1                137 (38.6%) 
breslow_thickness_at_diagnosis     Mean ± SD          5.0 ± 6.4 
primary_melanoma_mitotic_rate      Mean ± SD          7.2 ± 7.1 
age_at_diagnosis                   Mean ± SD        56.9 ± 15.7 
prospective_collection             1                238 (67.0%) 
                                   2                117 (33.0%) 
gender                             1                129 (36.3%) 
                                   2                226 (63.7%) 
tumor_status                       1                156 (43.9%) 
                                   2                199 (56.1%) 
clark_level_at_diagnosis           1                  14 (3.9%) 
                                   2                  18 (5.1%) 
                                   3                 86 (24.2%) 
                                   4                170 (47.9%) 
                                   5                 67 (18.9%) 
primary_melanoma_tumor_ulceration  1                190 (53.5%) 
                                   2                165 (46.5%) 
ajcc_tumor_pathologic_pt           Mean ± SD          9.1 ± 4.3 
ajcc_nodes_pathologic_pn           Mean ± SD          3.9 ± 3.4 
submitted_tumor_site.1             1                 49 (13.8%) 
                                   2                 58 (16.3%) 
                                   3                 64 (18.0%) 
                                   4                184 (51.8%) 
radiation_treatment_adjuvant       1                338 (95.2%) 
                                   2                  17 (4.8%) 
new_tumor_event_dx_indicator       1                 82 (23.1%) 
                                   2                273 (76.9%) 
BMI                                Mean ± SD          2.1 ± 0.3 
———————————————————————————————————————————————————————————————
代码
```{r}
#| fig-cap: 小提琴图
skcm |>
  filter(status == 1) |>
  ggplot(aes(x = gender, y = days, fill = gender)) +
  geom_violin()
```

小提琴图

3 预测模型

将数据集分为train与test:

代码
```{r}
library(rsample)
set.seed(101)
skcm_split <- initial_split(
  skcm,
  prop = 0.90
)
train <- training(skcm_split)
test <- testing(skcm_split)
test.y <- test$class
```

3.1 Kaplan-Meier估计

代码
```{r}
library(ggsurvfit)
library(tidycmprsk)
skcm$status <- as.numeric(skcm$status)
s <- Surv(skcm$days, skcm$status)
model1 <- survfit2(Surv(days, status) ~ 1, skcm)
model1 |> ggsurvfit() +
  labs(
    x = "Days",
    y = "Overall survival probability"
  ) +
  add_confidence_interval() +
  add_risktable()
```

风险函数曲线估计与平滑化:

代码
```{r}
library(muhaz)
result.pe1 <- pehaz(skcm$days, skcm$status, width = 100, max.time = 9000)
plot(result.pe1)
result.smooth <- muhaz(skcm$days, skcm$status,
  bw.smooth = 100,
  b.cor = "left", max.time = 9000
)
lines(result.smooth)
```

max.time= 9000
width= 100
nbins= 90

or:

代码
```{r}
p1 <- epiR::epi.insthaz(model1) %>%
  ggplot(aes(x = time, y = hest)) +
  geom_smooth(color = "red", method = "loess", formula = "y ~ x") +
  theme_light() +
  labs(
    title = "Kaplan-Meier Hazard Function Estimate",
    x = "Time", y = "Instantaneous Hazard"
  )
p1
```

3.2 log-rank test

代码
```{r}
survdiff(Surv(days, status) ~ strata(gender) + clark_level_at_diagnosis, data = skcm)
```
Call:
survdiff(formula = Surv(days, status) ~ strata(gender) + clark_level_at_diagnosis, 
    data = skcm)

                             N Observed Expected (O-E)^2/E (O-E)^2/V
clark_level_at_diagnosis=1  14        9     6.15    1.3190     1.405
clark_level_at_diagnosis=2  18        4     8.64    2.4930     2.678
clark_level_at_diagnosis=3  86       30    44.78    4.8783     7.581
clark_level_at_diagnosis=4 170       65    62.86    0.0728     0.139
clark_level_at_diagnosis=5  67       29    14.57   14.3006    16.476

 Chisq= 23.9  on 4 degrees of freedom, p= 8e-05

p值小,which is statistically significant at the 5 % level.

代码
```{r}
library(ggsurvfit)
library(tidycmprsk)
skcm$status <- as.numeric(skcm$status)
s <- Surv(skcm$days, skcm$status)
model1 <- survfit2(Surv(days, status) ~ tumor_status, skcm)
model1 |> ggsurvfit() +
  labs(
    x = "Days",
    y = "Overall survival probability"
  ) + labs(title = "tumor_status")
```

3.3 parametric mdoel

代码
```{r}
library(flexsurv)
par_fits <- tibble(
  dist_param = c(
    "exp", "weibull", "gompertz", "gamma", "lognormal", "llogis",
    "gengamma"
  ),
  dist_name = c(
    "Exponential", "Weibull (AFT)", "Gompertz", "Gamma",
    "Lognormal", "Log-logistic", "Generalized gamma"
  )
) %>%
  mutate(
    fit = map(dist_param, ~ flexsurvreg(Surv(time, status) ~ 1, data = df_lung, dist = .x)),
    fit_smry = map(fit, ~ summary(.x, type = "hazard", ci = FALSE, tidy = TRUE)),
    AIC = map_dbl(fit, ~ .x$AIC)
  )

p2 <- par_fits %>%
  select(-c(dist_param, fit)) %>%
  unnest(fit_smry) %>%
  ggplot(aes(x = time, y = est, color = dist_name)) +
  geom_line() +
  theme_light() +
  labs(title = "Parametric Distribution Fits to Lung Cancer Data.")

p1 + p2
```

代码
```{r}
par_fits %>%
  arrange(AIC) %>%
  select(dist_name, AIC)
```
# A tibble: 7 × 2
  dist_name           AIC
  <chr>             <dbl>
1 Weibull (AFT)     1185.
2 Generalized gamma 1186.
3 Gamma             1186.
4 Gompertz          1188.
5 Log-logistic      1199.
6 Exponential       1200.
7 Lognormal         1215.

比较AIC,最好的是Weibullmodel

代码
```{r}
weibull_fit <- flexsurvreg(Surv(days, status) ~ age_at_diagnosis + tumor_status + clark_level_at_diagnosis + submitted_tumor_site.1, data = skcm, dist = "weibull")
weibull_fit
```
Call:
flexsurvreg(formula = Surv(days, status) ~ age_at_diagnosis + 
    tumor_status + clark_level_at_diagnosis + submitted_tumor_site.1, 
    data = skcm, dist = "weibull")

Estimates: 
                           data mean  est        L95%       U95%     
shape                             NA   1.24e+00   1.10e+00   1.41e+00
scale                             NA   6.25e+04   2.30e+04   1.70e+05
age_at_diagnosis            5.69e+01  -2.15e-02  -3.20e-02  -1.10e-02
tumor_status2               5.61e-01  -2.20e+00  -2.86e+00  -1.53e+00
clark_level_at_diagnosis2   5.07e-02   4.68e-01  -5.13e-01   1.45e+00
clark_level_at_diagnosis3   2.42e-01   5.01e-01  -1.28e-01   1.13e+00
clark_level_at_diagnosis4   4.79e-01   2.68e-01  -3.44e-01   8.79e-01
clark_level_at_diagnosis5   1.89e-01  -2.96e-01  -9.86e-01   3.93e-01
submitted_tumor_site.12     1.63e-01  -1.02e+00  -1.66e+00  -3.75e-01
submitted_tumor_site.13     1.80e-01   3.36e-01  -1.35e-01   8.06e-01
submitted_tumor_site.14     5.18e-01  -4.86e-02  -4.13e-01   3.16e-01
                           se         exp(est)   L95%       U95%     
shape                       7.82e-02         NA         NA         NA
scale                       3.20e+04         NA         NA         NA
age_at_diagnosis            5.34e-03   9.79e-01   9.69e-01   9.89e-01
tumor_status2               3.38e-01   1.11e-01   5.74e-02   2.16e-01
clark_level_at_diagnosis2   5.00e-01   1.60e+00   5.99e-01   4.26e+00
clark_level_at_diagnosis3   3.21e-01   1.65e+00   8.80e-01   3.09e+00
clark_level_at_diagnosis4   3.12e-01   1.31e+00   7.09e-01   2.41e+00
clark_level_at_diagnosis5   3.52e-01   7.43e-01   3.73e-01   1.48e+00
submitted_tumor_site.12     3.29e-01   3.61e-01   1.89e-01   6.87e-01
submitted_tumor_site.13     2.40e-01   1.40e+00   8.74e-01   2.24e+00
submitted_tumor_site.14     1.86e-01   9.53e-01   6.61e-01   1.37e+00

N = 355,  Events: 137,  Censored: 218
Total time at risk: 618038
Log-likelihood = -1213.212, df = 11
AIC = 2448.424
代码
```{r}
new_dat <- expand.grid(
  age_at_diagnosis = median(skcm$age_at_diagnosis),
  tumor_status = levels(skcm$tumor_status),
  clark_level_at_diagnosis = levels(skcm$clark_level_at_diagnosis), submitted_tumor_site.1 = levels(skcm$submitted_tumor_site.1)
)
pre_weibull <- weibull_fit |> predict(type = "quantile", p = 0.5, newdata = new_dat)
head(bind_cols(pre_weibull, new_dat))
```
# A tibble: 6 × 6
  .quantile .pred_quantile age_at_diagnosis tumor_status clark_level_at_diagno…¹
      <dbl>          <dbl>            <dbl> <fct>        <fct>                  
1       0.5         13668.               57 1            1                      
2       0.5          1520.               57 2            1                      
3       0.5         21825.               57 1            2                      
4       0.5          2428.               57 2            2                      
5       0.5         22547.               57 1            3                      
6       0.5          2508.               57 2            3                      
# ℹ abbreviated name: ¹​clark_level_at_diagnosis
# ℹ 1 more variable: submitted_tumor_site.1 <fct>

3.4 建立cox regression model

3.4.1 the proportional hazards assumption

h1(t)=ψh0(t)

代码
```{r}
cox_fit <- coxph(Surv(days, status) ~ ., data = skcm)
cox_fit
(cox_tbl <- cox_fit %>% gtsummary::tbl_regression(exponentiate = TRUE, ))
```
Call:
coxph(formula = Surv(days, status) ~ ., data = skcm)

                                        coef exp(coef)  se(coef)      z
breslow_thickness_at_diagnosis      0.011759  1.011828  0.011842  0.993
primary_melanoma_mitotic_rate       0.001221  1.001222  0.014029  0.087
age_at_diagnosis                    0.023979  1.024269  0.007103  3.376
prospective_collection2             0.122019  1.129775  0.274397  0.445
gender2                             0.082628  1.086137  0.197449  0.418
tumor_status2                       2.718725 15.160975  0.393111  6.916
clark_level_at_diagnosis2          -0.368196  0.691981  0.703580 -0.523
clark_level_at_diagnosis3          -0.453592  0.635342  0.505061 -0.898
clark_level_at_diagnosis4          -0.172248  0.841770  0.489461 -0.352
clark_level_at_diagnosis5           0.537417  1.711580  0.492252  1.092
primary_melanoma_tumor_ulceration2  0.186059  1.204494  0.236393  0.787
ajcc_tumor_pathologic_pt           -0.009946  0.990104  0.038027 -0.262
ajcc_nodes_pathologic_pn            0.004887  1.004899  0.029744  0.164
submitted_tumor_site.12             1.009340  2.743790  0.500491  2.017
submitted_tumor_site.13            -0.512568  0.598956  0.311159 -1.647
submitted_tumor_site.14             0.035049  1.035671  0.246194  0.142
radiation_treatment_adjuvant2       0.332710  1.394743  0.519499  0.640
new_tumor_event_dx_indicator2      -0.046874  0.954208  0.302545 -0.155
BMI                                -0.253083  0.776403  0.409466 -0.618
                                          p
breslow_thickness_at_diagnosis     0.320719
primary_melanoma_mitotic_rate      0.930623
age_at_diagnosis                   0.000735
prospective_collection2            0.656552
gender2                            0.675599
tumor_status2                      4.65e-12
clark_level_at_diagnosis2          0.600753
clark_level_at_diagnosis3          0.369135
clark_level_at_diagnosis4          0.724903
clark_level_at_diagnosis5          0.274942
primary_melanoma_tumor_ulceration2 0.431237
ajcc_tumor_pathologic_pt           0.793672
ajcc_nodes_pathologic_pn           0.869501
submitted_tumor_site.12            0.043727
submitted_tumor_site.13            0.099499
submitted_tumor_site.14            0.886791
radiation_treatment_adjuvant2      0.521883
new_tumor_event_dx_indicator2      0.876875
BMI                                0.536522

Likelihood ratio test=151  on 19 df, p=< 2.2e-16
n= 355, number of events= 137
Characteristic HR 95% CI p-value
breslow_thickness_at_diagnosis 1.01 0.99, 1.04 0.3
primary_melanoma_mitotic_rate 1.00 0.97, 1.03 >0.9
age_at_diagnosis 1.02 1.01, 1.04 <0.001
prospective_collection


    1 — —
    2 1.13 0.66, 1.93 0.7
gender


    1 — —
    2 1.09 0.74, 1.60 0.7
tumor_status


    1 — —
    2 15.2 7.02, 32.8 <0.001
clark_level_at_diagnosis


    1 — —
    2 0.69 0.17, 2.75 0.6
    3 0.64 0.24, 1.71 0.4
    4 0.84 0.32, 2.20 0.7
    5 1.71 0.65, 4.49 0.3
primary_melanoma_tumor_ulceration


    1 — —
    2 1.20 0.76, 1.91 0.4
ajcc_tumor_pathologic_pt 0.99 0.92, 1.07 0.8
ajcc_nodes_pathologic_pn 1.00 0.95, 1.07 0.9
submitted_tumor_site.1


    1 — —
    2 2.74 1.03, 7.32 0.044
    3 0.60 0.33, 1.10 0.10
    4 1.04 0.64, 1.68 0.9
radiation_treatment_adjuvant


    1 — —
    2 1.39 0.50, 3.86 0.5
new_tumor_event_dx_indicator


    1 — —
    2 0.95 0.53, 1.73 0.9
BMI 0.78 0.35, 1.73 0.5
Abbreviations: CI = Confidence Interval, HR = Hazard Ratio

故从age_at_diagnosis,tumor_status与submitted_tumor_site.1开始创建model

代码
```{r}
cox_fit1 <- coxph(Surv(days, status) ~ age_at_diagnosis, data = skcm)
cox_fit2 <- coxph(Surv(days, status) ~ tumor_status, data = skcm)
cox_fit3 <- coxph(Surv(days, status) ~ submitted_tumor_site.1, data = skcm)
cox_fit1
cox_fit2
cox_fit3
```
Call:
coxph(formula = Surv(days, status) ~ age_at_diagnosis, data = skcm)

                    coef exp(coef) se(coef)    z        p
age_at_diagnosis 0.02412   1.02441  0.00597 4.04 5.34e-05

Likelihood ratio test=17.15  on 1 df, p=3.46e-05
n= 355, number of events= 137 
Call:
coxph(formula = Surv(days, status) ~ tumor_status, data = skcm)

                 coef exp(coef) se(coef)     z        p
tumor_status2  2.6398   14.0098   0.3885 6.794 1.09e-11

Likelihood ratio test=101.4  on 1 df, p=< 2.2e-16
n= 355, number of events= 137 
Call:
coxph(formula = Surv(days, status) ~ submitted_tumor_site.1, 
    data = skcm)

                           coef exp(coef) se(coef)      z       p
submitted_tumor_site.12  1.1186    3.0605   0.3972  2.816 0.00487
submitted_tumor_site.13 -0.2647    0.7675   0.2767 -0.956 0.33885
submitted_tumor_site.14 -0.1226    0.8846   0.2186 -0.561 0.57489

Likelihood ratio test=9.85  on 3 df, p=0.01987
n= 355, number of events= 137

由p值大小先将tumor_status加入model

逐步AIC确定最终model:

代码
```{r}
result <- step(cox_fit, scope = list(upper = ~., lower = ~tumor_status))
```
Start:  AIC=1218.24
Surv(days, status) ~ breslow_thickness_at_diagnosis + primary_melanoma_mitotic_rate + 
    age_at_diagnosis + prospective_collection + gender + tumor_status + 
    clark_level_at_diagnosis + primary_melanoma_tumor_ulceration + 
    ajcc_tumor_pathologic_pt + ajcc_nodes_pathologic_pn + submitted_tumor_site.1 + 
    radiation_treatment_adjuvant + new_tumor_event_dx_indicator + 
    BMI

                                    Df    AIC
- primary_melanoma_mitotic_rate      1 1216.2
- new_tumor_event_dx_indicator       1 1216.3
- ajcc_nodes_pathologic_pn           1 1216.3
- ajcc_tumor_pathologic_pt           1 1216.3
- gender                             1 1216.4
- prospective_collection             1 1216.4
- BMI                                1 1216.6
- radiation_treatment_adjuvant       1 1216.6
- primary_melanoma_tumor_ulceration  1 1216.9
- breslow_thickness_at_diagnosis     1 1217.1
<none>                                 1218.2
- clark_level_at_diagnosis           4 1218.3
- submitted_tumor_site.1             3 1222.7
- age_at_diagnosis                   1 1228.1

Step:  AIC=1216.25
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    prospective_collection + gender + tumor_status + clark_level_at_diagnosis + 
    primary_melanoma_tumor_ulceration + ajcc_tumor_pathologic_pt + 
    ajcc_nodes_pathologic_pn + submitted_tumor_site.1 + radiation_treatment_adjuvant + 
    new_tumor_event_dx_indicator + BMI

                                    Df    AIC
- new_tumor_event_dx_indicator       1 1214.3
- ajcc_nodes_pathologic_pn           1 1214.3
- ajcc_tumor_pathologic_pt           1 1214.3
- gender                             1 1214.4
- prospective_collection             1 1214.4
- BMI                                1 1214.6
- radiation_treatment_adjuvant       1 1214.7
- primary_melanoma_tumor_ulceration  1 1214.9
- breslow_thickness_at_diagnosis     1 1215.2
<none>                                 1216.2
- clark_level_at_diagnosis           4 1216.6
+ primary_melanoma_mitotic_rate      1 1218.2
- submitted_tumor_site.1             3 1221.0
- age_at_diagnosis                   1 1226.1

Step:  AIC=1214.28
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    prospective_collection + gender + tumor_status + clark_level_at_diagnosis + 
    primary_melanoma_tumor_ulceration + ajcc_tumor_pathologic_pt + 
    ajcc_nodes_pathologic_pn + submitted_tumor_site.1 + radiation_treatment_adjuvant + 
    BMI

                                    Df    AIC
- ajcc_nodes_pathologic_pn           1 1212.3
- ajcc_tumor_pathologic_pt           1 1212.3
- gender                             1 1212.5
- prospective_collection             1 1212.5
- BMI                                1 1212.7
- radiation_treatment_adjuvant       1 1212.7
- primary_melanoma_tumor_ulceration  1 1212.9
- breslow_thickness_at_diagnosis     1 1213.3
<none>                                 1214.3
- clark_level_at_diagnosis           4 1214.6
+ new_tumor_event_dx_indicator       1 1216.2
+ primary_melanoma_mitotic_rate      1 1216.3
- submitted_tumor_site.1             3 1219.0
- age_at_diagnosis                   1 1224.3

Step:  AIC=1212.32
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    prospective_collection + gender + tumor_status + clark_level_at_diagnosis + 
    primary_melanoma_tumor_ulceration + ajcc_tumor_pathologic_pt + 
    submitted_tumor_site.1 + radiation_treatment_adjuvant + BMI

                                    Df    AIC
- ajcc_tumor_pathologic_pt           1 1210.4
- prospective_collection             1 1210.6
- gender                             1 1210.6
- BMI                                1 1210.7
- radiation_treatment_adjuvant       1 1210.8
- primary_melanoma_tumor_ulceration  1 1211.0
- breslow_thickness_at_diagnosis     1 1211.4
<none>                                 1212.3
- clark_level_at_diagnosis           4 1212.6
+ ajcc_nodes_pathologic_pn           1 1214.3
+ new_tumor_event_dx_indicator       1 1214.3
+ primary_melanoma_mitotic_rate      1 1214.3
- submitted_tumor_site.1             3 1217.1
- age_at_diagnosis                   1 1222.4

Step:  AIC=1210.39
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    prospective_collection + gender + tumor_status + clark_level_at_diagnosis + 
    primary_melanoma_tumor_ulceration + submitted_tumor_site.1 + 
    radiation_treatment_adjuvant + BMI

                                    Df    AIC
- prospective_collection             1 1208.6
- gender                             1 1208.6
- BMI                                1 1208.8
- radiation_treatment_adjuvant       1 1208.8
- primary_melanoma_tumor_ulceration  1 1209.0
- breslow_thickness_at_diagnosis     1 1209.4
<none>                                 1210.4
- clark_level_at_diagnosis           4 1212.3
+ ajcc_tumor_pathologic_pt           1 1212.3
+ new_tumor_event_dx_indicator       1 1212.3
+ ajcc_nodes_pathologic_pn           1 1212.3
+ primary_melanoma_mitotic_rate      1 1212.4
- submitted_tumor_site.1             3 1215.2
- age_at_diagnosis                   1 1220.5

Step:  AIC=1208.61
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    gender + tumor_status + clark_level_at_diagnosis + primary_melanoma_tumor_ulceration + 
    submitted_tumor_site.1 + radiation_treatment_adjuvant + BMI

                                    Df    AIC
- gender                             1 1206.9
- BMI                                1 1206.9
- radiation_treatment_adjuvant       1 1207.1
- primary_melanoma_tumor_ulceration  1 1207.3
- breslow_thickness_at_diagnosis     1 1207.8
<none>                                 1208.6
- clark_level_at_diagnosis           4 1210.3
+ prospective_collection             1 1210.4
+ new_tumor_event_dx_indicator       1 1210.5
+ ajcc_nodes_pathologic_pn           1 1210.5
+ ajcc_tumor_pathologic_pt           1 1210.6
+ primary_melanoma_mitotic_rate      1 1210.6
- submitted_tumor_site.1             3 1213.7
- age_at_diagnosis                   1 1218.6

Step:  AIC=1206.89
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    tumor_status + clark_level_at_diagnosis + primary_melanoma_tumor_ulceration + 
    submitted_tumor_site.1 + radiation_treatment_adjuvant + BMI

                                    Df    AIC
- BMI                                1 1205.2
- radiation_treatment_adjuvant       1 1205.3
- primary_melanoma_tumor_ulceration  1 1205.5
- breslow_thickness_at_diagnosis     1 1206.1
<none>                                 1206.9
- clark_level_at_diagnosis           4 1208.3
+ gender                             1 1208.6
+ prospective_collection             1 1208.6
+ new_tumor_event_dx_indicator       1 1208.7
+ ajcc_nodes_pathologic_pn           1 1208.8
+ ajcc_tumor_pathologic_pt           1 1208.8
+ primary_melanoma_mitotic_rate      1 1208.9
- submitted_tumor_site.1             3 1212.0
- age_at_diagnosis                   1 1216.9

Step:  AIC=1205.22
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    tumor_status + clark_level_at_diagnosis + primary_melanoma_tumor_ulceration + 
    submitted_tumor_site.1 + radiation_treatment_adjuvant

                                    Df    AIC
- radiation_treatment_adjuvant       1 1203.8
- primary_melanoma_tumor_ulceration  1 1203.9
- breslow_thickness_at_diagnosis     1 1204.3
<none>                                 1205.2
- clark_level_at_diagnosis           4 1206.4
+ BMI                                1 1206.9
+ gender                             1 1206.9
+ prospective_collection             1 1207.1
+ new_tumor_event_dx_indicator       1 1207.1
+ ajcc_nodes_pathologic_pn           1 1207.1
+ ajcc_tumor_pathologic_pt           1 1207.2
+ primary_melanoma_mitotic_rate      1 1207.2
- submitted_tumor_site.1             3 1210.5
- age_at_diagnosis                   1 1214.9

Step:  AIC=1203.8
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    tumor_status + clark_level_at_diagnosis + primary_melanoma_tumor_ulceration + 
    submitted_tumor_site.1

                                    Df    AIC
- primary_melanoma_tumor_ulceration  1 1202.3
- breslow_thickness_at_diagnosis     1 1203.0
<none>                                 1203.8
- clark_level_at_diagnosis           4 1205.2
+ radiation_treatment_adjuvant       1 1205.2
+ BMI                                1 1205.3
+ ajcc_nodes_pathologic_pn           1 1205.5
+ gender                             1 1205.5
+ prospective_collection             1 1205.6
+ new_tumor_event_dx_indicator       1 1205.7
+ primary_melanoma_mitotic_rate      1 1205.8
+ ajcc_tumor_pathologic_pt           1 1205.8
- submitted_tumor_site.1             3 1208.8
- age_at_diagnosis                   1 1213.2

Step:  AIC=1202.27
Surv(days, status) ~ breslow_thickness_at_diagnosis + age_at_diagnosis + 
    tumor_status + clark_level_at_diagnosis + submitted_tumor_site.1

                                    Df    AIC
- breslow_thickness_at_diagnosis     1 1201.6
<none>                                 1202.3
+ BMI                                1 1203.7
+ primary_melanoma_tumor_ulceration  1 1203.8
+ ajcc_nodes_pathologic_pn           1 1203.9
+ radiation_treatment_adjuvant       1 1203.9
+ prospective_collection             1 1204.1
+ gender                             1 1204.1
+ new_tumor_event_dx_indicator       1 1204.1
+ primary_melanoma_mitotic_rate      1 1204.2
+ ajcc_tumor_pathologic_pt           1 1204.2
- clark_level_at_diagnosis           4 1204.6
- submitted_tumor_site.1             3 1207.3
- age_at_diagnosis                   1 1213.2

Step:  AIC=1201.59
Surv(days, status) ~ age_at_diagnosis + tumor_status + clark_level_at_diagnosis + 
    submitted_tumor_site.1

                                    Df    AIC
<none>                                 1201.6
+ breslow_thickness_at_diagnosis     1 1202.3
+ ajcc_nodes_pathologic_pn           1 1203.0
+ primary_melanoma_tumor_ulceration  1 1203.0
+ radiation_treatment_adjuvant       1 1203.2
+ prospective_collection             1 1203.2
+ new_tumor_event_dx_indicator       1 1203.2
+ BMI                                1 1203.2
+ primary_melanoma_mitotic_rate      1 1203.3
+ gender                             1 1203.4
+ ajcc_tumor_pathologic_pt           1 1203.4
- clark_level_at_diagnosis           4 1205.1
- submitted_tumor_site.1             3 1211.3
- age_at_diagnosis                   1 1212.1
代码
```{r}
result
```
Call:
coxph(formula = Surv(days, status) ~ age_at_diagnosis + tumor_status + 
    clark_level_at_diagnosis + submitted_tumor_site.1, data = skcm)

                               coef exp(coef)  se(coef)      z        p
age_at_diagnosis           0.023002  1.023268  0.006639  3.465 0.000531
tumor_status2              2.731072 15.349333  0.391476  6.976 3.03e-12
clark_level_at_diagnosis2 -0.563113  0.569434  0.621023 -0.907 0.364539
clark_level_at_diagnosis3 -0.588794  0.554996  0.403541 -1.459 0.144546
clark_level_at_diagnosis4 -0.300938  0.740124  0.389901 -0.772 0.440214
clark_level_at_diagnosis5  0.363886  1.438910  0.440077  0.827 0.408312
submitted_tumor_site.12    1.375597  3.957438  0.421945  3.260 0.001114
submitted_tumor_site.13   -0.431385  0.649609  0.300168 -1.437 0.150677
submitted_tumor_site.14    0.071662  1.074292  0.233326  0.307 0.758743

Likelihood ratio test=147.7  on 9 df, p=< 2.2e-16
n= 355, number of events= 137
代码
```{r}
termplot(result, se = T, terms = 1, ylabs = "Log hazard")
```

代码
```{r}
cox_fit <- result
# Predictions will be for all levels of sex and ph.ecog, but only at the median
# age.
new_dat <- expand.grid(
  age_at_diagnosis = median(skcm$age_at_diagnosis),
  tumor_status = levels(skcm$tumor_status),
  clark_level_at_diagnosis = levels(skcm$clark_level_at_diagnosis), submitted_tumor_site.1 = levels(skcm$submitted_tumor_site.1)
) %>%
  # strata is our key to join back to the fitted values.
  mutate(strata = as.factor(row_number()))

# Create fitted survival curves at the covariate presets.
fit_curves <- survfit(cox_fit, newdata = new_dat, data = skcm)

# `surv_summary()` is like `summary()` except that it includes risk table info,
# confidence interval attributes, and pivots the strata longer.
surv_summary <- surv_summary(fit_curves) %>%
  # The cases are labeled "strata", but `survsummary()` doesn't label what the
  # strata are! Get it from new_dat.
  inner_join(new_dat, by = "strata")

# Now use ggplot() just like normal.
median_line <- surv_summary %>%
  filter(surv >= .5) %>%
  group_by(tumor_status, clark_level_at_diagnosis, submitted_tumor_site.1) %>%
  summarize(.groups = "drop", max_t = max(time))
surv_summary %>%
  ggplot(aes(x = time, y = surv)) +
  geom_line(aes(color = clark_level_at_diagnosis)) +
  geom_ribbon(aes(ymin = lower, ymax = upper, color = clark_level_at_diagnosis, fill = clark_level_at_diagnosis),
    alpha = 0.4
  ) +
  geom_segment(
    data = median_line,
    aes(x = 0, xend = max_t, y = .5, yend = .5),
    linetype = 2, color = "gray40"
  ) +
  geom_segment(
    data = median_line,
    aes(x = max_t, xend = max_t, y = 0, yend = .5, color = clark_level_at_diagnosis),
    linetype = 2
  ) +
  facet_wrap(facets = vars(tumor_status, submitted_tumor_site.1)) +
  scale_y_continuous(labels = percent_format(1)) +
  theme_light() +
  theme(legend.position = "top") +
  labs(
    X = "days", y = "Survival Probability", color = NULL, fill = NULL,
    title = "Cox fitted model",
    subtitle = "Age held at median."
  )
```

gender 分1,2 submitted_tumor_site.1分1,2,3,4

4 模型验证

4.1 Assessing Goodness of Fit Using Residuals

4.1.1 Martingale and Deviance Residuals

代码
```{r}
smoothSEcurve <- function(yy, xx) {
  # use after a call to "plot"
  # fit a lowess curve and 95% confidence interval curve
  # make list of x values
  xx.list <- min(xx) + ((0:100) / 100) * (max(xx) - min(xx))
  # Then fit loess function through the points (xx, yy)
  # at the listed values
  yy.xx <- predict(loess(yy ~ xx),
    se = T,
    newdata = data.frame(xx = xx.list)
  )
  lines(yy.xx$fit ~ xx.list, lwd = 2)
  lines(yy.xx$fit -
    qt(0.975, yy.xx$df) * yy.xx$se.fit ~ xx.list, lty = 2)
  lines(yy.xx$fit +
    qt(0.975, yy.xx$df) * yy.xx$se.fit ~ xx.list, lty = 2)
}
skcm1 <- select(skcm, age_at_diagnosis, tumor_status, clark_level_at_diagnosis, submitted_tumor_site.1)
rr.final <- residuals(cox_fit, type = "martingale")
par(mfrow = c(1, 2))
plot(rr.final ~ skcm1$age_at_diagnosis)
smoothSEcurve(rr.final, skcm1$age_at_diagnosis)
title("Martingale residuals\nversus age")
plot(rr.final ~ skcm1$tumor_status)
title("Martingale residuals\nversus tumor_status")
```

代码
```{r}
par(mfrow = c(1, 2))
plot(rr.final ~ skcm1$clark_level_at_diagnosis)
title("Martingale residuals\nversus clark_level_at_diagnosis")
plot(rr.final ~ skcm1$submitted_tumor_site.1)
title("Martingale residuals\nversus submitted_tumor_site.1")
```

残差大致均匀分布在0周围,可见model拟合效果不错

4.1.2 Case Deletion Residuals

代码
```{r}
ggcoxdiagnostics(cox_fit, type = "dfbeta", linear.predictions = FALSE)
```

4.2 Checking the Proportion Hazards Assumption

4.2.1 Log Cumulative Hazard Plots

代码
```{r}
fit <- survfit(Surv(days, status) ~ tumor_status, data = skcm)
fit$label <-
  c(rep(1, fit$strata[1]), rep(2, fit$strata[2]))
t1 <- c(0, subset(fit$time, fit$label == 1))
t2 <- c(0, subset(fit$time, fit$label == 2))

St1 <- c(1, subset(fit$surv, fit$label == 1))
St2 <- c(1, subset(fit$surv, fit$label == 2))
plot(t1, -log(-log(St1)), col = 2, type = "s", xlab = "Time", ylab = expression(-ln(-
  ln(hat(S)))), main = "Log-Log KM Curves for each tumor_status")

lines(t2, -log(-log(St2)), col = 3, type = "s")
legend("topright", c("tumor_status=1", "tumor_status=2"), col = 2:3, lty = 1, cex = 0.8)
```

4.2.2 Schoenfeld Residuals

代码
```{r}
(cox_test_ph <- cox.zph(cox_fit))
ggcoxzph(cox_test_ph)
```
                         chisq df    p
age_at_diagnosis         0.314  1 0.58
tumor_status             0.930  1 0.33
clark_level_at_diagnosis 0.338  4 0.99
submitted_tumor_site.1   1.490  3 0.68
GLOBAL                   2.782  9 0.97

5 结果可视化

5.1 ROC曲线

代码
```{r}
train_model <- cox_fit
train <- skcm
library(timeROC)

train$lp <- predict(train_model, newdata = skcm1, type = "lp")
time_roc_train <- timeROC(
  T = train$days, # 结局时间
  delta = train$status, # 生存结局
  marker = train$lp, # 预测变量
  cause = 1, # 阳性结局指标值
  weighting = "marginal", # 权重计算方法，默认marginal，采用KM计算删失分布
  times = c(365, 730, 1100),
  ROC = TRUE, # 保存sensitivities 和 specificties值
  iid = TRUE # iid = TRUE 才会保存置信区间，但是样本量大了后，耗时耗资源
)

# 计算train-1年的AUC
train1y <- paste0(
  "train 1 year AUC (95%CI) = ",
  sprintf("%.3f", time_roc_train$AUC[1]), " (",
  sprintf("%.3f", confint(time_roc_train, level = 0.95)$CI_AUC[1, 1] / 100), " - ",
  sprintf("%.3f", confint(time_roc_train, level = 0.95)$CI_AUC[1, 2] / 100), ")"
)

# 计算train-3年的AUC
train3y <- paste0(
  "train 3 year AUC (95%CI) = ",
  sprintf("%.3f", time_roc_train$AUC[2]), " (",
  sprintf("%.3f", confint(time_roc_train, level = 0.95)$CI_AUC[2, 1] / 100), " - ",
  sprintf("%.3f", confint(time_roc_train, level = 0.95)$CI_AUC[2, 2] / 100), ")"
)

# 计算train-5年的AUC
train5y <- paste0(
  "train 5 year AUC (95%CI) = ",
  sprintf("%.3f", time_roc_train$AUC[3]), " (",
  sprintf("%.3f", confint(time_roc_train, level = 0.95)$CI_AUC[3, 1] / 100), " - ",
  sprintf("%.3f", confint(time_roc_train, level = 0.95)$CI_AUC[3, 2] / 100), ")"
)

# train AUC 多时间点合一
plot(title = "", time_roc_train, col = "DodgerBlue", time = 365, lty = 1, lwd = 2, family = "serif") # 绘制ROC曲线
plot(time_roc_train, time = 730, lty = 1, lwd = 2, add = TRUE, col = "LightSeaGreen", family = "serif") # add=TRUE指在前一条曲线基础上新增
plot(time_roc_train, time = 1100, lty = 1, lwd = 2, add = TRUE, col = "DarkOrange", family = "serif")
legend("bottomright", c(train1y, train3y, train5y), col = c("DodgerBlue", "LightSeaGreen", "DarkOrange"), lty = 1, lwd = 2)
```

效果不是很好😅

5.2 nomogram

代码
```{r}
skcm2 <- skcm |> mutate(months = floor(days / 30))

d <- datadist(skcm2)
options(datadist = "d")
f <- cph(formula = Surv(months, status) ~ age_at_diagnosis + tumor_status +
  clark_level_at_diagnosis + submitted_tumor_site.1, data = skcm2, surv = T)
med <- Quantile(f)
surv <- Survival(f) # This would also work if f was from cph
plot(nomogram(f, fun = function(x) med(lp = x), funlabel = "Median Survival Time"))
nom <- nomogram(f,
  fun = list(
    function(x) surv(6, x),
    function(x) surv(12, x),
    function(x) surv(24, x),
    function(x) surv(36, x),
    function(x) surv(72, x)
  ),
  funlabel = c(
    "6-Month Survival Probability",
    "12-month Survival Probability",
    "24-month Survival Probability",
    "36-month Survival Probability",
    "72-month Survival Probability"
  )
)
plot(nom, xfrac = .7)
```

5.3 calibration plot

代码
```{r}
#| fig-cap: Calibration curve for 2 year survival. X-axis shows the nomogram predicted probability, while the Y-axis is actual 2-year survival as estimated by the Kaplan-Meier method.

d <- datadist(skcm)
options(datadist = "d")
f <- cph(formula = Surv(days, status) ~ age_at_diagnosis + tumor_status +
  clark_level_at_diagnosis + submitted_tumor_site.1, data = skcm, x = TRUE, y = TRUE, surv = TRUE, time.inc = 365)
pa <- requireNamespace("polspline")
if (pa) {
  cal <- calibrate(f, u = 365, B = 20) # cmethod='hare'
  plot(cal)
}
```
Using Cox survival estimates at   365 Days

Calibration curve for 2 year survival. X-axis shows the nomogram predicted probability, while the Y-axis is actual 2-year survival as estimated by the Kaplan-Meier method.

6 总结

代码
```{r}
#| label: tbl1
#| tbl-cap: summary model

(cox_tbl <- cox_fit %>% gtsummary::tbl_regression(exponentiate = TRUE, ))
```
Characteristic HR 95% CI p-value
age_at_diagnosis 1.02 1.01, 1.04 <0.001
tumor_status


    1 — —
    2 15.3 7.13, 33.1 <0.001
clark_level_at_diagnosis


    1 — —
    2 0.57 0.17, 1.92 0.4
    3 0.55 0.25, 1.22 0.14
    4 0.74 0.34, 1.59 0.4
    5 1.44 0.61, 3.41 0.4
submitted_tumor_site.1


    1 — —
    2 3.96 1.73, 9.05 0.001
    3 0.65 0.36, 1.17 0.2
    4 1.07 0.68, 1.70 0.8
Abbreviations: CI = Confidence Interval, HR = Hazard Ratio
summary model

如表所示,根据model,年龄与肿瘤情况诊断为2,clark_level_at_diagnosis 诊断为5,肿瘤位置为2会增大风险

附录

代码
```{r}
sessionInfo()
```
R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=Chinese (Simplified)_China.utf8 
[2] LC_CTYPE=Chinese (Simplified)_China.utf8   
[3] LC_MONETARY=Chinese (Simplified)_China.utf8
[4] LC_NUMERIC=C                               
[5] LC_TIME=Chinese (Simplified)_China.utf8    

time zone: Asia/Shanghai
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] timeROC_0.4        flexsurv_2.3.2     muhaz_1.2.6.4      tidycmprsk_1.1.0  
 [5] ggsurvfit_1.1.0    rsample_1.2.1      mice_3.17.0        scales_1.3.0      
 [9] conflicted_1.2.0   rms_7.0-0          Hmisc_5.2-2        gtsummary_2.1.0   
[13] lubridate_1.9.4    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
[17] purrr_1.0.2        tidyr_1.3.1        tibble_3.2.1       tidyverse_2.0.0   
[21] see_0.10.0         report_0.6.1       parameters_0.24.1  performance_0.13.0
[25] modelbased_0.9.0   insight_1.0.2      effectsize_1.0.0   datawizard_1.0.0  
[29] correlation_0.8.6  bayestestR_0.15.2  easystats_0.7.4    survMisc_0.5.6    
[33] readr_2.1.5        autoReg_0.3.3      survminer_0.5.0    ggpubr_0.6.0      
[37] survival_3.8-3     ellmer_0.1.1       showtext_0.9-7     showtextdb_3.0    
[41] sysfonts_0.8.9     ggthemr_1.1.0      ggplot2_3.5.1     

loaded via a namespace (and not attached):
  [1] sf_1.0-19               numDeriv_2016.8-1.1     tools_4.4.2            
  [4] backports_1.5.0         sjlabelled_1.2.0        utf8_1.2.4             
  [7] R6_2.6.1                lazyeval_0.2.2          mgcv_1.9-1             
 [10] jomo_2.7-6              withr_3.0.2             gridExtra_2.3          
 [13] quantreg_6.00           cli_3.6.2               textshaping_1.0.0      
 [16] gt_0.11.1               exactRankTests_0.8-35   officer_0.6.7          
 [19] sandwich_3.1-1          sass_0.4.9              labeling_0.4.3         
 [22] mvtnorm_1.3-3           S7_0.2.0                polspline_1.1.25       
 [25] proxy_0.4-27            askpass_1.2.1           commonmark_1.9.2       
 [28] systemfonts_1.2.1       foreign_0.8-87          parallelly_1.42.0      
 [31] labelled_2.14.0         rstudioapi_0.17.1       generics_0.1.3         
 [34] shape_1.4.6.1           vroom_1.6.5             car_3.1-3              
 [37] zip_2.3.2               Matrix_1.7-1            abind_1.4-8            
 [40] lifecycle_1.0.4         multcomp_1.4-28         yaml_2.3.10            
 [43] carData_3.0-5           mstate_0.3.3            grid_4.4.2             
 [46] crayon_1.5.3            mitml_0.4-5             lattice_0.22-6         
 [49] haven_2.5.4             pillar_1.10.1           knitr_1.49             
 [52] epiR_2.0.80             boot_1.3-31             estimability_1.5.1     
 [55] future.apply_1.11.3     codetools_0.2-20        pan_1.9                
 [58] glue_1.8.0              fontLiberation_0.1.0    data.table_1.16.4      
 [61] broom.helpers_1.19.0    vctrs_0.6.5             Rdpack_2.6.2           
 [64] gtable_0.3.6            assertthat_0.2.1        cachem_1.1.0           
 [67] xfun_0.50               rbibutils_2.3           prodlim_2024.06.25     
 [70] coda_0.19-4.1           reformulas_0.4.0        iterators_1.0.14       
 [73] KMsurv_0.1-5            units_0.8-5             lava_1.8.1             
 [76] statmod_1.5.0           TH.data_1.1-3           nlme_3.1-166           
 [79] bit64_4.6.0-1           fontquiver_0.2.1        maxstat_0.7-25         
 [82] KernSmooth_2.23-24      rpart_4.1.23            colorspace_2.1-0       
 [85] DBI_1.2.3               nnet_7.3-19             tidyselect_1.2.1       
 [88] emmeans_1.10.7          bit_4.5.0.1             compiler_4.4.2         
 [91] glmnet_4.1-8            httr2_1.1.0             htmlTable_2.4.3        
 [94] BiasedUrn_2.0.12        flextable_0.9.7         SparseM_1.84-2         
 [97] pammtools_0.5.93        xml2_1.3.6              fontBitstreamVera_0.1.1
[100] checkmate_2.3.2         classInt_0.4-11         pec_2023.04.12         
[103] quadprog_1.5-8          rappdirs_0.3.3          digest_0.6.35          
[106] minqa_1.2.8             rmarkdown_2.29          htmltools_0.5.8.1      
[109] pkgconfig_2.0.3         base64enc_0.1-3         coro_1.1.0             
[112] lme4_1.1-36             fastmap_1.2.0           rlang_1.1.4            
[115] htmlwidgets_1.6.4       farver_2.1.2            zoo_1.8-13             
[118] jsonlite_1.9.0          magrittr_2.0.3          Formula_1.2-5          
[121] patchwork_1.3.0         munsell_0.5.1           Rcpp_1.0.13-1          
[124] gdtools_0.4.1           furrr_0.3.1             stringi_1.8.4          
[127] MASS_7.3-61             parallel_4.4.2          listenv_0.9.1          
[130] sjmisc_2.8.10           splines_4.4.2           pander_0.6.5           
[133] hms_1.1.3               timereg_2.0.6           uuid_1.2-1             
[136] markdown_1.13           ggsignif_0.6.4          evaluate_1.0.3         
[139] deSolve_1.40            nloptr_2.1.1            tzdb_0.4.0             
[142] foreach_1.5.2           MatrixModels_0.5-3      moonBook_0.3.1         
[145] cards_0.5.0             openssl_2.3.2           future_1.34.0          
[148] km.ci_0.5-6             broom_1.0.7             xtable_1.8-4           
[151] e1071_1.7-16            rstatix_0.7.2           class_7.3-22           
[154] ragg_1.3.3              memoise_2.0.1           cluster_2.1.6          
[157] timechange_0.3.0        globals_0.16.3