Inflammatory bowel disease

Author
Affiliations

Nathan Constantine-Cooke

Institute of Genetics and Cancer

Published

January 13, 2026

Introduction

Code
source("Survival/utils.R")

# Setup analysis environment
analysis_setup <- setup_analysis()
paths <- analysis_setup$paths
demo <- analysis_setup$demo

demo$FC <- log(demo$FC)

flare.df <- readRDS(paste0(paths$outdir, "flares-demographics.RDS"))
flare.cd.df <- readRDS(paste0(paths$outdir, "flares-demographics-cd.RDS"))
flare.uc.df <- readRDS(paste0(paths$outdir, "flares-demographics-uc.RDS"))

cd.clin.forest <- readRDS(paste0(paths$outdir, "cd-clin-demographics.RDS"))
cd.hard.forest <- readRDS(paste0(paths$outdir, "cd-hard-demographics.RDS"))
uc.clin.forest <- readRDS(paste0(paths$outdir, "uc-clin-demographics.RDS"))
uc.hard.forest <- readRDS(paste0(paths$outdir, "uc-hard-demographics.RDS"))

IBD Control-8

Crohn’s disease

Code
# Categorize IBD Control-8
flare.cd.df <- categorize_variable(flare.cd.df, "control_8",
                                  breaks = c(0, 13, 16),
                                  labels = c("0-12", "13-16"))

# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "control_8",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "IBD Control-8",
  plot_base_path = "plots/cd/soft-flare/ibd/control-8",
  break_time_by = 200,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for continuous control_8 variable
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + IMD + cat + control_8 + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.clin.forest <- rbind(cd.clin.forest, get_HR(fit.me, "control_8"))

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/control-8.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/control-8"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.9227 1.4759 2.5049 0.0000
IMD2 0.8074 0.4922 1.3246 0.3970
IMD3 0.8700 0.5315 1.4241 0.5796
IMD4 0.8528 0.5286 1.3758 0.5141
IMD5 0.9833 0.6177 1.5655 0.9435
catFC 50-250 1.4397 1.0987 1.8866 0.0082
catFC > 250 2.4014 1.7740 3.2506 0.0000
control_8 0.8902 0.8606 0.9207 0.0000
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0012 0.9903 0.9717
IMD 6.0925 3.9280 0.1855
cat 1.3577 1.9739 0.5009
control_8 8.4939 0.9670 0.0034
GLOBAL 15.1377 14.9466 0.4376 
Warning: `gather_()` was deprecated in tidyr 1.2.0.
ℹ Please use `gather()` instead.
ℹ The deprecated feature was likely used in the survminer package.
  Please report the issue at <https://github.com/kassambara/survminer/issues>.
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/control-8-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/control-8-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "control_8",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "IBD Control-8",
  plot_base_path = "plots/cd/hard-flare/ibd/control-8",
  break_time_by = 500,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for continuous control_8 variable
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + IMD + cat + control_8 + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.hard.forest <- rbind(cd.hard.forest, get_HR(fit.me, "control_8"))

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/control-8.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/control-8"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3756 1.0112 1.8714 0.0423
IMD2 0.6980 0.3790 1.2855 0.2485
IMD3 0.8591 0.4697 1.5715 0.6221
IMD4 0.7891 0.4392 1.4178 0.4282
IMD5 0.7981 0.4509 1.4124 0.4387
catFC 50-250 1.9417 1.3669 2.7584 0.0002
catFC > 250 3.5422 2.4220 5.1804 0.0000
control_8 0.9792 0.9359 1.0245 0.3621
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 1.5133 0.9850 0.2149
IMD 8.1252 3.9413 0.0840
cat 7.2016 1.9860 0.0269
control_8 2.8005 0.9806 0.0918
GLOBAL 20.3170 15.1753 0.1677 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/control-8-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/control-8-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Ulcerative colitis

Code
# Categorize IBD Control-8
flare.uc.df <- categorize_variable(flare.uc.df, "control_8",
                                  breaks = c(0, 13, 16),
                                  labels = c("0-12", "13-16"))

# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "control_8",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "IBD Control-8",
  plot_base_path = "plots/uc/soft-flare/ibd/control-8",
  break_time_by = 200,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for continuous control_8 variable
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + IMD + cat + control_8 + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

uc.clin.forest <- rbind(uc.clin.forest, get_HR(fit.me, "control_8"))

