── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(survival)library(survminer)
Loading required package: ggpubr
Attaching package: 'survminer'
The following object is masked from 'package:survival':
myeloma
library(car)
Loading required package: carData
Attaching package: 'car'
The following object is masked from 'package:dplyr':
recode
The following object is masked from 'package:purrr':
some
data(trial)
trial <- trial %>%# manually re-name columns# NEW name # OLD namerename(Treatment = trt) %>%# de-duplicatedistinct()trial
# A tibble: 200 × 8
Treatment age marker stage grade response death ttdeath
<chr> <dbl> <dbl> <fct> <fct> <int> <int> <dbl>
1 Drug A 23 0.16 T1 II 0 0 24
2 Drug B 9 1.11 T2 I 1 0 24
3 Drug A 31 0.277 T1 II 0 0 24
4 Drug A NA 2.07 T3 III 1 1 17.6
5 Drug A 51 2.77 T4 III 1 1 16.4
6 Drug B 39 0.613 T4 I 0 1 15.6
7 Drug A 37 0.354 T1 II 0 0 24
8 Drug A 32 1.74 T1 I 0 1 18.4
9 Drug A 31 0.144 T1 II 0 0 24
10 Drug B 34 0.205 T3 I 0 1 10.5
# ℹ 190 more rows
The following objects are masked from 'package:dplyr':
src, summarize
The following objects are masked from 'package:base':
format.pval, units
Attaching package: 'rms'
The following objects are masked from 'package:car':
Predict, vif
dd <-datadist(trial) options(datadist ="dd")
logistic1 <-glm(response ~rcs(age, 3) + Treatment + marker + stage + grade, data = trial, family ="binomial")
logistic2 <-glm(response ~ age + Treatment +rcs(marker, 3) + stage + grade, data = trial, family ="binomial")
Comparing the cubic spline model for the marker with the full model
anova(model_response_multi_full, logistic2, test ="LRT")
Analysis of Deviance Table
Model 1: response ~ Treatment + age + marker + stage + grade
Model 2: response ~ age + Treatment + rcs(marker, 3) + stage + grade
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 164 205.89
2 163 205.89 1 0.0054495 0.9412
The difference in deviance between the models is very small (0.0054). The non-linear spline term for marker does not significantly improve the model (p = 0.9412), so the linearity assumption for marker holds.
Comparing the cubic spline model for age with the full model
anova(model_response_multi_full,logistic1, test ="LRT")
Analysis of Deviance Table
Model 1: response ~ Treatment + age + marker + stage + grade
Model 2: response ~ rcs(age, 3) + Treatment + marker + stage + grade
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 164 205.89
2 163 205.88 1 0.010616 0.9179
The difference in deviance between the models is very small (0.010616). The non-linear spline term for age does not significantly improve the model (p = 0.9179), so the linearity assumption for age holds.
Backward elimination (threshold p > 0.10) to validate the analysis of deviance
model_response_multi_step <-step(model_response_multi_full, direction ="backward", trace =FALSE)summary(model_response_multi_step)
Call:
glm(formula = response ~ age + marker, family = "binomial", data = trial)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.95041 0.62357 -3.128 0.00176 **
age 0.01880 0.01172 1.604 0.10882
marker 0.27647 0.18640 1.483 0.13801
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 214.80 on 172 degrees of freedom
Residual deviance: 209.97 on 170 degrees of freedom
(27 observations deleted due to missingness)
AIC: 215.97
Number of Fisher Scoring iterations: 4
Multi-collinearity for the Cox model
vif(model_survival_multi_full)
TreatmentDrug B age marker stageT2 stageT3
1.010961 1.057127 1.069657 1.845739 1.682669
stageT4 gradeII gradeIII
1.949996 1.392637 1.371761
Analysis of deviance to assess predictor variables that improve the Cox regression model
(anova(coxph(surv_object ~ Treatment + age + marker + stage + grade, data = trial)))
tbl_merge_1 <-tbl_merge(tbls =list(tbl1, tbl2),tab_spanner =c("**Tumor Response**", "**Time to Death**")) %>%bold_labels()%>%modify_spanning_header(all_stat_cols() ~"**Treatment group**") %>%as_gt() %>% gt::tab_header(title = gt::md("**Table 2. Associations of Clinical and Treatment Factors with Tumor Response and Time to Death**"),subtitle = gt::md("_Highly Confidential_") ) %>% gt::tab_source_note("Concordance= 0.637 (se = 0.03 )Likelihood ratio test= 25.18 on 8 df, p=0.001Wald test = 25.43 on 8 df, p=0.001Score (logrank) test = 26.7 on 8 df, p=8e-04")tbl_merge_1
Table 2. Associations of Clinical and Treatment Factors with Tumor Response and Time to Death
Highly Confidential
Characteristic
Tumor Response
Time to Death
OR
95% CI
p-value
HR
95% CI
p-value
Chemotherapy Treatment
0.4
0.11
Drug A
—
—
—
—
Drug B
1.34
0.69, 2.64
1.39
0.93, 2.10
Age
1.02
1.00, 1.04
0.11
1.01
1.00, 1.03
0.2
Grade
>0.9
0.086
I
—
—
—
—
II
1.04
0.45, 2.42
1.15
0.68, 1.95
III
1.04
0.47, 2.32
1.69
1.05, 2.74
T Stage
0.3
<0.001
T1
—
—
—
—
T2
0.45
0.17, 1.11
1.52
0.82, 2.81
T3
0.86
0.33, 2.22
1.42
0.73, 2.77
T4
0.64
0.25, 1.62
3.23
1.79, 5.86
Marker Level (ng/mL)
1.41
0.96, 2.09
0.080
0.89
0.69, 1.14
0.3
Abbreviations: CI = Confidence Interval, HR = Hazard Ratio, OR = Odds Ratio
Concordance= 0.637 (se = 0.03 ) Likelihood ratio test= 25.18 on 8 df, p=0.001 Wald test = 25.43 on 8 df, p=0.001 Score (logrank) test = 26.7 on 8 df, p=8e-04
Schoenfeld residuals for Cox regression assumptions
