1 Raw distributions: Demographics/CBC

2 Log distributions: Demo/CBC variables

3 KP/Cytokine inspection (RAW)

4 KP/Cytokines (log transformations)

5 Merging with log-transformed

6 Clinical outcome inspection (HAMD17)

SG_df_new_long <- read_excel("~/Desktop/SG_SII_SIRI/SG_df_new_long_05082023.xlsx")
SG_df_new_long <- SG_df_new_long %>% subset(Pt_Group=="TRBDD")
SG_df_new_long$Pt_Group<-as.factor(SG_df_new_long$Pt_Group)
SG_df_new_long$Treatment<-as.factor(SG_df_new_long$Treatment)
SG_df_new_long$Sex<-as.factor(SG_df_new_long$Sex)
SG_df_new_long$Timepoint<-as.factor(SG_df_new_long$Timepoint)

ggplot(SG_df_new_long, aes(x = Timepoint, y = log(HAMD17_total)))+
  geom_boxplot(aes(fill=Timepoint))+
 geom_jitter(width = 0.1)+
  facet_wrap(~Treatment)+
  theme_bw()+   
  theme(legend.position = "none")+
  ggpubr::stat_compare_means(method="t.test", label.y=4)
## Warning: Removed 2 rows containing non-finite values (stat_boxplot).
## Warning: Removed 2 rows containing non-finite values (stat_compare_means).
## Warning: Removed 2 rows containing missing values (geom_point).

7 TABLE 1A

#TABLE 1A

cont_vars<-c(
  "log_MONO_BL")


combined_df %>%
  dplyr::select(Cohort,
                Sex,
                Age,
                log_BMI,
                log_NEUT_BL,
                log_MONO_BL,
                log_LYMPH_BL,
                log_SIRI_BL) %>%
  tbl_summary(by = Cohort, 
              type=list(cont_vars~"continuous"),
              missing="no") %>%
  add_p(pvalue_fun = ~ style_pvalue(.x, digits = 2)) %>%
  add_overall() %>%
  add_n() %>%
  modify_header(label ~ "**Variable**") %>%
  # modify_spanning_header(c("stat_1", "stat_2") ~ "**Treatment Arm**") %>%
  modify_footnote(
    all_stat_cols() ~ "Median (IQR) or Frequency (%)") %>%
  modify_caption("**Table 1. Demographics and CBC markers by Group**") %>%
  bold_labels() %>% 
  bold_p(t=0.1) %>% 
  add_q(method="bonferroni")