# Display plot and model summary
knitr::include_graphics("plots/uc/soft-flare/ibd/control-8.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/soft-flare/ibd/control-8"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4904 1.1840 1.8761 0.0007
IMD2 1.3212 0.7897 2.2106 0.2888
IMD3 1.2137 0.7394 1.9922 0.4437
IMD4 1.4949 0.9289 2.4059 0.0977
IMD5 1.3596 0.8507 2.1729 0.1992
catFC 50-250 1.5424 1.1893 2.0002 0.0011
catFC > 250 2.0139 1.5095 2.6868 0.0000
control_8 0.9225 0.8909 0.9553 0.0000
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 2.9678 0.9953 0.0844
IMD 4.2823 3.9667 0.3644
cat 4.5769 1.9835 0.1000
control_8 5.2699 0.9948 0.0215
GLOBAL 15.8660 12.7158 0.2390 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/control-8-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/control-8-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "control_8",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "IBD Control-8",
  plot_base_path = "plots/uc/hard-flare/ibd/control-8",
  break_time_by = 500,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for continuous control_8 variable
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + IMD + cat + control_8 + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

uc.hard.forest <- rbind(uc.hard.forest, get_HR(fit.me, "control_8"))

# Display plot and model summary
knitr::include_graphics("plots/uc/hard-flare/ibd/control-8.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/hard-flare/ibd/control-8"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3913 1.0274 1.8840 0.0328
IMD2 1.3477 0.6476 2.8046 0.4248
IMD3 1.5031 0.7484 3.0187 0.2520
IMD4 2.1278 1.0942 4.1375 0.0261
IMD5 1.6178 0.8327 3.1432 0.1557
catFC 50-250 2.0148 1.4131 2.8728 0.0001
catFC > 250 3.0295 2.0721 4.4293 0.0000
control_8 0.9453 0.9026 0.9899 0.0169
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1400 0.9860 0.7027
IMD 2.7617 3.9396 0.5890
cat 3.5264 1.9665 0.1671
control_8 1.7488 0.9845 0.1826
GLOBAL 8.2255 23.1839 0.9982 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/control-8-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/control-8-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

IBD Control visual analogue scale

Crohn’s disease

Code
# Handle VAS control (already categorized as character)
flare.cd.df$vas_control_cat <- factor(flare.cd.df$vas_control)

# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "vas_control",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "IBD Control VAS",
  plot_base_path = "plots/cd/soft-flare/ibd/control-vas",
  break_time_by = 200,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for vas_control variable
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + IMD + cat + vas_control + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/control-vas.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/control-vas"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.0996 1.6160 2.7280 0.0000
IMD2 0.7695 0.4706 1.2582 0.2962
IMD3 0.9004 0.5525 1.4672 0.6736
IMD4 0.8423 0.5237 1.3546 0.4789
IMD5 0.9688 0.6121 1.5333 0.8924
catFC 50-250 1.3678 1.0415 1.7963 0.0243
catFC > 250 2.3897 1.7645 3.2362 0.0000
vas_control85+ 0.5710 0.4500 0.7246 0.0000
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0298 0.9934 0.8610
IMD 5.8028 3.9533 0.2095
cat 1.3945 1.9837 0.4940
vas_control 6.4668 0.9893 0.0108
GLOBAL 13.2734 12.1071 0.3576 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/control-vas-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/control-vas-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "vas_control",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "IBD Control VAS",
  plot_base_path = "plots/cd/hard-flare/ibd/control-vas",
  break_time_by = 500,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for continuous vas_control variable
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + IMD + cat + vas_control + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/control-vas.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/control-vas"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3751 1.0101 1.8721 0.0430
IMD2 0.6859 0.3715 1.2665 0.2282
IMD3 0.8599 0.4699 1.5734 0.6244
IMD4 0.7481 0.4148 1.3491 0.3347
IMD5 0.7802 0.4412 1.3800 0.3937
catFC 50-250 1.8928 1.3283 2.6973 0.0004
catFC > 250 3.4895 2.3810 5.1142 0.0000
vas_control85+ 0.9359 0.6894 1.2705 0.6708
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 1.6756 0.9837 0.1918
IMD 7.9021 3.9395 0.0918
cat 8.1082 1.9864 0.0171
vas_control 0.2524 0.9877 0.6102
GLOBAL 20.3223 15.8147 0.1970 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/control-vas-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/control-vas-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Ulcerative colitis

Code
# Handle VAS control (already categorized as character)
flare.uc.df$vas_control_cat <- factor(flare.uc.df$vas_control)

# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "vas_control",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "IBD Control VAS",
  plot_base_path = "plots/uc/soft-flare/ibd/control-vas",
  break_time_by = 200,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for vas_control variable
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + IMD + cat + vas_control + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

# Display plot and model summary
knitr::include_graphics("plots/uc/soft-flare/ibd/control-vas.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/soft-flare/ibd/control-vas"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.5419 1.2247 1.9413 0.0002
IMD2 1.3325 0.7893 2.2496 0.2826
IMD3 1.2056 0.7280 1.9968 0.4675
IMD4 1.4696 0.9046 2.3874 0.1199
IMD5 1.3493 0.8371 2.1749 0.2188
catFC 50-250 1.5798 1.2174 2.0502 0.0006
catFC > 250 1.9587 1.4652 2.6185 0.0000
vas_control85+ 0.6861 0.5423 0.8679 0.0017
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 2.9736 0.9943 0.0840
IMD 4.3116 3.9618 0.3601
cat 5.7100 1.9817 0.0566
vas_control 3.5219 0.9923 0.0599
GLOBAL 16.5023 13.4488 0.2495 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/control-vas-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/control-vas-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "vas_control",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "IBD Control VAS",
  plot_base_path = "plots/uc/hard-flare/ibd/control-vas",
  break_time_by = 200,
  palette = c("#1A8FE3", "#E76D83")
)

# Extract hazard ratio for vas_control variable
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + IMD + cat + vas_control + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

# Display plot and model summary
knitr::include_graphics("plots/uc/hard-flare/ibd/control-vas.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/hard-flare/ibd/control-vas"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3848 1.0220 1.8763 0.0357
IMD2 1.2918 0.6208 2.6877 0.4935
IMD3 1.4220 0.7087 2.8532 0.3218
IMD4 1.9936 1.0297 3.8599 0.0407
IMD5 1.5530 0.8031 3.0034 0.1908
catFC 50-250 2.0459 1.4357 2.9155 0.0001
catFC > 250 2.9077 1.9851 4.2592 0.0000
vas_control85+ 0.6559 0.4815 0.8933 0.0075
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1785 0.9852 0.6666
IMD 2.8643 3.9383 0.5711
cat 3.3111 1.9665 0.1862
vas_control 0.2632 0.9852 0.6017
GLOBAL 7.3842 23.0999 0.9992 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/control-vas-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/control-vas-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Biologic use

Crohn’s disease

Code
# Transform biologic variable and create categorized version
flare.cd.df <- flare.cd.df %>%
  mutate(Biologic = plyr::mapvalues(Biologic,
                                    from = c("Current",
                                             "Previously",
                                             "Never prescribed"),
                                    to = c("Prescribed",
                                           "Not prescribed",
                                           "Not prescribed")))

# Create categorized version for survival analysis
flare.cd.df$Biologic_cat <- flare.cd.df$Biologic

# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Biologic",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Biologic",
  plot_base_path = "plots/cd/soft-flare/ibd/biologic",
  break_time_by = 200,
  palette = c("#1A8FE3", "#E76D83")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + IMD + cat + Biologic + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/biologic.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/biologic"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.9618 1.5445 2.4920 0.0000
IMD2 0.9342 0.5961 1.4640 0.7665
IMD3 0.8935 0.5657 1.4112 0.6293
IMD4 0.9413 0.6065 1.4608 0.7872
IMD5 0.9803 0.6415 1.4982 0.9269
catFC 50-250 1.5753 1.2209 2.0325 0.0005
catFC > 250 2.4212 1.8260 3.2105 0.0000
BiologicNot prescribed 1.1577 0.9238 1.4509 0.2034
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.3184 0.9938 0.5700
IMD 5.8410 3.9587 0.2071
cat 2.2735 1.9848 0.3178
Biologic 6.5065 0.9707 0.0102
GLOBAL 15.8608 13.0633 0.2607 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/biologic-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/biologic-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Biologic",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Biologic",
  plot_base_path = "plots/cd/hard-flare/ibd/biologic",
  break_time_by = 500,
  palette = c("#1A8FE3", "#E76D83")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + IMD + cat + Biologic + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/biologic.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/biologic"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3883 1.0557 1.8257 0.0189
IMD2 0.9219 0.5363 1.5847 0.7685
IMD3 0.9676 0.5566 1.6820 0.9070
IMD4 0.8948 0.5220 1.5340 0.6862
IMD5 0.9034 0.5369 1.5200 0.7019
catFC 50-250 2.0215 1.4724 2.7754 0.0000
catFC > 250 3.3369 2.3692 4.6999 0.0000
BiologicNot prescribed 1.0022 0.7633 1.3158 0.9876
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2572 0.9867 0.6065
IMD 4.2175 3.9404 0.3689
cat 8.8719 1.9846 0.0116
Biologic 0.0442 0.9583 0.8196
GLOBAL 14.1651 20.5809 0.8464 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/biologic-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/biologic-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Ulcerative colitis

Code
# Transform biologic variable and create categorized version
flare.uc.df <- flare.uc.df %>%
  mutate(Biologic = plyr::mapvalues(Biologic,
                                    from = c("Current",
                                             "Previously",
                                             "Never prescribed"),
                                    to = c("Prescribed",
                                           "Not prescribed",
                                           "Not prescribed")))

# Create categorized version for survival analysis
flare.uc.df$Biologic_cat <- flare.uc.df$Biologic

# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "Biologic",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Biologic",
  plot_base_path = "plots/uc/soft-flare/ibd/biologic",
  break_time_by = 200,
  palette = c("#1A8FE3", "#E76D83")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + IMD + cat + Biologic + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

# Display plot and model summary
knitr::include_graphics("plots/uc/soft-flare/ibd/biologic.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/soft-flare/ibd/biologic"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.5457 1.2492 1.9126 0.0001
IMD2 1.2605 0.7970 1.9934 0.3223
IMD3 1.0597 0.6766 1.6596 0.8001
IMD4 1.4046 0.9150 2.1560 0.1202
IMD5 1.1438 0.7503 1.7438 0.5322
catFC 50-250 1.5415 1.2061 1.9702 0.0005
catFC > 250 2.1682 1.6622 2.8283 0.0000
BiologicNot prescribed 1.6853 1.2700 2.2363 0.0003
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 1.2814 0.9920 0.2554
IMD 4.0168 3.9495 0.3963
cat 5.5881 1.9739 0.0597
Biologic 0.8627 0.9374 0.3312
GLOBAL 12.2828 17.5375 0.8105 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/biologic-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/biologic-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "Biologic",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Biologic",
  plot_base_path = "plots/uc/hard-flare/ibd/biologic",
  break_time_by = 500,
  palette = c("#1A8FE3", "#E76D83")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + IMD + cat + Biologic + frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

# Display plot and model summary
knitr::include_graphics("plots/uc/hard-flare/ibd/biologic.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/hard-flare/ibd/biologic"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3299 1.0239 1.7273 0.0326
IMD2 1.4037 0.7829 2.5167 0.2549
IMD3 1.3959 0.7936 2.4554 0.2471
IMD4 1.7610 1.0200 3.0403 0.0422
IMD5 1.3247 0.7702 2.2783 0.3095
catFC 50-250 2.0523 1.5024 2.8035 0.0000
catFC > 250 3.2198 2.3238 4.4613 0.0000
BiologicNot prescribed 0.8410 0.6128 1.1542 0.2838
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1497 0.9849 0.6927
IMD 2.5933 3.9365 0.6182
cat 4.3484 1.9676 0.1106
Biologic 3.8871 0.9212 0.0431
GLOBAL 11.4032 24.4735 0.9884 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/biologic-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/biologic-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Variables only relevant to Crohn’s disease

Code
demo.cd <- readRDS(paste0(paths$outdir, "demo-cd.RDS"))
demo.cd <- demo.cd[, c(
  "ParticipantNo",
  "Surgery",
  "Location",
  "L4",
  "Behaviour",
  "Perianal",
  "HBI",
  "PRO2"
)]
flare.cd.df <- merge(flare.cd.df,
  demo.cd,
  by = "ParticipantNo",
  all.x = TRUE,
  all.y = FALSE
)

# Create categorized versions for survival analysis
flare.cd.df$Surgery_cat <- flare.cd.df$Surgery
flare.cd.df$Location_cat <- flare.cd.df$Location
flare.cd.df$L4_cat <- flare.cd.df$L4
flare.cd.df$Behaviour_cat <- flare.cd.df$Behaviour
flare.cd.df$Perianal_cat <- flare.cd.df$Perianal

Surgery

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Surgery",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Previous surgery",
  plot_base_path = "plots/cd/soft-flare/ibd/surgery",
  break_time_by = 200,
  palette = c("#3DDC97", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + cat + IMD + Surgery +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)
cd.clin.forest <- rbind(cd.clin.forest, get_HR(fit.me, "SurgeryYes"))

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/surgery.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/surgery"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.9970 1.5704 2.5394 0.0000
catFC 50-250 1.5774 1.2210 2.0379 0.0005
catFC > 250 2.3834 1.7891 3.1751 0.0000
IMD2 0.8921 0.5692 1.3983 0.6186
IMD3 0.8326 0.5264 1.3170 0.4336
IMD4 0.8750 0.5640 1.3575 0.5512
IMD5 0.9325 0.6106 1.4240 0.7462
SurgeryYes 1.0130 0.8111 1.2652 0.9092
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.3138 0.9943 0.5730
cat 2.4348 1.9874 0.2936
IMD 5.3005 3.9665 0.2539
Surgery 0.6320 0.9915 0.4234
GLOBAL 8.9435 11.8427 0.6958 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/surgery-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/surgery-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Surgery",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Previous surgery",
  plot_base_path = "plots/cd/hard-flare/ibd/surgery",
  break_time_by = 200,
  palette = c("#3DDC97", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + cat + IMD + Surgery +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)
cd.hard.forest <- rbind(cd.hard.forest, get_HR(fit.me, "SurgeryYes"))

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/surgery.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/surgery"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4133 1.0743 1.8592 0.0134
catFC 50-250 2.0584 1.4977 2.8292 0.0000
catFC > 250 3.2393 2.2905 4.5812 0.0000
IMD2 0.8768 0.5106 1.5058 0.6338
IMD3 0.8923 0.5114 1.5569 0.6882
IMD4 0.8420 0.4913 1.4430 0.5314
IMD5 0.8745 0.5201 1.4706 0.6132
SurgeryYes 0.7516 0.5744 0.9835 0.0374
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.3300 0.9858 0.5597
cat 8.5496 1.9843 0.0137
IMD 4.0006 3.9415 0.3973
Surgery 2.5667 0.9872 0.1073
GLOBAL 17.3483 20.3096 0.6492 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/surgery-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/surgery-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Montreal location

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Location",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Montreal location",
  plot_base_path = "plots/cd/soft-flare/ibd/location",
  break_time_by = 200,
  palette = c("#3B3561", "#5BC0EB", "#FFD166", "#E4572E")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + cat + IMD + Location +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.clin.forest <- rbind(
  cd.clin.forest,
  get_HR(fit.me, c("LocationL2", "LocationL3"))
)

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/location.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/location"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.0379 1.5673 2.6498 0.0000
catFC 50-250 1.6798 1.2692 2.2231 0.0003
catFC > 250 2.4535 1.8007 3.3429 0.0000
IMD2 0.9218 0.5727 1.4836 0.7372
IMD3 0.8187 0.5015 1.3365 0.4237
IMD4 0.8363 0.5225 1.3386 0.4564
IMD5 0.9443 0.6018 1.4816 0.8029
LocationL2 0.8397 0.6095 1.1570 0.2855
LocationL3 1.0929 0.8212 1.4543 0.5425
LocationL4 only 0.9482 0.2984 3.0126 0.9281
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2323 0.9980 0.6290
cat 3.2675 1.9968 0.1947
IMD 5.1401 3.9898 0.2720
Location 1.1604 2.9932 0.7614
GLOBAL 9.8488 10.7869 0.5251 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/location-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/location-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Location",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Montreal location",
  plot_base_path = "plots/cd/hard-flare/ibd/location",
  break_time_by = 200,
  palette = c("#3B3561", "#5BC0EB", "#FFD166", "#E4572E")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + cat + IMD + Location +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.hard.forest <- rbind(
  cd.hard.forest,
  get_HR(fit.me, c("LocationL2", "LocationL3"))
)

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/location.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/location"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4122 1.0506 1.8982 0.0222
catFC 50-250 1.9337 1.3740 2.7214 0.0002
catFC > 250 3.1807 2.2094 4.5789 0.0000
IMD2 0.8819 0.4941 1.5742 0.6707
IMD3 0.9504 0.5256 1.7185 0.8663
IMD4 0.9099 0.5129 1.6144 0.7470
IMD5 0.9891 0.5687 1.7202 0.9691
LocationL2 1.2964 0.8834 1.9025 0.1846
LocationL3 1.3328 0.9285 1.9132 0.1193
LocationL4 only 1.2113 0.2871 5.1108 0.7941
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1424 0.9830 0.6991
cat 6.3620 1.9837 0.0409
IMD 3.8329 3.9412 0.4203
Location 0.5903 2.9550 0.8941
GLOBAL 11.3705 21.5180 0.9629 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/location-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/location-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Montreal L4

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "L4",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Montreal L4 presence",
  plot_base_path = "plots/cd/soft-flare/ibd/l4",
  break_time_by = 200,
  palette = c("#3DDC97", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + cat + IMD + L4 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)
cd.clin.forest <- rbind(
  cd.clin.forest,
  get_HR(fit.me, "L4Present")
)

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/l4.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/l4"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.0268 1.5589 2.6350 0.0000
catFC 50-250 1.6892 1.2767 2.2352 0.0002
catFC > 250 2.3738 1.7423 3.2342 0.0000
IMD2 0.8815 0.5472 1.4199 0.6040
IMD3 0.7833 0.4794 1.2799 0.3296
IMD4 0.8069 0.5031 1.2941 0.3733
IMD5 0.9046 0.5765 1.4193 0.6625
L4Present 1.2936 0.9172 1.8244 0.1423
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2367 0.9981 0.6258
cat 3.4174 1.9965 0.1806
IMD 5.2745 3.9899 0.2591
L4 2.9817 0.9968 0.0838
GLOBAL 12.5997 8.8065 0.1700 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/l4-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/l4-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "L4",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Montreal L4 presence",
  plot_base_path = "plots/cd/hard-flare/ibd/l4",
  break_time_by = 200,
  palette = c("#3DDC97", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + cat + IMD + L4 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.hard.forest <- rbind(
  cd.hard.forest,
  get_HR(fit.me, "L4Present")
)

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/l4.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/l4"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4360 1.0684 1.9301 0.0165
catFC 50-250 1.9358 1.3756 2.7240 0.0002
catFC > 250 3.1320 2.1758 4.5085 0.0000
IMD2 0.8257 0.4627 1.4736 0.5170
IMD3 0.8677 0.4791 1.5715 0.6396
IMD4 0.8321 0.4680 1.4794 0.5314
IMD5 0.9026 0.5188 1.5705 0.7169
L4Present 1.5676 1.0693 2.2981 0.0213
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1595 0.9826 0.6825
cat 6.4740 1.9836 0.0386
IMD 4.0148 3.9425 0.3956
L4 0.0375 0.9606 0.8338
GLOBAL 10.7469 19.2888 0.9386 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/l4-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/l4-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Montreal Behaviour

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Behaviour",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Montreal behaviour",
  plot_base_path = "plots/cd/soft-flare/ibd/behaviour",
  break_time_by = 200,
  palette = c("#3DDC97", "#FFD166", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + cat + IMD + Behaviour +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.clin.forest <- rbind(
  cd.clin.forest,
  get_HR(fit.me, c("BehaviourB2", "BehaviourB3"))
)

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/behaviour.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/behaviour"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.0640 1.5766 2.7020 0.0000
catFC 50-250 1.6340 1.2279 2.1746 0.0008
catFC > 250 2.3258 1.6953 3.1909 0.0000
IMD2 0.9230 0.5693 1.4965 0.7452
IMD3 0.7493 0.4518 1.2430 0.2637
IMD4 0.8865 0.5480 1.4342 0.6236
IMD5 0.9870 0.6239 1.5614 0.9554
BehaviourB2 1.0662 0.7953 1.4294 0.6683
BehaviourB3 0.7879 0.4947 1.2549 0.3155
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1946 0.9978 0.6582
cat 3.9115 1.9961 0.1410
IMD 4.3369 3.9877 0.3606
Behaviour 0.9019 1.9953 0.6359
GLOBAL 9.4596 10.0044 0.4895 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/behaviour-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/behaviour-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Behaviour",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Montreal behaviour",
  plot_base_path = "plots/cd/hard-flare/ibd/behaviour",
  break_time_by = 200,
  palette = c("#3DDC97", "#FFD166", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + cat + IMD + Behaviour +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.hard.forest <- rbind(
  cd.hard.forest,
  get_HR(fit.me, c("BehaviourB2", "BehaviourB3"))
)

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/behaviour.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/behaviour"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.5091 1.1174 2.0383 0.0073
catFC 50-250 2.0811 1.4681 2.9499 0.0000
catFC > 250 3.4647 2.3857 5.0316 0.0000
IMD2 0.9419 0.5228 1.6971 0.8421
IMD3 0.8970 0.4872 1.6518 0.7272
IMD4 0.9794 0.5448 1.7605 0.9444
IMD5 1.0317 0.5866 1.8145 0.9137
BehaviourB2 0.8913 0.6267 1.2677 0.5221
BehaviourB3 0.8505 0.5012 1.4431 0.5483
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0070 0.9845 0.9303
cat 5.6993 1.9835 0.0570
IMD 2.8336 3.9427 0.5771
Behaviour 2.7982 1.9761 0.2427
GLOBAL 12.1681 20.1919 0.9154 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/behaviour-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/behaviour-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Perinanal disease

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Perianal",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Perianal disease",
  plot_base_path = "plots/cd/soft-flare/ibd/perinal",
  break_time_by = 200,
  palette = c("#3DDC97", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + cat + IMD + Perianal +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.clin.forest <- rbind(
  cd.clin.forest,
  get_HR(fit.me, "PerianalYes")
)

# Display plot and model summary
knitr::include_graphics("plots/cd/soft-flare/ibd/perinal.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/soft-flare/ibd/perinal"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.0572 1.5795 2.6794 0.0000
catFC 50-250 1.6277 1.2282 2.1573 0.0007
catFC > 250 2.3910 1.7543 3.2587 0.0000
IMD2 0.9586 0.5934 1.5485 0.8628
IMD3 0.8177 0.4975 1.3441 0.4274
IMD4 0.8333 0.5158 1.3464 0.4564
IMD5 0.8992 0.5679 1.4237 0.6503
PerianalYes 1.3022 1.0174 1.6667 0.0360
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0684 0.9950 0.7920
cat 4.5984 1.9892 0.0994
IMD 5.5262 3.9719 0.2343
Perianal 0.2214 0.9919 0.6346
GLOBAL 10.2702 10.9257 0.4997 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/perinal-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/soft-flare/ibd/perinal-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.cd.df,
  var_name = "Perianal",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Perianal disease",
  plot_base_path = "plots/cd/hard-flare/ibd/perinal",
  break_time_by = 200,
  palette = c("#3DDC97", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + cat + IMD + Perianal +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

cd.hard.forest <- rbind(
  cd.hard.forest,
  get_HR(fit.me, "PerianalYes")
)

# Display plot and model summary
knitr::include_graphics("plots/cd/hard-flare/ibd/perinal.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/cd/hard-flare/ibd/perinal"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4036 1.0483 1.8795 0.0228
catFC 50-250 1.8951 1.3518 2.6569 0.0002
catFC > 250 3.1094 2.1650 4.4659 0.0000
IMD2 0.9969 0.5566 1.7855 0.9918
IMD3 1.0334 0.5690 1.8768 0.9140
IMD4 0.9261 0.5160 1.6621 0.7969
IMD5 0.9379 0.5339 1.6477 0.8235
PerianalYes 1.1428 0.8537 1.5299 0.3697
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1937 0.9841 0.6533
cat 7.0942 1.9838 0.0283
IMD 2.3289 3.9479 0.6677
Perianal 1.9313 0.9841 0.1614
GLOBAL 11.5538 19.8863 0.9279 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/perinal-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/cd/hard-flare/ibd/perinal-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Variables only relevant to UC/IBDU

Code
demo.uc <- readRDS(paste0(paths$outdir, "demo-uc.RDS"))
demo.uc <- demo.uc[, c("ParticipantNo", "Extent", "Mayo")]
flare.uc.df <- merge(flare.uc.df,
  demo.uc,
  by = "ParticipantNo",
  all.x = TRUE,
  all.y = FALSE
)

# Create categorized version for survival analysis
flare.uc.df$Extent_cat <- flare.uc.df$Extent

Montreal extent

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "Extent",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Montreal extent",
  plot_base_path = "plots/uc/soft-flare/ibd/extent",
  break_time_by = 200,
  palette = c("#3DDC97", "orange", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + cat + IMD + Extent +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

uc.clin.forest <- rbind(
  uc.clin.forest,
  get_HR(fit.me, c("ExtentE2", "ExtentE3"))
)

# Display plot and model summary
knitr::include_graphics("plots/uc/soft-flare/ibd/extent.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/soft-flare/ibd/extent"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4261 1.1158 1.8227 0.0046
catFC 50-250 1.6346 1.2328 2.1673 0.0006
catFC > 250 2.2553 1.6710 3.0440 0.0000
IMD2 1.2814 0.7812 2.1019 0.3262
IMD3 1.0697 0.6587 1.7371 0.7854
IMD4 1.4373 0.9127 2.2636 0.1175
IMD5 1.1749 0.7483 1.8449 0.4838
ExtentE2 0.9196 0.6607 1.2801 0.6196
ExtentE3 0.6967 0.4810 1.0090 0.0558
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0215 0.9891 0.8804
cat 4.6159 1.9587 0.0959
IMD 2.9011 3.9332 0.5640
Extent 4.1788 1.9657 0.1202
GLOBAL 12.0732 21.9250 0.9547 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/extent-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/extent-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "Extent",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Montreal extent",
  plot_base_path = "plots/uc/hard-flare/ibd/extent",
  break_time_by = 200,
  palette = c("#3DDC97", "orange", "#DD7373")
)

# Cox model
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + cat + IMD + Extent +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

uc.hard.forest <- rbind(
  uc.hard.forest,
  get_HR(fit.me, c("ExtentE2", "ExtentE3"))
)

# Display plot and model summary
knitr::include_graphics("plots/uc/hard-flare/ibd/extent.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/hard-flare/ibd/extent"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4242 1.0492 1.9334 0.0234
catFC 50-250 1.9574 1.3666 2.8037 0.0002
catFC > 250 3.2458 2.2460 4.6907 0.0000
IMD2 1.6852 0.8594 3.3047 0.1288
IMD3 1.6813 0.8824 3.2033 0.1142
IMD4 2.2154 1.1918 4.1182 0.0119
IMD5 1.6230 0.8722 3.0201 0.1264
ExtentE2 2.1493 1.2632 3.6570 0.0048
ExtentE3 1.8507 1.0529 3.2532 0.0324
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.6040 0.9857 0.4315
cat 5.9334 1.9616 0.0496
IMD 2.5189 3.9364 0.6315
Extent 0.7654 1.9746 0.6760
GLOBAL 9.5946 21.8882 0.9891 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/extent-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/extent-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Mayo score

Code
flare.uc.df$Mayo_cat <- cut(flare.uc.df$Mayo,
  c(0, 2, 5, 9),
  labels = c("Inactive", "Mild", "Moderate"),
  right = TRUE
)
Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "Mayo",
  outcome_time = "softflare_time",
  outcome_event = "softflare",
  legend_title = "Mayo score",
  plot_base_path = "plots/uc/soft-flare/ibd/mayo",
  break_time_by = 200,
  palette = c("#3DDC97", "orange", "#DD7373")
)

# Cox model using continuous Mayo variable
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex + cat + Mayo + Extent +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

uc.clin.forest <- rbind(
  uc.clin.forest,
  get_HR(fit.me, "Mayo")
)

# Display plot and model summary
knitr::include_graphics("plots/uc/soft-flare/ibd/mayo.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/soft-flare/ibd/mayo"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4049 1.0761 1.8343 0.0125
catFC 50-250 1.6646 1.2283 2.2557 0.0010
catFC > 250 2.4320 1.7451 3.3892 0.0000
Mayo 1.1117 1.0094 1.2243 0.0316
ExtentE2 1.0050 0.6957 1.4518 0.9789
ExtentE3 0.7477 0.4960 1.1274 0.1652
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2221 0.9943 0.6351
cat 3.2655 1.9791 0.1924
Mayo 2.8510 0.9265 0.0823
Extent 4.4854 1.9674 0.1032
GLOBAL 10.5608 13.9227 0.7149 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/mayo-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/soft-flare/ibd/mayo-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
# Run survival analysis using utility function
analysis_result <- run_survival_analysis(
  data = flare.uc.df,
  var_name = "Mayo",
  outcome_time = "hardflare_time",
  outcome_event = "hardflare",
  legend_title = "Mayo score",
  plot_base_path = "plots/uc/hard-flare/ibd/mayo",
  break_time_by = 200,
  palette = c("#3DDC97", "orange", "#DD7373")
)

# Cox model using continuous Mayo variable
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex + cat + IMD + Mayo +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

uc.hard.forest <- rbind(
  uc.hard.forest,
  get_HR(fit.me, "Mayo")
)

# Display plot and model summary
knitr::include_graphics("plots/uc/hard-flare/ibd/mayo.png")

Code
invisible(cox_summary(fit.me, plot_base_path = "plots/uc/hard-flare/ibd/mayo"))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.2672 0.9438 1.7015 0.1152
catFC 50-250 1.8912 1.3331 2.6830 0.0004
catFC > 250 3.1128 2.1355 4.5373 0.0000
IMD2 1.6242 0.8553 3.0844 0.1383
IMD3 1.1995 0.6373 2.2574 0.5729
IMD4 1.5787 0.8645 2.8828 0.1373
IMD5 1.2277 0.6754 2.2314 0.5011
Mayo 1.0890 0.9753 1.2161 0.1297
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0043 0.9878 0.9457
cat 3.1681 1.9834 0.2026
IMD 1.4273 3.9518 0.8342
Mayo 0.1658 0.9504 0.6632
GLOBAL 4.6901 17.1518 0.9987 
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/mayo-dfbeta.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

[1] “plots/uc/hard-flare/ibd/mayo-martingale.png” attr(,“class”) [1] “knit_image_paths” “knit_asis”

Code
saveRDS(flare.df, paste0(paths$outdir, "flares-IBD.RDS"))
saveRDS(flare.cd.df, paste0(paths$outdir, "flares-IBD-cd.RDS"))
saveRDS(flare.uc.df, paste0(paths$outdir, "flares-IBD-uc.RDS"))

saveRDS(cd.clin.forest, paste0(paths$outdir, "cd-clin-IBD.RDS"))
saveRDS(cd.hard.forest, paste0(paths$outdir, "cd-hard-IBD.RDS"))
saveRDS(uc.clin.forest, paste0(paths$outdir, "uc-clin-IBD.RDS"))
saveRDS(uc.hard.forest, paste0(paths$outdir, "uc-hard-IBD.RDS"))
Control for additional covariates

Crohn’s disease

Patient-reported flare

IBD Control-8
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    control_8 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.0592 1.5443 2.7457 0.0000
catFC 50-250 1.3701 1.0310 1.8208 0.0300
catFC > 250 2.2602 1.6428 3.1097 0.0000
IMD2 0.7494 0.4529 1.2397 0.2613
IMD3 0.7689 0.4605 1.2837 0.3149
IMD4 0.8471 0.5214 1.3762 0.5027
IMD5 0.9123 0.5686 1.4637 0.7036
IBD Duration 0.9874 0.9760 0.9991 0.0345
BMI 0.9987 0.9761 1.0218 0.9109
TreatmentMono biologic 0.9608 0.6604 1.3978 0.8344
TreatmentCombo therapy 0.7384 0.4587 1.1886 0.2118
Treatment5-ASA 1.4528 0.8007 2.6359 0.2192
TreatmentNone reported 1.0019 0.7064 1.4210 0.9916
Age 1.0112 1.0018 1.0207 0.0195
control_8 0.8928 0.8609 0.9258 0.0000
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0122 0.9924 0.9104
cat 0.7938 1.9833 0.6684
IMD 6.4286 3.9480 0.1648
IBD Duration 2.0041 0.9953 0.1560
BMI 1.5485 0.9898 0.2109
Treatment 5.9298 3.9194 0.1964
Age 0.5089 0.9905 0.4718
control_8 8.9564 0.9790 0.0027
GLOBAL 29.0747 18.3377 0.0529 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

IBD Control VAS
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    vas_control +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.2982 1.7321 3.0491 0.0000
catFC 50-250 1.3127 0.9852 1.7489 0.0631
catFC > 250 2.2600 1.6390 3.1162 0.0000
IMD2 0.7181 0.4360 1.1829 0.1934
IMD3 0.8002 0.4821 1.3283 0.3887
IMD4 0.8597 0.5319 1.3895 0.5373
IMD5 0.9224 0.5802 1.4664 0.7327
IBD Duration 0.9885 0.9768 1.0004 0.0577
BMI 1.0023 0.9797 1.0255 0.8405
TreatmentMono biologic 1.0721 0.7395 1.5544 0.7132
TreatmentCombo therapy 0.7799 0.4860 1.2517 0.3031
Treatment5-ASA 1.5093 0.8385 2.7170 0.1699
TreatmentNone reported 0.9824 0.6911 1.3965 0.9211
Age 1.0074 0.9981 1.0169 0.1193
vas_control85+ 0.6000 0.4663 0.7720 0.0001
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0028 0.9989 0.9574
cat 0.7396 1.9974 0.6903
IMD 5.9711 3.9910 0.2004
IBD Duration 2.3800 0.9991 0.1228
BMI 1.0571 0.9985 0.3034
Treatment 6.4094 3.9868 0.1694
Age 0.5997 0.9982 0.4380
vas_control 5.6480 0.9978 0.0174
GLOBAL 23.0400 15.4505 0.0958 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Previous surgery
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Surgery +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.1614 1.6662 2.8038 0.0000
catFC 50-250 1.5087 1.1533 1.9736 0.0027
catFC > 250 2.1917 1.6122 2.9796 0.0000
IMD2 0.8105 0.5092 1.2902 0.3757
IMD3 0.7676 0.4740 1.2429 0.2821
IMD4 0.8583 0.5444 1.3532 0.5106
IMD5 0.8803 0.5681 1.3641 0.5683
IBD Duration 0.9923 0.9804 1.0044 0.2141
BMI 1.0061 0.9842 1.0284 0.5886
TreatmentMono biologic 0.9472 0.6639 1.3513 0.7649
TreatmentCombo therapy 0.6983 0.4532 1.0761 0.1036
Treatment5-ASA 1.3343 0.7494 2.3758 0.3272
TreatmentNone reported 0.9486 0.6820 1.3193 0.7537
Age 1.0052 0.9962 1.0143 0.2572
SurgeryYes 1.1234 0.8754 1.4416 0.3606
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.3388 0.9914 0.5570
cat 1.5116 1.9845 0.4660
IMD 5.5692 3.9560 0.2289
IBD Duration 3.4573 0.9954 0.0626
BMI 2.1444 0.9904 0.1414
Treatment 3.3457 3.9107 0.4878
Age 2.2642 0.9915 0.1310
Surgery 0.9573 0.9917 0.3251
GLOBAL 18.2683 19.7881 0.5561 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal location
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Location +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.1798 1.6498 2.8800 0.0000
catFC 50-250 1.5990 1.1978 2.1346 0.0014
catFC > 250 2.2763 1.6429 3.1539 0.0000
IMD2 0.8669 0.5308 1.4156 0.5680
IMD3 0.8129 0.4875 1.3557 0.4274
IMD4 0.8865 0.5467 1.4376 0.6254
IMD5 0.9033 0.5665 1.4402 0.6691
IBD Duration 0.9912 0.9789 1.0036 0.1628
BMI 1.0056 0.9826 1.0290 0.6378
TreatmentMono biologic 0.9615 0.6648 1.3905 0.8347
TreatmentCombo therapy 0.6733 0.4271 1.0613 0.0884
Treatment5-ASA 1.3604 0.7381 2.5073 0.3238
TreatmentNone reported 0.9027 0.6358 1.2818 0.5673
Age 1.0078 0.9983 1.0175 0.1091
LocationL2 0.8025 0.5737 1.1226 0.1989
LocationL3 1.1264 0.8342 1.5210 0.4371
LocationL4 only 0.9506 0.2968 3.0448 0.9320
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.3013 0.9948 0.5809
cat 2.3312 1.9920 0.3102
IMD 6.2866 3.9762 0.1766
IBD Duration 3.1256 0.9966 0.0767
BMI 1.9769 0.9937 0.1585
Treatment 5.9104 3.9477 0.2007
Age 1.4629 0.9962 0.2255
Location 1.3563 2.9811 0.7126
GLOBAL 23.6618 19.0505 0.2118 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal L4 presence
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    L4 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.1709 1.6432 2.8680 0.0000
catFC 50-250 1.5972 1.1963 2.1323 0.0015
catFC > 250 2.1855 1.5764 3.0300 0.0000
IMD2 0.8134 0.4973 1.3304 0.4106
IMD3 0.7575 0.4543 1.2632 0.2871
IMD4 0.8463 0.5213 1.3739 0.4996
IMD5 0.8452 0.5297 1.3485 0.4804
IBD Duration 0.9917 0.9795 1.0041 0.1871
BMI 1.0060 0.9828 1.0296 0.6158
TreatmentMono biologic 0.9471 0.6547 1.3702 0.7731
TreatmentCombo therapy 0.6733 0.4272 1.0613 0.0884
Treatment5-ASA 1.2960 0.7056 2.3804 0.4032
TreatmentNone reported 0.9088 0.6401 1.2905 0.5931
Age 1.0074 0.9978 1.0170 0.1302
L4Present 1.4018 0.9480 2.0730 0.0906
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2963 0.9946 0.5840
cat 2.4313 1.9900 0.2946
IMD 6.4263 3.9729 0.1671
IBD Duration 3.0585 0.9973 0.0800
BMI 2.0835 0.9940 0.1478
Treatment 5.6189 3.9432 0.2233
Age 1.5146 0.9948 0.2171
L4 4.3385 0.9926 0.0368
GLOBAL 26.8201 17.4663 0.0703 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal behaviour
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Behaviour +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.2528 1.6914 3.0005 0.0000
catFC 50-250 1.5728 1.1721 2.1105 0.0025
catFC > 250 2.1176 1.5163 2.9573 0.0000
IMD2 0.8569 0.5202 1.4116 0.5442
IMD3 0.7345 0.4337 1.2438 0.2509
IMD4 0.9233 0.5614 1.5185 0.7533
IMD5 0.9321 0.5790 1.5007 0.7724
IBD Duration 0.9906 0.9780 1.0034 0.1485
BMI 1.0040 0.9800 1.0285 0.7475
TreatmentMono biologic 0.9358 0.6456 1.3562 0.7258
TreatmentCombo therapy 0.6465 0.4055 1.0307 0.0668
Treatment5-ASA 1.2600 0.6823 2.3270 0.4602
TreatmentNone reported 0.8363 0.5854 1.1949 0.3263
Age 1.0056 0.9958 1.0154 0.2631
BehaviourB2 1.1609 0.8572 1.5722 0.3351
BehaviourB3 0.7813 0.4760 1.2826 0.3292
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2640 0.9936 0.6047
cat 2.5321 1.9892 0.2799
IMD 5.2296 3.9699 0.2610
IBD Duration 2.7026 0.9956 0.0996
BMI 1.3795 0.9924 0.2382
Treatment 2.7420 3.9287 0.5907
Age 2.3574 0.9941 0.1237
Behaviour 2.0665 1.9901 0.3538
GLOBAL 19.7745 18.8200 0.3969 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Perianal disease
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Perianal +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 2.2044 1.6637 2.9209 0.0000
catFC 50-250 1.5646 1.1685 2.0950 0.0027
catFC > 250 2.1640 1.5614 2.9990 0.0000
IMD2 0.9120 0.5569 1.4936 0.7145
IMD3 0.8242 0.4927 1.3788 0.4615
IMD4 0.8757 0.5356 1.4318 0.5967
IMD5 0.8869 0.5536 1.4207 0.6175
IBD Duration 0.9850 0.9727 0.9975 0.0190
BMI 0.9953 0.9714 1.0199 0.7064
TreatmentMono biologic 0.9470 0.6530 1.3734 0.7739
TreatmentCombo therapy 0.6016 0.3784 0.9563 0.0316
Treatment5-ASA 1.4594 0.7768 2.7418 0.2400
TreatmentNone reported 0.9591 0.6741 1.3644 0.8162
Age 1.0074 0.9977 1.0171 0.1374
PerianalYes 1.4876 1.1455 1.9319 0.0029
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1272 0.9958 0.7197
cat 3.3291 1.9921 0.1882
IMD 6.1545 3.9768 0.1857
IBD Duration 3.1836 0.9968 0.0741
BMI 2.4785 0.9941 0.1145
Treatment 2.0844 3.9488 0.7129
Age 0.5700 0.9954 0.4485
Perianal 0.0856 0.9950 0.7679
GLOBAL 18.3080 16.8199 0.3580 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Objective flare

IBD Control-8
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    control_8 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.6541 1.1776 2.3236 0.0037
catFC 50-250 1.9360 1.3264 2.8257 0.0006
catFC > 250 3.6730 2.4449 5.5180 0.0000
IMD2 0.7537 0.3989 1.4240 0.3837
IMD3 0.7606 0.3986 1.4515 0.4066
IMD4 0.8737 0.4728 1.6147 0.6665
IMD5 0.9585 0.5301 1.7330 0.8884
IBD Duration 0.9831 0.9674 0.9991 0.0381
BMI 1.0159 0.9871 1.0455 0.2818
TreatmentMono biologic 0.9843 0.6334 1.5294 0.9438
TreatmentCombo therapy 0.7176 0.4058 1.2688 0.2538
Treatment5-ASA 1.3625 0.5994 3.0975 0.4603
TreatmentNone reported 0.7444 0.4872 1.1373 0.1723
Age 0.9937 0.9822 1.0054 0.2918
control_8 0.9878 0.9402 1.0378 0.6252
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 1.5368 0.9967 0.2143
cat 8.2858 1.9973 0.0158
IMD 3.7556 3.9851 0.4378
IBD Duration 0.1345 0.9986 0.7132
BMI 3.3280 0.9983 0.0679
Treatment 3.6997 3.9751 0.4444
Age 2.6202 0.9980 0.1052
control_8 2.5673 0.9953 0.1084
GLOBAL 24.4616 15.7546 0.0740 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

IBD Control VAS
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    vas_control +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.6595 1.1823 2.3294 0.0034
catFC 50-250 1.8986 1.2964 2.7804 0.0010
catFC > 250 3.5997 2.3917 5.4178 0.0000
IMD2 0.7476 0.3939 1.4191 0.3737
IMD3 0.7574 0.3966 1.4464 0.4000
IMD4 0.8255 0.4440 1.5348 0.5444
IMD5 0.9431 0.5225 1.7022 0.8458
IBD Duration 0.9830 0.9671 0.9992 0.0396
BMI 1.0159 0.9868 1.0459 0.2865
TreatmentMono biologic 1.0043 0.6436 1.5671 0.9850
TreatmentCombo therapy 0.7420 0.4178 1.3180 0.3087
Treatment5-ASA 1.3934 0.6111 3.1769 0.4302
TreatmentNone reported 0.7494 0.4886 1.1494 0.1862
Age 0.9933 0.9816 1.0051 0.2667
vas_control85+ 0.9782 0.7022 1.3625 0.8961
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 1.6320 0.9949 0.2002
cat 9.3063 1.9967 0.0095
IMD 3.6734 3.9797 0.4489
IBD Duration 0.0282 0.9980 0.8662
BMI 3.9058 0.9977 0.0480
Treatment 3.9109 3.9673 0.4133
Age 2.6417 0.9972 0.1037
vas_control 0.4292 0.9969 0.5111
GLOBAL 25.8796 16.0653 0.0570 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Previous surgery
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Surgery +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.7063 1.2621 2.3070 0.0005
catFC 50-250 2.0496 1.4628 2.8717 0.0000
catFC > 250 3.2225 2.2238 4.6697 0.0000
IMD2 0.9673 0.5460 1.7135 0.9092
IMD3 0.8549 0.4699 1.5555 0.6077
IMD4 0.9144 0.5137 1.6277 0.7610
IMD5 1.0694 0.6188 1.8479 0.8101
IBD Duration 0.9879 0.9725 1.0035 0.1287
BMI 1.0160 0.9900 1.0427 0.2310
TreatmentMono biologic 0.9737 0.6504 1.4575 0.8968
TreatmentCombo therapy 0.6387 0.3846 1.0606 0.0832
Treatment5-ASA 1.2006 0.5485 2.6277 0.6474
TreatmentNone reported 0.7371 0.4986 1.0896 0.1261
Age 0.9884 0.9778 0.9992 0.0350
SurgeryYes 0.8828 0.6539 1.1920 0.4159
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.3484 0.9859 0.5492
cat 7.4677 1.9869 0.0236
IMD 2.2132 3.9502 0.6893
IBD Duration 0.1496 0.9966 0.6976
BMI 4.5520 0.9899 0.0324
Treatment 3.4629 3.8998 0.4680
Age 4.4954 0.9917 0.0336
Surgery 3.2762 0.9926 0.0696
GLOBAL 25.5592 21.0557 0.2263 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal location
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Location +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.6164 1.1772 2.2193 0.0030
catFC 50-250 1.9638 1.3758 2.8029 0.0002
catFC > 250 3.1927 2.1725 4.6921 0.0000
IMD2 0.9011 0.4947 1.6413 0.7335
IMD3 0.8958 0.4817 1.6656 0.7279
IMD4 0.9706 0.5345 1.7624 0.9218
IMD5 1.1257 0.6357 1.9934 0.6847
IBD Duration 0.9853 0.9697 1.0011 0.0681
BMI 1.0171 0.9901 1.0448 0.2160
TreatmentMono biologic 0.9848 0.6488 1.4948 0.9428
TreatmentCombo therapy 0.6232 0.3696 1.0509 0.0761
Treatment5-ASA 1.3422 0.6047 2.9791 0.4694
TreatmentNone reported 0.6601 0.4350 1.0016 0.0509
Age 0.9903 0.9791 1.0017 0.0951
LocationL2 1.1934 0.8012 1.7774 0.3844
LocationL3 1.2251 0.8404 1.7860 0.2911
LocationL4 only 1.1608 0.2746 4.9062 0.8393
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2665 0.9830 0.5986
cat 5.0495 1.9873 0.0792
IMD 2.9329 3.9494 0.5612
IBD Duration 0.1756 0.9954 0.6734
BMI 3.1934 0.9905 0.0729
Treatment 3.5121 3.8864 0.4584
Age 6.5989 0.9917 0.0101
Location 0.6136 2.9546 0.8885
GLOBAL 23.0040 23.1768 0.4710 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal L4 presence
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    L4 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.6360 1.1925 2.2445 0.0023
catFC 50-250 1.9695 1.3797 2.8115 0.0002
catFC > 250 3.1553 2.1455 4.6404 0.0000
IMD2 0.8635 0.4736 1.5742 0.6319
IMD3 0.8420 0.4526 1.5666 0.5872
IMD4 0.9260 0.5092 1.6838 0.8010
IMD5 1.0554 0.5953 1.8713 0.8535
IBD Duration 0.9848 0.9693 1.0006 0.0589
BMI 1.0177 0.9907 1.0454 0.2007
TreatmentMono biologic 0.9643 0.6343 1.4661 0.8650
TreatmentCombo therapy 0.6258 0.3711 1.0553 0.0787
Treatment5-ASA 1.3677 0.6196 3.0189 0.4382
TreatmentNone reported 0.6540 0.4315 0.9914 0.0454
Age 0.9915 0.9803 1.0028 0.1419
L4Present 1.3791 0.8798 2.1617 0.1610
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.2820 0.9843 0.5888
cat 5.1952 1.9856 0.0735
IMD 3.1652 3.9484 0.5225
IBD Duration 0.1867 0.9963 0.6641
BMI 3.1317 0.9900 0.0757
Treatment 3.6132 3.8867 0.4435
Age 6.6664 0.9912 0.0097
L4 0.0262 0.9780 0.8650
GLOBAL 23.7261 21.7383 0.3470 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal behaviour
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Behaviour +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.7552 1.2730 2.4202 0.0006
catFC 50-250 2.1619 1.5036 3.1084 0.0000
catFC > 250 3.4054 2.2930 5.0575 0.0000
IMD2 0.9987 0.5421 1.8398 0.9966
IMD3 0.8851 0.4677 1.6749 0.7076
IMD4 1.0399 0.5618 1.9247 0.9009
IMD5 1.1860 0.6610 2.1281 0.5673
IBD Duration 0.9842 0.9683 1.0004 0.0561
BMI 1.0176 0.9900 1.0459 0.2132
TreatmentMono biologic 0.9511 0.6225 1.4531 0.8167
TreatmentCombo therapy 0.6567 0.3881 1.1112 0.1171
Treatment5-ASA 1.3831 0.6206 3.0824 0.4277
TreatmentNone reported 0.6731 0.4416 1.0259 0.0656
Age 0.9909 0.9797 1.0022 0.1143
BehaviourB2 1.1275 0.7814 1.6269 0.5211
BehaviourB3 0.9689 0.5583 1.6815 0.9105
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.1413 0.9845 0.7008
cat 4.7787 1.9867 0.0906
IMD 2.2188 3.9518 0.6885
IBD Duration 0.4967 0.9944 0.4787
BMI 2.9289 0.9912 0.0860
Treatment 3.5840 3.8782 0.4465
Age 4.8157 0.9909 0.0278
Behaviour 3.2280 1.9833 0.1966
GLOBAL 22.9547 21.8830 0.3974 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Perianal disease
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Perianal +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.cd.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.6224 1.1862 2.2191 0.0025
catFC 50-250 1.9451 1.3659 2.7699 0.0002
catFC > 250 2.9845 2.0352 4.3766 0.0000
IMD2 1.0469 0.5737 1.9104 0.8813
IMD3 0.9910 0.5317 1.8473 0.9774
IMD4 0.9513 0.5170 1.7503 0.8724
IMD5 1.0728 0.6018 1.9124 0.8118
IBD Duration 0.9822 0.9666 0.9980 0.0270
BMI 1.0124 0.9847 1.0408 0.3847
TreatmentMono biologic 0.9765 0.6412 1.4872 0.9119
TreatmentCombo therapy 0.6243 0.3687 1.0572 0.0796
Treatment5-ASA 1.8879 0.8480 4.2030 0.1196
TreatmentNone reported 0.7692 0.5092 1.1619 0.2124
Age 0.9893 0.9780 1.0006 0.0643
PerianalYes 1.2412 0.9115 1.6902 0.1702
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.5730 0.9848 0.4432
cat 5.3420 1.9880 0.0684
IMD 1.7887 3.9568 0.7689
IBD Duration 0.2294 0.9950 0.6299
BMI 3.8347 0.9898 0.0494
Treatment 4.4126 3.8877 0.3375
Age 7.7106 0.9904 0.0054
Perianal 2.9925 0.9900 0.0825
GLOBAL 23.5249 20.5111 0.2903 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Ulcerative colitis

Patient-reported flare

IBD Control-8
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    control_8 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4705 1.1513 1.8782 0.0020
catFC 50-250 1.6275 1.2386 2.1383 0.0005
catFC > 250 2.0514 1.5185 2.7713 0.0000
IMD2 1.3719 0.8068 2.3328 0.2430
IMD3 1.0876 0.6552 1.8052 0.7453
IMD4 1.4622 0.9019 2.3705 0.1232
IMD5 1.2493 0.7723 2.0211 0.3644
IBD Duration 0.9970 0.9837 1.0105 0.6660
BMI 0.9897 0.9662 1.0138 0.3995
TreatmentMono biologic 0.8073 0.5129 1.2706 0.3548
TreatmentCombo therapy 0.4098 0.1995 0.8416 0.0151
Treatment5-ASA 1.3270 0.9419 1.8697 0.1057
TreatmentNone reported 1.0908 0.7708 1.5437 0.6236
Age 0.9898 0.9807 0.9989 0.0281
control_8 0.9224 0.8890 0.9572 0.0000
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 3.1822 0.9990 0.0743
cat 3.1246 1.9961 0.2091
IMD 3.3548 3.9906 0.4988
IBD Duration 1.4940 0.9985 0.2212
BMI 0.2071 0.9981 0.6483
Treatment 1.7696 3.9821 0.7757
Age 0.0124 0.9954 0.9103
control_8 5.0310 0.9988 0.0248
GLOBAL 17.2676 15.8783 0.3604 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

IBD Control VAS
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    vas_control +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.5394 1.2059 1.9652 0.0005
catFC 50-250 1.6653 1.2666 2.1897 0.0003
catFC > 250 2.0204 1.4942 2.7318 0.0000
IMD2 1.3921 0.8107 2.3903 0.2304
IMD3 1.1061 0.6603 1.8529 0.7016
IMD4 1.4684 0.8974 2.4028 0.1262
IMD5 1.2553 0.7680 2.0518 0.3645
IBD Duration 0.9966 0.9833 1.0101 0.6195
BMI 0.9917 0.9680 1.0160 0.4995
TreatmentMono biologic 0.7984 0.5071 1.2569 0.3308
TreatmentCombo therapy 0.3892 0.1897 0.7987 0.0101
Treatment5-ASA 1.2486 0.8864 1.7588 0.2041
TreatmentNone reported 1.0493 0.7418 1.4842 0.7856
Age 0.9886 0.9796 0.9977 0.0138
vas_control85+ 0.7444 0.5813 0.9534 0.0194
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 3.0657 0.9987 0.0798
cat 4.1583 1.9963 0.1247
IMD 3.6741 3.9909 0.4505
IBD Duration 1.8560 0.9986 0.1728
BMI 0.4419 0.9983 0.5055
Treatment 2.1132 3.9812 0.7123
Age 0.1468 0.9952 0.6997
vas_control 3.2175 0.9984 0.0727
GLOBAL 18.8708 15.8740 0.2681 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal extent
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Extent +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3767 1.0672 1.7758 0.0139
catFC 50-250 1.6148 1.2099 2.1551 0.0011
catFC > 250 2.2051 1.6238 2.9943 0.0000
IMD2 1.2645 0.7614 2.1001 0.3646
IMD3 0.9215 0.5585 1.5204 0.7490
IMD4 1.3271 0.8327 2.1149 0.2340
IMD5 1.0640 0.6676 1.6957 0.7943
IBD Duration 1.0004 0.9867 1.0142 0.9584
BMI 0.9778 0.9542 1.0019 0.0706
TreatmentMono biologic 0.7052 0.4300 1.1567 0.1666
TreatmentCombo therapy 0.4102 0.2024 0.8314 0.0134
Treatment5-ASA 1.2768 0.8824 1.8474 0.1949
TreatmentNone reported 1.1744 0.8077 1.7074 0.3999
Age 0.9904 0.9810 1.0000 0.0490
ExtentE2 1.0307 0.7316 1.4520 0.8629
ExtentE3 0.7579 0.5150 1.1154 0.1597
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0569 0.9882 0.8075
cat 5.1167 1.9643 0.0749
IMD 2.5677 3.9426 0.6237
IBD Duration 3.4467 0.9871 0.0622
BMI 0.3924 0.9896 0.5267
Treatment 2.8080 3.8644 0.5690
Age 0.9359 0.9668 0.3222
Extent 3.7382 1.9695 0.1506
GLOBAL 18.1260 24.7328 0.8269 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Mayo score
Code
fit.me <- coxph(
  Surv(softflare_time, softflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Mayo +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4398 1.1189 1.8528 0.0046
catFC 50-250 1.5945 1.1950 2.1276 0.0015
catFC > 250 2.2595 1.6425 3.1083 0.0000
IMD2 1.3074 0.7584 2.2539 0.3346
IMD3 0.9626 0.5670 1.6342 0.8878
IMD4 1.3455 0.8164 2.2175 0.2444
IMD5 1.0652 0.6483 1.7502 0.8030
IBD Duration 1.0001 0.9867 1.0136 0.9915
BMI 0.9854 0.9625 1.0088 0.2198
TreatmentMono biologic 0.5915 0.3572 0.9796 0.0413
TreatmentCombo therapy 0.4295 0.2134 0.8644 0.0179
Treatment5-ASA 1.2812 0.8937 1.8366 0.1776
TreatmentNone reported 0.9580 0.6649 1.3804 0.8179
Age 0.9889 0.9792 0.9987 0.0267
Mayo 1.1328 1.0375 1.2367 0.0054
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 1.3866 0.9922 0.2369
cat 4.6541 1.9880 0.0966
IMD 2.8186 3.9673 0.5835
IBD Duration 2.3938 0.9918 0.1205
BMI 0.8275 0.9925 0.3604
Treatment 9.9460 3.8974 0.0385
Age 0.3657 0.9809 0.5374
Mayo 2.2088 0.9552 0.1295
GLOBAL 24.8425 18.9834 0.1651 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Objective flare

IBD Control-8
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    control_8 +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4361 1.0471 1.9695 0.0247
catFC 50-250 2.1254 1.4710 3.0711 0.0001
catFC > 250 2.9427 1.9938 4.3433 0.0000
IMD2 1.3172 0.6253 2.7747 0.4686
IMD3 1.3871 0.6873 2.7994 0.3611
IMD4 2.0933 1.0736 4.0816 0.0301
IMD5 1.5672 0.7967 3.0828 0.1930
IBD Duration 0.9950 0.9768 1.0135 0.5910
BMI 0.9895 0.9584 1.0216 0.5162
TreatmentMono biologic 1.1950 0.7038 2.0288 0.5096
TreatmentCombo therapy 0.8357 0.3995 1.7482 0.6336
Treatment5-ASA 1.1050 0.7089 1.7226 0.6593
TreatmentNone reported 0.7765 0.4906 1.2290 0.2803
Age 0.9880 0.9760 1.0002 0.0548
control_8 0.9486 0.9041 0.9953 0.0315
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0637 0.9892 0.7969
cat 3.4772 1.9688 0.1716
IMD 1.8773 3.9418 0.7505
IBD Duration 0.7156 0.9840 0.3917
BMI 0.0077 0.9829 0.9269
Treatment 20.9342 3.8558 0.0003
Age 0.8171 0.9731 0.3565
control_8 1.8868 0.9842 0.1663
GLOBAL 31.3141 25.4810 0.1966 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

IBD Control VAS
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    vas_control +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.4290 1.0409 1.9616 0.0272
catFC 50-250 2.1609 1.4984 3.1165 0.0000
catFC > 250 2.8718 1.9454 4.2392 0.0000
IMD2 1.2678 0.6018 2.6709 0.5325
IMD3 1.3359 0.6624 2.6944 0.4184
IMD4 1.9847 1.0210 3.8580 0.0433
IMD5 1.5216 0.7745 2.9894 0.2231
IBD Duration 0.9950 0.9768 1.0135 0.5912
BMI 0.9928 0.9616 1.0251 0.6584
TreatmentMono biologic 1.1935 0.7038 2.0240 0.5116
TreatmentCombo therapy 0.8317 0.3978 1.7392 0.6244
Treatment5-ASA 1.0667 0.6869 1.6566 0.7736
TreatmentNone reported 0.7580 0.4792 1.1988 0.2361
Age 0.9880 0.9759 1.0002 0.0532
vas_control85+ 0.7035 0.5116 0.9675 0.0305
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0865 0.9894 0.7648
cat 3.3658 1.9713 0.1818
IMD 1.9472 3.9420 0.7375
IBD Duration 0.6698 0.9851 0.4075
BMI 0.0145 0.9841 0.9003
Treatment 20.9874 3.8597 0.0003
Age 1.0160 0.9736 0.3050
vas_control 0.1703 0.9876 0.6748
GLOBAL 29.9586 24.8833 0.2210 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Montreal extent
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Extent +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.3650 0.9935 1.8753 0.0549
catFC 50-250 1.9498 1.3437 2.8292 0.0004
catFC > 250 3.1515 2.1596 4.5989 0.0000
IMD2 1.5489 0.7695 3.1178 0.2203
IMD3 1.5502 0.7933 3.0293 0.1996
IMD4 2.0892 1.0952 3.9852 0.0254
IMD5 1.6045 0.8339 3.0872 0.1568
IBD Duration 0.9946 0.9767 1.0128 0.5578
BMI 0.9789 0.9488 1.0100 0.1817
TreatmentMono biologic 1.1601 0.6890 1.9531 0.5764
TreatmentCombo therapy 0.8549 0.4262 1.7151 0.6590
Treatment5-ASA 1.1528 0.7407 1.7944 0.5287
TreatmentNone reported 0.8009 0.4979 1.2883 0.3599
Age 0.9900 0.9781 1.0021 0.1056
ExtentE2 2.1019 1.2247 3.6075 0.0070
ExtentE3 1.8335 1.0303 3.2631 0.0393
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.8817 0.9846 0.3425
cat 4.4095 1.9650 0.1070
IMD 2.7392 3.9419 0.5933
IBD Duration 0.0285 0.9837 0.8613
BMI 0.1230 0.9875 0.7209
Treatment 21.7625 3.8321 0.0002
Age 0.6202 0.9745 0.4212
Extent 0.7788 1.9770 0.6720
GLOBAL 32.5691 25.2863 0.1512 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Mayo score
Code
fit.me <- coxph(
  Surv(hardflare_time, hardflare) ~
    Sex +
    cat +
    IMD +
    `IBD Duration` +
    BMI +
    Treatment +
    Age +
    Mayo +
    frailty(SiteNo),
  control = coxph.control(outer.max = 20),
  data = flare.uc.df
)

invisible(cox_summary(fit.me))

Cox model summary:

Variable HR Lower 95% Upper 95% P-value
SexFemale 1.2152 0.8974 1.6454 0.2076
catFC 50-250 1.9271 1.3448 2.7617 0.0004
catFC > 250 3.0280 2.0635 4.4433 0.0000
IMD2 1.4994 0.7699 2.9201 0.2336
IMD3 1.1379 0.5902 2.1942 0.6997
IMD4 1.4786 0.7901 2.7670 0.2213
IMD5 1.2697 0.6779 2.3780 0.4558
IBD Duration 0.9998 0.9827 1.0173 0.9834
BMI 0.9887 0.9605 1.0178 0.4432
TreatmentMono biologic 1.0596 0.6310 1.7793 0.8268
TreatmentCombo therapy 1.1105 0.5812 2.1217 0.7510
Treatment5-ASA 1.0246 0.6646 1.5797 0.9122
TreatmentNone reported 0.7025 0.4458 1.1070 0.1281
Age 0.9868 0.9748 0.9990 0.0344
Mayo 1.1144 1.0001 1.2418 0.0498
Proportional hazards assumption test
Chi-squared statistic DF P-value
Sex 0.0003 0.9900 0.9849
cat 2.1046 1.9869 0.3464
IMD 1.5180 3.9641 0.8193
IBD Duration 0.1208 0.9907 0.7245
BMI 0.0302 0.9898 0.8591
Treatment 24.2076 3.8735 0.0001
Age 0.1736 0.9876 0.6719
Mayo 0.1053 0.9618 0.7306
GLOBAL 30.0448 19.2924 0.0561 
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

Reproduction and reproducibility

Session info

R version 4.4.0 (2024-04-24)

Platform: aarch64-unknown-linux-gnu

locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: gtsummary(v.1.7.2), DescTools(v.0.99.54), finalfit(v.1.0.7), coxme(v.2.2-20), bdsmatrix(v.1.3-7), pander(v.0.6.5), survminer(v.0.4.9), ggpubr(v.0.6.0), survival(v.3.5-8), datefixR(v.1.6.1), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1), tidyverse(v.2.0.0) and readxl(v.1.4.3)

loaded via a namespace (and not attached): gridExtra(v.2.3), gld(v.2.6.6), rlang(v.1.1.3), magrittr(v.2.0.3), e1071(v.1.7-14), compiler(v.4.4.0), mgcv(v.1.9-1), systemfonts(v.1.3.1), vctrs(v.0.6.5), pkgconfig(v.2.0.3), shape(v.1.4.6.1), fastmap(v.1.2.0), backports(v.1.5.0), labeling(v.0.4.3), KMsurv(v.0.1-5), utf8(v.1.2.4), rmarkdown(v.2.27), markdown(v.1.12), tzdb(v.0.4.0), nloptr(v.2.0.3), ragg(v.1.3.2), xfun(v.0.44), glmnet(v.4.1-8), jomo(v.2.7-6), jsonlite(v.1.8.8), pan(v.1.9), broom(v.1.0.6), R6(v.2.5.1), stringi(v.1.8.4), car(v.3.1-2), boot(v.1.3-30), rpart(v.4.1.23), cellranger(v.1.1.0), Rcpp(v.1.0.12), iterators(v.1.0.14), knitr(v.1.47), zoo(v.1.8-12), Matrix(v.1.7-0), splines(v.4.4.0), nnet(v.7.3-19), timechange(v.0.3.0), tidyselect(v.1.2.1), rstudioapi(v.0.16.0), abind(v.1.4-5), yaml(v.2.3.8), ggtext(v.0.1.2), codetools(v.0.2-20), plyr(v.1.8.9), lattice(v.0.22-6), withr(v.3.0.0), evaluate(v.0.23), proxy(v.0.4-27), xml2(v.1.3.6), survMisc(v.0.5.6), pillar(v.1.9.0), carData(v.3.0-5), mice(v.3.16.0), foreach(v.1.5.2), generics(v.0.1.3), hms(v.1.1.3), commonmark(v.1.9.1), munsell(v.0.5.1), scales(v.1.3.0), rootSolve(v.1.8.2.4), minqa(v.1.2.7), xtable(v.1.8-4), class(v.7.3-22), glue(v.1.7.0), lmom(v.3.0), tools(v.4.4.0), data.table(v.1.15.4), lme4(v.1.1-35.3), ggsignif(v.0.6.4), Exact(v.3.2), mvtnorm(v.1.2-5), grid(v.4.4.0), colorspace(v.2.1-0), nlme(v.3.1-164), cli(v.3.6.2), km.ci(v.0.5-6), textshaping(v.0.4.0), fansi(v.1.0.6), expm(v.0.999-9), broom.helpers(v.1.15.0), gt(v.0.10.1), gtable(v.0.3.5), rstatix(v.0.7.2), digest(v.0.6.35), farver(v.2.1.2), htmlwidgets(v.1.6.4), htmltools(v.0.5.8.1), lifecycle(v.1.0.4), httr(v.1.4.7), mitml(v.0.4-5), gridtext(v.0.1.5) and MASS(v.7.3-60.2)

Licensed by CC BY unless otherwise stated.