-
Notifications
You must be signed in to change notification settings - Fork 0
/
figure4.R
138 lines (115 loc) · 3.45 KB
/
figure4.R
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
128
129
130
131
132
133
134
135
136
137
138
# --------------------------------------------------------------------------------- #
# Neighborhood-Level Deprivation and Survival in Lung Cancer
# --------------------------------------------------------------------------------- #
#
# Figure 4: High NDI is associated with increased risk of mortality among lung
# cancer patients
#
# Created by: Ian D. Buller, Ph.D., M.A. (@idblr)
# Created on: 2022-11-07
#
# Most recently modified by: @idblr
# Most recently modified on: 2024-08-06
#
# Notes:
# A) 2022-10-30 (@idblr): Initial script created by Ignacio Jusué-Torres, MD
# B) 2023-04-26 (@idblr): Updated script created by Ignacio Jusué-Torres, MD
# C) 2023-07-09 (@idblr): Relabeled variables used in Cox models for figures
# D) 2024-07-10 (@idblr): Re-formatted code
# E) 2024-08-06 (@idblr): Re-formatted code
# --------------------------------------------------------------------------------- #
# ---------------- #
# DATA IMPORTATION #
# ---------------- #
source(file.path('code', 'preparation.R'))
# -------------------- #
# ADDITIONAL LIBRARIES #
# -------------------- #
loadedPackages <- c('grDevices', 'forestmodel', 'stats', 'survival')
suppressMessages(invisible(lapply(loadedPackages, library, character.only = TRUE)))
# -------- #
# FIGURE 4 #
# -------- #
# FIGURE 4A
## Multivariable Cox Proportional Hazards
## NDI (Messer) as factor (quartiles)
## Age-adjusted
summary(CANCER$SurvTime)
table(CANCER$SurvivalStatus, CANCER$Death)
# Label preparation
CANCER <- CANCER %>% rename(`Age at surgery` = 'Age.at.surgery')
CANCER <- CANCER %>% rename(NDI = 'NDImesser_qt')
# Model
f4a <- coxph(
Surv(SurvTime, SurvivalStatus) ~ `Age at surgery` + NDI,
data = CANCER)
summary(f4a)
# Test for PH assumptions
## Scaled Schoenfeld residuals independent of time?
cox.zph(f4a)
# Test for linearity (p-trend)
f4a_ <- coxph(
Surv(SurvTime, SurvivalStatus) ~ `Age at surgery` + NDImesser_qt.1,
data = CANCER
)
summary(f4a_) # numeric variable is significant
anova(f4a, f4a_) # not linear
f4a__ <- coxph(
Surv(SurvTime, SurvivalStatus) ~ `Age at surgery` + NDImesser_tr,
data = CANCER
)
summary(f4a__) # numeric variable is *not* significant
# Plot
p4a <- forest_model(
f4a,
format_options = list(
colour = 'black',
shape = 15,
banded = TRUE,
point_size = 4
)
)
png(file.path('figures', 'figure4a.png'), height = 800, width = 1100)
p4a
dev.off()
# FIGURE 4B
## Multivariable Cox Proportional Hazards
## Stage + Race + HOXA7 + NDI (Messer) dichotomized (1-3, 4)
## Age-adjusted
# Label preparation
CANCER <- CANCER %>% rename(NDImesser_qt = 'NDI')
CANCER <- CANCER %>%
mutate(Race = fct_relevel(Race, c('W', 'B', 'O')))
levels(CANCER$Race) <- c('white', 'Black', 'other')
CANCER$Stage <- droplevels(CANCER$Stage)
CANCER <- CANCER %>% rename(NDI = 'NDImesser_qt_d')
levels(CANCER$NDI) <- c('Low deprivation', 'High deprivation')
CANCER <- CANCER %>% rename(`HOXA7 methylation` = 'HOXA7.Plasma')
# Model
f4b <- coxph(
Surv(SurvTime, SurvivalStatus) ~
`Age at surgery` +
Race +
Stage +
`HOXA7 methylation` +
NDI,
data = CANCER
)
summary(f4b)
# Test for PH assumptions
## Scaled Schoenfeld residuals independent of time?
cox.zph(f4b)
# Plot
p4b <- forest_model(
f4b,
format_options = list(
colour = 'black',
shape = 15,
banded = TRUE,
point_size = 4
)
)
png(file.path('figures', 'figure4b.png'), height = 800, width = 1100)
p4b
dev.off()
# ---------------------------------- END OF CODE ---------------------------------- #