## 2 observations missing `Cohort` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `Cohort` column before passing to `tbl_summary()`.
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(cont_vars)` instead of `cont_vars` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## add_q: Adjusting p-values with
## `stats::p.adjust(x$table_body$p.value, method = "bonferroni")`
Table 1. Demographics and CBC markers by Group
Variable N Overall, N = 841 CBX_ESC, N = 291 HC, N = 321 PBO_ESC, N = 231 p-value2 q-value3
Sex 83 0.016 0.11
Male 44 (53%) 17 (59%) 11 (34%) 16 (73%)
Female 39 (47%) 12 (41%) 21 (66%) 6 (27%)
Age 74 40 (31, 52) 38 (31, 44) 37 (26, 53) 47 (35, 58) 0.083 0.58
log_BMI 71 3.37 (3.23, 3.49) 3.39 (3.29, 3.57) 3.20 (3.12, 3.34) 3.44 (3.33, 3.60) <0.001 0.005
log_NEUT_BL 81 1.25 (1.03, 1.57) 1.28 (1.08, 1.49) 1.15 (1.02, 1.53) 1.50 (1.06, 1.68) 0.30 >0.99
log_MONO_BL 82 -0.92 (-0.92, -0.51) -0.69 (-0.92, -0.51) -0.92 (-0.92, -0.65) -0.69 (-1.13, -0.40) 0.81 >0.99
log_LYMPH_BL 82 0.59 (0.47, 0.79) 0.59 (0.52, 0.79) 0.59 (0.52, 0.83) 0.59 (0.47, 0.79) 0.90 >0.99
log_SIRI_BL 81 -0.16 (-0.51, 0.24) -0.03 (-0.52, 0.20) -0.28 (-0.54, 0.19) 0.04 (-0.32, 0.40) 0.41 >0.99
1 Median (IQR) or Frequency (%)
2 Pearson's Chi-squared test; Kruskal-Wallis rank sum test
3 Bonferroni correction for multiple testing 
# %>%
#   gtsummary::as_tibble() %>%
#   writexl::write_xlsx(., "Table 1. Patient Characteristics (Whole Cohort)")

#Table 1B

#TABLE 1B

cont_vars<-c(
"log_IL1A_BL",
"log_IL1B_BL",
"log_IFNG_BL")

combined_df %>%
  dplyr::select(-Subject_ID, 
                -Group, 
                -Cohort, 
                -log_SII_BL, 
                -log_NEUT_BL, 
                -log_PLT_BL,
                -log_MONO_BL, 
                -log_LYMPH_BL,
                -contains("WK8"), 
                -contains("KP"),
                log_KP_Trp_BL,
                log_KP_Kyn_BL,
                log_KP_Quin_BL,
                log_KP_Pic_BL,
                log_KP_KYN_TRP_BL,
                log_KP_QUIN_PIC_BL,
                HAMD17_WK8) %>%
  tbl_summary(by = Arm, 
              type=list(cont_vars~"continuous"),
              missing="no") %>%
  add_p(pvalue_fun = ~ style_pvalue(.x, digits = 2)) %>%
  add_overall() %>%
  add_n() %>%
  modify_header(label ~ "**Variable**") %>%
  # modify_spanning_header(c("stat_1", "stat_2") ~ "**Treatment Arm**") %>%
  modify_footnote(
    all_stat_cols() ~ "Median (IQR) or Frequency (%)") %>%
  modify_caption("**Table 1. Patient Characteristics (Whole Cohort)**") %>%
  bold_labels() %>% 
  bold_p(t=0.1) %>% 
  add_q(method="bonferroni")
## 34 observations missing `Arm` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `Arm` column before passing to `tbl_summary()`.
## Warning for variable 'Age':
## simpleWarning in wilcox.test.default(x = c(61, 32, 64, 64, 36, 50, 57, 59, 62, : cannot compute exact p-value with ties
## Warning for variable 'HAMD17_BL':
## simpleWarning in wilcox.test.default(x = c(21, 23, 33, 27, 26, 35, 23, 15, 26, : cannot compute exact p-value with ties
## Warning for variable 'log_BMI':
## simpleWarning in wilcox.test.default(x = c(3.68637632389582, 3.61899332664977, : cannot compute exact p-value with ties
## Warning for variable 'log_SIRI_BL':
## simpleWarning in wilcox.test.default(x = c(-0.341749293722057, -1.01291756756323, : cannot compute exact p-value with ties
## Warning for variable 'log_IL1A_BL':
## simpleWarning in wilcox.test.default(x = c(0.792992515529661, 0.693147180559945, : cannot compute exact p-value with ties
## Warning for variable 'log_IL1B_BL':
## simpleWarning in wilcox.test.default(x = c(0.989541193613748, 0.693147180559945, : cannot compute exact p-value with ties
## Warning for variable 'log_IL6_BL':
## simpleWarning in wilcox.test.default(x = c(1.68268837417369, 1.72988406550997, : cannot compute exact p-value with ties
## Warning for variable 'log_TNFA_BL':
## simpleWarning in wilcox.test.default(x = c(1.39871688111845, 1.43983512804792, : cannot compute exact p-value with ties
## Warning for variable 'log_IFNG_BL':
## simpleWarning in wilcox.test.default(x = c(1.19694818938897, 0.693147180559945, : cannot compute exact p-value with ties
## Warning for variable 'log_FGF_BL':
## simpleWarning in wilcox.test.default(x = c(1.3609765531356, 1.58923520511658, : cannot compute exact p-value with ties
## Warning for variable 'log_EGF_BL':
## simpleWarning in wilcox.test.default(x = c(1.32175583998232, 1.53901544813755, : cannot compute exact p-value with ties
## Warning for variable 'log_KP_Quin_BL':
## simpleWarning in wilcox.test.default(x = c(4.32809829264833, 3.8177123259569, : cannot compute exact p-value with ties
## Warning for variable 'log_KP_Pic_BL':
## simpleWarning in wilcox.test.default(x = c(3.45315712059287, 3.16968558067743, : cannot compute exact p-value with ties
## Warning for variable 'HAMD17_WK8':
## simpleWarning in wilcox.test.default(x = c(19, 11, 33, 33, 9, 8, 17, 9, 10, 21, : cannot compute exact p-value with ties
## add_q: Adjusting p-values with
## `stats::p.adjust(x$table_body$p.value, method = "bonferroni")`
Table 1. Patient Characteristics (Whole Cohort)
Variable N Overall, N = 521 PBO_ESC, N = 231 CBX_ESC, N = 291 p-value2 q-value3
Sex 51 0.30 >0.99
Male 33 (65%) 16 (73%) 17 (59%)
Female 18 (35%) 6 (27%) 12 (41%)
Age 51 40 (34, 50) 47 (35, 58) 38 (31, 44) 0.042 0.92
HAMD17_BL 45 24.0 (20.0, 29.0) 23.5 (21.0, 26.0) 24.0 (20.0, 30.0) 0.67 >0.99
log_BMI 50 3.40 (3.32, 3.57) 3.44 (3.33, 3.60) 3.39 (3.29, 3.57) 0.45 >0.99
log_SIRI_BL 49 -0.03 (-0.45, 0.32) 0.04 (-0.32, 0.40) -0.03 (-0.52, 0.20) 0.51 >0.99
log_IL1A_BL 21 0.69 (0.69, 0.79) 0.74 (0.69, 0.80) 0.69 (0.69, 0.69) 0.24 >0.99
log_IL1B_BL 22 0.69 (0.69, 0.88) 0.69 (0.69, 0.96) 0.69 (0.69, 0.69) 0.46 >0.99
log_IL6_BL 21 1.25 (1.00, 1.44) 1.56 (1.01, 1.70) 1.22 (1.00, 1.34) 0.19 >0.99
log_TNFA_BL 21 1.18 (1.11, 1.40) 1.21 (1.12, 1.41) 1.18 (1.08, 1.30) 0.69 >0.99
log_IFNG_BL 21 0.69 (0.69, 0.69) 0.69 (0.69, 0.82) 0.69 (0.69, 0.69) >0.99 >0.99
log_CRPSet1_ug_ml_BL 15 1.61 (1.16, 1.95) 1.58 (1.22, 1.65) 1.95 (1.19, 2.58) 0.34 >0.99
log_MCP1_BL 21 4.59 (4.48, 4.85) 4.76 (4.47, 4.90) 4.56 (4.48, 4.71) 0.50 >0.99
log_FGF_BL 15 1.18 (0.88, 1.58) 1.36 (1.08, 1.57) 1.11 (0.76, 1.58) 0.77 >0.99
log_VEGF-Elisa_BL 31 3.49 (3.36, 3.81) 3.41 (3.30, 3.49) 3.63 (3.48, 3.92) 0.059 >0.99
log_EGF_BL 21 1.28 (0.69, 1.71) 1.43 (1.13, 1.61) 1.20 (0.69, 1.71) 0.48 >0.99
log_KP_Trp_BL 41 9.70 (9.51, 9.79) 9.59 (9.45, 9.76) 9.73 (9.58, 9.80) 0.16 >0.99
log_KP_Kyn_BL 41 5.76 (5.56, 5.94) 5.74 (5.64, 6.07) 5.76 (5.55, 5.86) 0.72 >0.99
log_KP_Quin_BL 41 4.02 (3.76, 4.29) 4.04 (3.83, 4.37) 3.98 (3.64, 4.25) 0.39 >0.99
log_KP_Pic_BL 41 3.17 (2.84, 3.33) 3.20 (2.69, 3.36) 3.00 (2.84, 3.33) 0.98 >0.99
log_KP_KYN_TRP_BL 41 0.703 (0.701, 0.706) 0.705 (0.702, 0.707) 0.702 (0.701, 0.706) 0.26 >0.99
log_KP_QUIN_PIC_BL 41 1.56 (1.41, 1.76) 1.55 (1.47, 1.71) 1.56 (1.39, 1.76) 0.64 >0.99
HAMD17_WK8 45 10 (7, 17) 12 (9, 18) 8 (5, 13) 0.008 0.17
1 Median (IQR) or Frequency (%)
2 Pearson's Chi-squared test; Wilcoxon rank sum test; Wilcoxon rank sum exact test
3 Bonferroni correction for multiple testing 
combined_df$Cohort<-factor(combined_df$Cohort, levels=c("HC", "PBO_ESC", "CBX_ESC"))

8 Covariate for SIRI_BL

list_vars<-combined_df %>% dplyr::select(-Subject_ID, 
                -Group, 
                -Cohort, 
                -log_SII_BL, 
                -log_NEUT_BL, 
                -log_PLT_BL,
                -log_MONO_BL, 
                -log_LYMPH_BL,
                -contains("WK8"), 
                -contains("KP"),
                HAMD17_WK8,
                log_KP_Trp_BL,
                log_KP_Kyn_BL,
                log_KP_Quin_BL,
                log_KP_Pic_BL,
                log_KP_KYN_TRP_BL,
                log_KP_QUIN_PIC_BL) %>% 
  names()

heatmap_df<-combined_df %>% dplyr::select(all_of(list_vars))
heatmap_df <- heatmap_df[which(heatmap_df$Arm!=c("HC")), ] %>% select(-Arm)

gtsummary::tbl_uvregression(
  heatmap_df,
  lm,
  log_SIRI_BL,
  conf.int=TRUE
) %>% 
  bold_p(t=0.1) %>% 
  add_global_p() %>% 
  add_q(method="bonferroni")
## add_q: Adjusting p-values with
## `stats::p.adjust(x$table_body$p.value, method = "bonferroni")`
Characteristic N Beta 95% CI1 p-value q-value2
Sex 48 0.8 >0.9
Male — —
Female 0.04 -0.35, 0.42
Age 48 -0.01 -0.02, 0.01 0.2 >0.9
HAMD17_BL 43 -0.03 -0.06, 0.00 0.055 >0.9
log_BMI 47 -0.55 -1.5, 0.44 0.3 >0.9
log_IL1A_BL 21 -2.6 -6.4, 1.2 0.2 >0.9
log_IL1B_BL 22 -0.80 -2.3, 0.66 0.3 >0.9
log_IL6_BL 21 -0.30 -0.93, 0.33 0.3 >0.9
log_TNFA_BL 21 0.06 -0.27, 0.38 0.7 >0.9
log_IFNG_BL 21 -0.01 -0.79, 0.77 >0.9 >0.9
log_CRPSet1_ug_ml_BL 15 0.15 -0.35, 0.66 0.5 >0.9
log_MCP1_BL 21 -0.92 -1.5, -0.35 0.003 0.067
log_FGF_BL 15 0.16 -0.63, 0.95 0.7 >0.9
log_VEGF-Elisa_BL 31 -0.31 -0.94, 0.32 0.3 >0.9
log_EGF_BL 21 -0.17 -0.69, 0.35 0.5 >0.9
HAMD17_WK8 43 0.01 -0.01, 0.03 0.4 >0.9
log_KP_Trp_BL 39 -0.11 -1.0, 0.76 0.8 >0.9
log_KP_Kyn_BL 39 0.27 -0.35, 0.90 0.4 >0.9
log_KP_Quin_BL 39 0.32 -0.12, 0.76 0.15 >0.9
log_KP_Pic_BL 39 0.26 -0.18, 0.69 0.2 >0.9
log_KP_KYN_TRP_BL 39 16 -31, 63 0.5 >0.9
log_KP_QUIN_PIC_BL 39 0.11 -0.40, 0.61 0.7 >0.9
1 CI = Confidence Interval
2 Bonferroni correction for multiple testing

9 MODEL 1: HAMD17_WK8 by SIRI_BL

SIRI_model<-lm(HAMD17_WK8~
                 Sex+
                 Age+
                 log_BMI+
                 Arm+
                 HAMD17_BL+
                 log_SIRI_BL
               , data=combined_df)
sjPlot::tab_model(SIRI_model)
  HAMD 17 WK 8
Predictors Estimates CI p
(Intercept) -15.64 -56.10 – 24.83 0.438
Sex [Female] 4.39 0.15 – 8.63 0.043
Age 0.12 -0.06 – 0.29 0.198
log BMI 4.61 -6.91 – 16.14 0.422
Arm [CBX ESC] -5.30 -9.89 – -0.71 0.025
HAMD17 BL 0.34 -0.02 – 0.69 0.061
log SIRI BL 2.48 -1.34 – 6.29 0.196
Observations 43
R2 / R2 adjusted 0.352 / 0.244

10 MODEL 2: HAMD17_WK8 by SIRI_BL (interaction screen)

SIRI_model_interaction<-lm(HAMD17_WK8~
                             Sex+
                             Age+
                             log_BMI+
                             Arm+
                             HAMD17_BL*log_SIRI_BL, data=combined_df)

sjPlot::tab_model(SIRI_model_interaction)
  HAMD 17 WK 8
Predictors Estimates CI p
(Intercept) -10.29 -47.75 – 27.17 0.581
Sex [Female] 6.02 1.93 – 10.11 0.005
Age 0.05 -0.12 – 0.22 0.525
log BMI 3.31 -7.35 – 13.96 0.533
Arm [CBX ESC] -5.13 -9.36 – -0.90 0.019
HAMD17 BL 0.42 0.09 – 0.75 0.014
log SIRI BL -17.85 -33.26 – -2.44 0.024
HAMD17 BL * log SIRI BL 0.91 0.24 – 1.59 0.009
Observations 43
R2 / R2 adjusted 0.467 / 0.361 
interactions::interact_plot(SIRI_model_interaction, pred = log_SIRI_BL, modx = HAMD17_BL, jitter=0.1, plot.points = TRUE,  main.title = "HAMD17 (week 8) linked to baseline SIRI-to-HAMD17_BL interaction")

# base::plot(SII_model_interaction, which=c(2,6))

#Post-hoc power analysis (valid N=43, power=0.8)

 pwr::pwr.r.test(n=43, r=NULL, power=0.8)
## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 43
##               r = 0.4124439
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided

#Power calc from Model 2 (using adj R2)

 pwr::pwr.r.test(n=43, r=0.361, power=NULL)
## 
##      approximate correlation power calculation (arctangh transformation) 
## 
##               n = 43
##               r = 0.361
##       sig.level = 0.05
##           power = 0.6754509
##     alternative = two.sided

```