master_thesis_code/05_coxph_rcs.R at main · lict99/master_thesis_code · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
# %%
# Attaching packages and functions
library("readr")
library("dplyr")
library("ggplot2")
library("showtext")

source("functions/font_config.R", local = TRUE)
source("functions/rcs_coxph.R", local = TRUE)

showtext_auto()

# %%
# Setting up output directory
output_dir <- "results/05"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

# %%
# Reading UK Biobank data
ukb_data <- read_csv("results/00/ukb_data.csv") |>
  mutate(
    sex = factor(sex, levels = c("female", "male")),
    smoking_status = factor(smoking_status, levels = c("never", "ever")),
    alcohol_drinker_status = factor(
      alcohol_drinker_status,
      levels = c("never", "ever")
    )
  ) |>
  as.data.frame()

# %%
# Plotting restricted cubic splines for UK Biobank data
for (event in c("os", "css")) {
  for (covar_set in c("set1", "set2")) {
    char_event <- switch(
      event,
      os = "总体生存期",
      css = "癌症特异性生存期",
      stop("Unknown event", call. = FALSE)
    )
    covars <- switch(
      covar_set,
      set1 = c("age_at_diagnosis", "sex"),
      set2 = c(
        "age_at_diagnosis",
        "sex",
        "body_mass_index",
        "smoking_status",
        "alcohol_drinker_status"
      ),
      stop("Unknown covariate set", call. = FALSE)
    )
    ukb_p <- plot_coxph_rcs(
      data = ukb_data,
      time = paste0(event, "_time"),
      event = event,
      target = "platelet_count",
      covariates = covars,
      font_family = font_zh,
      xlab = "血小板计数 (10<sup>9</sup>/L)",
      ylab = sprintf("风险比 (%s)", char_event)
    )
    ggsave(
      file.path(output_dir, sprintf("ukb_%s_%s_rcs.pdf", event, covar_set)),
      plot = ukb_p,
      width = 5,
      height = 5
    )
  }
}

# %%
# Reading West China data
hx_data <- read_csv("results/01/hx_data.csv") |>
  mutate(
    sex = factor(sex, levels = c("female", "male")),
    neo_adjuvant_therapy = factor(
      neo_adjuvant_therapy,
      levels = c("no", "yes")
    ),
    smoking = factor(smoking, levels = c("never", "ever")),
    alcohol = factor(alcohol, levels = c("never", "ever"))
  ) |>
  as.data.frame()

# %%
# Plotting restricted cubic splines for West China data
for (event in c("os", "css", "dfs")) {
  for (covar_set in c("set1", "set2")) {
    char_event <- switch(
      event,
      os = "总体生存期",
      css = "癌症特异性生存期",
      dfs = "无疾病生存期",
      stop("Unknown event", call. = FALSE)
    )
    covars <- switch(
      covar_set,
      set1 = c("age", "sex"),
      set2 = c(
        "age",
        "sex",
        "body_mass_index",
        "smoking",
        "alcohol",
        "neo_adjuvant_therapy"
      ),
      stop("Unknown covariate set", call. = FALSE)
    )
    hx_p <- plot_coxph_rcs(
      data = hx_data,
      time = paste0(event, "_time"),
      event = event,
      target = "platelet_count",
      covariates = covars,
      font_family = font_zh,
      xlab = "血小板计数 (10<sup>9</sup>/L)",
      ylab = sprintf("风险比 (%s)", char_event)
    )
    ggsave(
      file.path(output_dir, sprintf("hx_%s_%s_rcs.pdf", event, covar_set)),
      plot = hx_p,
      width = 5,
      height = 5
    )
  }
}