The survival analysis reference used was writen by Zabor
Transfer libraries
The survival analysis reference used was Jaeger.
To add p-values to plots, I used the Kassambara 2017 and 2020 references
library(maditr)
## Warning: package 'maditr' was built under R version 4.1.2
##
## To modify variables or add new variables:
## let(mtcars, new_var = 42, new_var2 = new_var*hp) %>% head()
Load the GSVA object back to the R environment. This object now contains all the metadata information.
r2_gse62564_GSVA_Metadata <- readRDS( "../results/r2_gse62564_GSVA_Metadata.rds")
In neuroblastoma, the main risk classification system classifies patients using features such as age, stage, histology and MYCN amplification status.
MYCN is a gene that encodes a transcription factor that is very important in tumor prognosis, the association between tumor features and patient survival.
Let’s take a look at the MYCN expression compared by the neuroblastoma risk group.
First, let’s see which type of variable is the the MYCN status and the MYCN expression
class(r2_gse62564_GSVA_Metadata$mycn_status)
## [1] "character"
class(r2_gse62564_GSVA_Metadata$MYCN)
## [1] "character"
They are both of type “character.” In order to plot, we have to make the MYCN expression to be of type numeric
r2_gse62564_GSVA_Metadata$MYCN <- as.numeric(r2_gse62564_GSVA_Metadata$MYCN)
class(r2_gse62564_GSVA_Metadata$MYCN)
## [1] "numeric"
r2_gse62564_GSVA_Metadata$age_at_diagnosis <- as.numeric(r2_gse62564_GSVA_Metadata$age_at_diagnosis)
age_diagnosis_days_gse62564_498 <- mean(r2_gse62564_GSVA_Metadata$age_at_diagnosis)
Let’s try and see the boxplot of the MYCN expression vs. MYCN amplification status.
Using the boxplot function that we tried before, we will have an error.
boxplot(data=r2_gse62564_GSVA_Metadata, r2_gse62564_GSVA_Metadata$mycn_status, r2_gse62564_GSVA_Metadata$MYCN)
Maybe we could try the boxplot function making the mycn_status a numeric variable.
boxplot(data=r2_gse62564_GSVA_Metadata, as.numeric(r2_gse62564_GSVA_Metadata$mycn_status), r2_gse62564_GSVA_Metadata$MYCN)
Let’s try and see if ggplot, which we have also used before, takes the mycn_status variable as categorical as it is right now.
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.2
ggplot(r2_gse62564_GSVA_Metadata, aes(x=mycn_status, y=MYCN, fill=mycn_status)) +
geom_boxplot()
In this figure, we see that group with clinical MYCN amplification, also shows high MYCN expression
What is the difference, in measurement, between MYCN amplification and MYCN expression?
How are the MYCN amplification and MYCN gene expression correlated?
For now, I am not bothering neither with formating of the X and Y labels nor with the N/A values of the MYCN amplification status.
We can plot the high risk status vs. MYCN amplification and leave the fill argument as mycn_status
library(ggplot2)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=MYCN, fill=mycn_status)) +
geom_boxplot()
Or we can make it the high_risk status instead
library(ggplot2)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=MYCN, fill=high_risk)) +
geom_boxplot()
Here we see that the high_risk status correlates positively with MYCN expression.
Make the OS (overall survival) status numeric, instead of character. Do the same for the EFS.
## OS
r2_gse62564_GSVA_Metadata <- r2_gse62564_GSVA_Metadata %>%
dplyr::mutate(os_bin_numeric =
ifelse(r2_gse62564_GSVA_Metadata$os_bin == "event", "1", "0"))
r2_gse62564_GSVA_Metadata["os_bin_numeric"] <- as.numeric(r2_gse62564_GSVA_Metadata$os_bin_numeric)
## EFS
r2_gse62564_GSVA_Metadata <- r2_gse62564_GSVA_Metadata %>%
dplyr::mutate(efs_bin_numeric =
ifelse(r2_gse62564_GSVA_Metadata$efs_bin == "event", "1", "0"))
r2_gse62564_GSVA_Metadata["efs_bin_numeric"] <- as.numeric(r2_gse62564_GSVA_Metadata$efs_bin_numeric)
library(survival)
## Warning: package 'survival' was built under R version 4.1.2
fit_os_high_risk <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata$high_risk,
data = r2_gse62564_GSVA_Metadata)
library(survminer)
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 4.1.2
##
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
##
## myeloma
ggsurvplot(fit_os_high_risk,
data = r2_gse62564_GSVA_Metadata,
risk.table = TRUE,
# risk.table = "abs_pct",
pval = TRUE,
pval.coord = c(2000, 1),
)
ggsurvplot(fit_os_high_risk,
data = r2_gse62564_GSVA_Metadata,
# risk.table = TRUE,
risk.table = "abs_pct",
pval = TRUE,
pval.coord = c(2000, 1),
)
Let’s now establish a relationship between neuroblastoma high risk status and hypoxia, one of the novelties of our work.
r2_gse62564_GSVA_Metadata$HIF1A <- as.numeric(r2_gse62564_GSVA_Metadata$HIF1A)
compare_means(data = r2_gse62564_GSVA_Metadata, HIF1A ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HIF1A, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test", label.x = 1.4)+
theme_linedraw()
We also must make HALLMARK_HYPOXIA, numeric
r2_gse62564_GSVA_Metadata$HALLMARK_HYPOXIA <- as.numeric(r2_gse62564_GSVA_Metadata$HALLMARK_HYPOXIA)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA, fill=high_risk)) +
geom_boxplot()
In the next chunk, we construct a data-frame with statistics calculations for our plot.
The reference used here was the Kassambara 2023 article:
library(rstatix)
## Warning: package 'rstatix' was built under R version 4.1.2
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
# Statistical test
stat.test <- r2_gse62564_GSVA_Metadata %>%
t_test(HALLMARK_HYPOXIA ~ high_risk) %>%
add_significance()
stat.test
In the next chunk, we create the object that will show the statistical results. After, we use the stat_pvalue_manual function.
We also change the theme of the ggplot2 figure plot.
library(ggpubr)
r2_gse62564_GSVA_Metadata$high_risk <- as.factor(r2_gse62564_GSVA_Metadata$high_risk)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA)) +
geom_boxplot()+
stat_pvalue_manual(stat.test, label = "p", y.position = 0.6) + # y.position is mandatory
scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))
r2_gse62564_GSVA_Metadata$high_risk <- as.factor(r2_gse62564_GSVA_Metadata$high_risk)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA, fill="high_risk")) +
geom_boxplot()+
stat_pvalue_manual(stat.test, label = "p", y.position = 0.6) + # y.position is mandatory
scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))
compare_means(data = r2_gse62564_GSVA_Metadata, HALLMARK_HYPOXIA ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_HYPOXIA, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
r2_gse62564_GSVA_Metadata$ELVIDGE_HIF1A_TARGETS_UP <- as.numeric(r2_gse62564_GSVA_Metadata$ELVIDGE_HIF1A_TARGETS_UP)
compare_means(data = r2_gse62564_GSVA_Metadata, ELVIDGE_HIF1A_TARGETS_UP ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ELVIDGE_HIF1A_TARGETS_UP, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_UP <- as.numeric(r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_UP)
compare_means(data = r2_gse62564_GSVA_Metadata, GROSS_HYPOXIA_VIA_HIF1A_UP ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=GROSS_HYPOXIA_VIA_HIF1A_UP, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_DN <- as.numeric(r2_gse62564_GSVA_Metadata$GROSS_HYPOXIA_VIA_HIF1A_DN)
compare_means(data = r2_gse62564_GSVA_Metadata, GROSS_HYPOXIA_VIA_HIF1A_DN ~ high_risk, method = "t.test")
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=GROSS_HYPOXIA_VIA_HIF1A_DN, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt <- as.numeric(r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ADRN_Norm_vs_Hypo_Down_635.txt, fill=high_risk)) +
geom_boxplot()
compare_means(data = r2_gse62564_GSVA_Metadata, ADRN_Norm_vs_Hypo_Down_635.txt ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt <- as.numeric(r2_gse62564_GSVA_Metadata$ADRN_Norm_vs_Hypo_Down_635.txt)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ADRN_Norm_vs_Hypo_Down_635.txt, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
Reference paper for the use of the T cell exhaustion gene set in this section was Chen 2020:
r2_gse62564_GSVA_Metadata$`T cell exhaustion` <- as.numeric(r2_gse62564_GSVA_Metadata$`T cell exhaustion`)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=`T cell exhaustion`, fill=high_risk)) +
geom_boxplot()
The compare_means function can not deal with spaces in the T cell exhaustion variable name.
We will have to use a trick to deal with the spaces in “T cell exhaustion.”
We will use the package stringr to insert dots (.) replacing the spaces in the r2_gse62564_GSVA_Metadata dataframe.
To start, we have to create the r2_gse62564_GSVA_Metadata_Spaces dataframe, which will have structure names that contained spaces, changed to one of these new characteres: “.” , “,” or ““.
r2_gse62564_GSVA_Metadata_Spaces <- r2_gse62564_GSVA_Metadata
library(stringr)
## Warning: package 'stringr' was built under R version 4.1.2
names(r2_gse62564_GSVA_Metadata_Spaces)<-str_replace_all(names(r2_gse62564_GSVA_Metadata_Spaces), c(" " = "." , "," = "" ))
Then, instead of using the r2_gse62564_GSVA_Metadata dataframe, we use r2_gse62564_GSVA_Metadata_Spaces. Note that now we will use the variable name as T.cell.exhaustion, without the spaces.
compare_means(data = r2_gse62564_GSVA_Metadata_Spaces, T.cell.exhaustion ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata_Spaces$T.cell.exhaustion <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$T.cell.exhaustion)
ggplot(r2_gse62564_GSVA_Metadata_Spaces, aes(x=high_risk, y=T.cell.exhaustion, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
Hypoxia induces quiescence, the cell cycle arrest at the G0 phase (Barbara Muz 2015).
We found a quiescence signature on the work that discussed signatures of quiescent stem cells (Cheung and Rando 2013).
HIF1A activation in cancer cells augments T cell exhaustion and reduces T cell killing (Pilar Baldominos 2022).
In another example, we interrogate about the Quiescence gene signatures.
r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung)
compare_means(data = r2_gse62564_GSVA_Metadata_Spaces, Quiescence.Up.Cheung ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Up.Cheung)
ggplot(r2_gse62564_GSVA_Metadata_Spaces, aes(x=high_risk, y=Quiescence.Up.Cheung, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung)
compare_means(data = r2_gse62564_GSVA_Metadata_Spaces, Quiescence.Down.Cheung ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces$Quiescence.Down.Cheung)
ggplot(r2_gse62564_GSVA_Metadata_Spaces, aes(x=high_risk, y=Quiescence.Down.Cheung, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE <- as.numeric(r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_INFLAMMATORY_RESPONSE, fill=high_risk)) +
geom_boxplot()
compare_means(data = r2_gse62564_GSVA_Metadata, HALLMARK_INFLAMMATORY_RESPONSE ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE <- as.numeric(r2_gse62564_GSVA_Metadata$HALLMARK_INFLAMMATORY_RESPONSE)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=HALLMARK_INFLAMMATORY_RESPONSE, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
Reference paper for this section is Mabe 2022:
r2_gse62564_GSVA_Metadata$ST8SIA1 <- as.numeric(r2_gse62564_GSVA_Metadata$ST8SIA1)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ST8SIA1, fill=high_risk)) +
geom_boxplot()
compare_means(data = r2_gse62564_GSVA_Metadata, ST8SIA1 ~ high_risk, method = "t.test")
r2_gse62564_GSVA_Metadata$ST8SIA1 <- as.numeric(r2_gse62564_GSVA_Metadata$ST8SIA1)
ggplot(r2_gse62564_GSVA_Metadata, aes(x=high_risk, y=ST8SIA1, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
library(tidyverse)
r2_gse62564_GSVA_Metadata_survival <- r2_gse62564_GSVA_Metadata %>%
mutate(quantile = ntile(HIF1A, 2))
r2_gse62564_GSVA_Metadata_survival$quantile <- as.factor(r2_gse62564_GSVA_Metadata_survival$quantile )
qplot(HALLMARK_HYPOXIA, HIF1A,
data = r2_gse62564_GSVA_Metadata_survival,
colour=quantile,
ylab = "HIF1A",
xlab = "Hallmark Hypoxia")
library(survival)
fit_os_HIF1A <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata_survival$quantile,
data = r2_gse62564_GSVA_Metadata_survival)
library(survminer)
ggsurvplot(fit_os_HIF1A,
data = r2_gse62564_GSVA_Metadata_survival,
risk.table = TRUE,
# risk.table = "abs_pct",
pval = TRUE,
pval.coord = c(4000, 1),
)
library(tidyverse)
r2_gse62564_GSVA_Metadata_survival <- r2_gse62564_GSVA_Metadata %>%
mutate(quantile = ntile(HIF1A, 2))
library(survival)
fit_efs_HIF1A <- survfit(Surv(as.numeric(efs_day), efs_bin_numeric) ~ r2_gse62564_GSVA_Metadata_survival$quantile,
data = r2_gse62564_GSVA_Metadata_survival)
library(survminer)
ggsurvplot(fit_efs_HIF1A,
data = r2_gse62564_GSVA_Metadata_survival,
risk.table = TRUE,
# risk.table = "abs_pct",
pval = TRUE,
pval.coord = c(4000, 1),
)
library(tidyverse)
r2_gse62564_GSVA_Metadata_survival <- r2_gse62564_GSVA_Metadata %>%
mutate(quantile = ntile(ELVIDGE_HIF1A_TARGETS_UP, 4))
library(survival)
fit_os_ELVIDGE_HIF1A_TARGETS_UP <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata_survival$quantile,
data = r2_gse62564_GSVA_Metadata_survival)
library(survminer)
ggsurvplot(fit_os_ELVIDGE_HIF1A_TARGETS_UP,
data = r2_gse62564_GSVA_Metadata_survival,
risk.table = TRUE,
# risk.table = "abs_pct",
pval = TRUE,
pval.coord = c(4000, 1),
)
Another solution to this is using the ifelse statement to generate the gene set score as a categorical vector.
Once the categorical vector is generated, it can be used to be plotted in the Kaplan Meyer Survival curve.
library(tidyverse)
r2_gse62564_GSVA_Metadata_Spaces_survival <- r2_gse62564_GSVA_Metadata_Spaces %>%
mutate(quantile = ntile(T.cell.exhaustion, 4))
library(survival)
fit_os_T.cell.exhaustion <- survfit(Surv(as.numeric(os_day), os_bin_numeric) ~ r2_gse62564_GSVA_Metadata_Spaces_survival$quantile,
data = r2_gse62564_GSVA_Metadata_Spaces_survival)
library(survminer)
ggsurvplot(fit_os_T.cell.exhaustion,
data = r2_gse62564_GSVA_Metadata_Spaces_survival,
risk.table = TRUE,
# risk.table = "abs_pct",
pval = TRUE,
pval.coord = c(4000, 1),
)
library(tidyverse)
r2_gse62564_GSVA_Metadata_Spaces_survival <- r2_gse62564_GSVA_Metadata_Spaces %>%
mutate(quantile = ntile(T.cell.exhaustion, 2))
library(survival)
fit_efs_T.cell.exhaustion <- survfit(Surv(as.numeric(efs_day), efs_bin_numeric) ~ r2_gse62564_GSVA_Metadata_Spaces_survival$quantile,
data = r2_gse62564_GSVA_Metadata_Spaces_survival)
library(survminer)
ggsurvplot(fit_efs_T.cell.exhaustion,
data = r2_gse62564_GSVA_Metadata_Spaces_survival,
risk.table = TRUE,
# risk.table = "abs_pct",
pval = TRUE,
pval.coord = c(4000, 1),
)
r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis <- as.numeric(r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis)
max(r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis)
## [1] 8983
min(r2_gse62564_GSVA_Metadata_Spaces_survival$age_at_diagnosis)
## [1] 0
r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival %>%
mutate(age_qtls = ntile(age_at_diagnosis, 4))
r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls <- as.factor(r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls)
ggplot(r2_gse62564_GSVA_Metadata_Spaces_survival_age, aes(x=age_qtls, y=HIF1A, fill=high_risk)) +
geom_boxplot()+
stat_compare_means(method = "t.test")+
theme_linedraw()
r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival %>%
mutate(age_qtls = ntile(age_at_diagnosis, 4))
r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls <- as.factor(r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls)
ggplot(r2_gse62564_GSVA_Metadata_Spaces_survival_age, aes(x=age_qtls, y=HIF1A, fill=mycn_status)) +
geom_boxplot()+
stat_compare_means(method = "anova")+
theme_linedraw()
Remove the NAs from the mycn_status variable
r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival_age[which(
r2_gse62564_GSVA_Metadata_Spaces_survival_age$mycn_status != "na" ), ]
r2_gse62564_GSVA_Metadata_Spaces_survival_age <- r2_gse62564_GSVA_Metadata_Spaces_survival_age %>%
mutate(age_qtls = ntile(age_at_diagnosis, 4))
r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls <- as.factor(r2_gse62564_GSVA_Metadata_Spaces_survival_age$age_qtls)
ggplot(r2_gse62564_GSVA_Metadata_Spaces_survival_age, aes(x=age_qtls, y=HIF1A, fill=mycn_status)) +
geom_boxplot()+
stat_compare_means(method = "anova")+
theme_linedraw()
sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lubridate_1.9.2 forcats_1.0.0 dplyr_1.1.2 purrr_1.0.1
## [5] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 tidyverse_2.0.0
## [9] stringr_1.5.0 rstatix_0.7.2 survminer_0.4.9 ggpubr_0.6.0
## [13] survival_3.5-5 ggplot2_3.4.2 maditr_0.8.3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.10 lattice_0.21-8 zoo_1.8-12 digest_0.6.32
## [5] utf8_1.2.3 R6_2.5.1 backports_1.4.1 evaluate_0.21
## [9] highr_0.10 pillar_1.9.0 rlang_1.1.1 rstudioapi_0.14
## [13] data.table_1.14.8 car_3.1-2 jquerylib_0.1.4 Matrix_1.5-1
## [17] rmarkdown_2.22 labeling_0.4.2 splines_4.1.1 munsell_0.5.0
## [21] gridtext_0.1.5 broom_1.0.5 compiler_4.1.1 xfun_0.39
## [25] pkgconfig_2.0.3 htmltools_0.5.5 ggtext_0.1.2 tidyselect_1.2.0
## [29] gridExtra_2.3 km.ci_0.5-6 fansi_1.0.4 tzdb_0.4.0
## [33] withr_2.5.0 commonmark_1.9.0 grid_4.1.1 jsonlite_1.8.7
## [37] xtable_1.8-4 gtable_0.3.3 lifecycle_1.0.3 magrittr_2.0.3
## [41] KMsurv_0.1-5 scales_1.2.1 stringi_1.7.12 cli_3.6.1
## [45] cachem_1.0.8 carData_3.0-5 farver_2.1.1 ggsignif_0.6.4
## [49] xml2_1.3.4 bslib_0.5.0 survMisc_0.5.6 generics_0.1.3
## [53] vctrs_0.6.3 tools_4.1.1 glue_1.6.2 markdown_1.7
## [57] hms_1.1.3 abind_1.4-5 fastmap_1.1.1 yaml_2.3.7
## [61] timechange_0.2.0 colorspace_2.1-0 knitr_1.43 sass_0.4.6