Using Interpretable Machine Learning to Identify Baseline Predictive Factors of Remission and Drug Durability in Crohn’s Disease Patients on Ustekinumab

Citation: Chaparro, M.; Baston-Rey, I.; Fernández Salgado, E.; González García, J.; Ramos, L.; Diz-Lois Palomares, M.T.; Argüelles-Arias, F.; Iglesias Flores, E.; Cabello, M.; Rubio Iturria, S.; et al. Using Interpretable Machine Learning to Identify Baseline Predictive Factors of Remission and Drug Durability in Crohn’s Disease Patients on Ustekinumab. J. Clin. Med. 2022,11, 4518. https:// doi.org/10.3390/jcm11154518 Academic Editor: Gary Edward Wild Received: 10 June 2022 Accepted: 29 July 2022 Published: 3 August 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Journal of Clinical Medicine Article Using Interpretable Machine Learning to Identify Baseline Predictive Factors of Remission and Drug Durability in Crohn’s Disease Patients on Ustekinumab María Chaparro 1,2 , Iria Baston-Rey 3, Estela Fernández Salgado 4, Javier González García5, Laura Ramos 6, María Teresa Diz-Lois Palomares 7, Federico Argüelles-Arias 8,9 , Eva Iglesias Flores 10, Mercedes Cabello 11, Saioa Rubio Iturria 12, Andrea Núñez Ortiz 13, Mara Charro 14, Daniel Ginard 15, Carmen Dueñas Sadornil 16, Olga Merino Ochoa 17, David Busquets 18, Eduardo Iyo 19, Ana Gutiérrez Casbas 2,20, Patricia Ramírez de la Piscina 21, Marta Maia Boscá-Watts 22, Maite Arroyo 23 , María JoséGarcía24 , Esther Hinojosa 25, Jordi Gordillo 26, Pilar Martínez Montiel 27, Benito Velayos Jiménez 28, Cristina Quílez Ivorra 29, Juan María Vázquez Morón30, JoséMaría Huguet 31 , Yago González-Lama 32, Ana Isabel Muñagorri Santos 33, Víctor Manuel Amo 34, María Dolores Martín Arranz 35,36, Fernando Bermejo 36,37, Jesús Martínez Cadilla 38 , Cristina Rubín de Célix 1,2 , Paola Fradejas Salazar 39, Antonio López San Román40, Nuria Jiménez 41, Santiago García-López 42 , Anna Figuerola 43, Itxaso Jiménez 44, Francisco JoséMartínez Cerezo 45, Carlos Taxonera 46, Pilar Varela 47, Ruth de Francisco 48 , David Monfort 49, Gema Molina Arriero 50, Alejandro Hernández-Camba 51 , Francisco Javier García Alonso 52, Manuel Van Domselaar 53, Ramón Pajares-Villarroya 54, Alejandro Núñez 55, Francisco Rodríguez Moranta 56, Ignacio Marín-Jiménez 57, Virginia Robles Alonso 58, María del Mar Martín Rodríguez 59, Patricia Camo-Monterde 60, Iván García Tercero 61, Mercedes Navarro-Llavat 62 , Lara Arias García63, Daniel Hervías Cruz 64, Sebastian Kloss 65, Alun Passey 65, Cynthia Novella 66,*, Eugenia Vispo 66, Manuel Barreiro-de Acosta 3and Javier P. Gisbert 1,2 1Hospital Universitario de La Princesa, Instituto de Investigación Sanitaria Princesa (IIS-IP), Universidad Autónoma de Madrid, 28006 Madrid, Spain; [email protected] (M.C.); [email protected] (C.R.d.C.); javier[email protected] (J.P.G.) 2Centro de Investigación Biomédica en Red de Enfermedades Hepáticas y Digestivas (CIBERehd), 28029 Madrid, Spain 3Complejo Hospitalario Universitario de Santiago, 15706 Santiago de Compostela, Spain; [email protected] (I.B.-R.); manubarreir[email protected] (M.B.-d.A.) 4Complejo Hospitalario de Pontevedra, 36164 Pontevedra, Spain; [email protected] 5Hospital Público Comarcal la Inmaculada, 04600 Almería, Spain; [email protected] 6Hospital Universitario de Canarias, 38320 Tenerife, Spain; [email protected] 7Hospital Universitario A Coruña, 15006 A Coruña, Spain; [email protected] 8Hospital Universitario Virgen Macarena, 41009 Seville, Spain; [email protected] 9Facultad de Medicina, Universidad de Sevilla, 41009 Seville, Spain 10 Hospital Universitario Reina Sofía, 14004 Córdoba, Spain; evaiflor[email protected] 11 Hospital Universitario Virgen de Valme, 41014 Seville, Spain; [email protected] 12 Complejo Hospitalario de Navarra, 31008 Pamplona, Spain; [email protected] 13 Hospital Universitario Virgen del Rocío, 41013 Seville, Spain; andr[email protected] 14 Hospital de Barbastro, 22300 Barbastro, Spain; [email protected] 15 Hospital Universitario Son Espases, 07120 Palma, Spain; [email protected] 16 Hospital San Pedro de Alcántara, 10003 Cáceres, Spain; [email protected] 17 Hospital Universitario Cruces, 48903 Barakaldo, Spain; [email protected] 18 Hospital Universitari de Girona Doctor Josep Trueta, 17007 Girona, Spain; dbusquets.gir[email protected] 19 Hospital Comarcal de Inca, 07300 Inca, Spain; [email protected] 20 Hospital General Universitario de Alicante, Instituto de Investigación Sanitaria y Biomédica de Alicante (ISABIAL), 03010 Alicante, Spain; [email protected] 21 Hospital Universitario de Araba-Txagorritxu, 01004 Vitoria-Gasteiz, Spain; [email protected] 22 Hospital Clínico Universitario de Valencia, 46010 Valencia, Spain; [email protected] 23 Hospital Clínico Universitario Lozano Blesa, 50009 Zaragoza, Spain; [email protected] 24 Hospital Universitario Marqués de Valdecilla, IDIVAL, 39008 Santander, Spain; [email protected] 25 Hospital de Manises, 46940 Manises, Spain; [email protected] 26 Hospital de la Santa Creu i Sant Pau, 08041 Barcelona, Spain; [email protected] 27 Hospital Universitario 12 de Octubre, 28041 Madrid, Spain; [email protected] 28 Hospital Clínico Universitario de Valladolid, 47003 Valladolid, Spain; [email protected] 29 Hospital Marina Baixa, 03570 Villajoyosa, Spain; [email protected] 30 Hospital Universitario Juan Ramón Jiménez, 21002 Huelva, Spain; [email protected] J. Clin. Med. 2022,11, 4518. https://doi.org/10.3390/jcm11154518 https://www.mdpi.com/journal/jcm J. Clin. Med. 2022,11, 4518 2 of 17 31 Hospital General Universitario de Valencia, 46014 Valencia, Spain; [email protected] 32 Hospital Universitario Puerta de Hierro, 28222 Majadahonda, Spain; [email protected] 33 Hospital Universitario de Donostia, 20014 Donostia-San Sebastián, Spain; [email protected] 34 Hospital Regional de Málaga, 29010 Málaga, Spain; [email protected] 35 Hospital Universitario de La Paz, 28046 Madrid, Spain; [email protected] 36 Instituto de Investigación Sanitaria del Hospital La Paz (IdiPaz), 28046 Madrid, Spain; [email protected] 37 Hospital Universitario de Fuenlabrada, 28942 Fuenlabrada, Spain 38 Hospital Alvaro Cunqueiro, 36213 Vigo, Spain; [email protected] 39 Hospital Virgen de la Concha, 49022 Zamora, Spain; [email protected] 40 Hospital Universitario Ramón y Cajal, 28034 Madrid, Spain; [email protected]g 41 Hospital General Universitario de Elche, 03203 Elche, Spain; [email protected] 42 Hospital Universitario Miguel Servet, 50009 Zaragoza, Spain; [email protected] 43 Hospital General Universitario de Castellón, 12004 Castellón de la Plana, Spain; [email protected] 44 Hospital Universitario de Galdakao-Usansolo, 48960 Galdakao, Spain; [email protected] 45 Hospital Universitario Sant Joan de Reus, 43204 Reus, Spain; [email protected] 46 Hospital Clínico Universitario San Carlos, Instituto de Investigación del Hospital Clínico San Carlos [IdISSC], 28040 Madrid, Spain; [email protected] 47 Hospital Universitario de Cabueñes, 33203 Gijón, Spain; [email protected] 48 Hospital Universitario Central de Asturias, Instituto de Investigación Biosanitaria del Principado de Asturias [ISPA], 33011 Oviedo, Spain; [email protected] 49 Consorci Sanitari de Terrassa, 08227 Terrassa, Spain; [email protected] 50 Complejo Hospitalario Universitario de Ferrol, 15405 Ferrol, Spain; [email protected] 51 Hospital Universitario Nuestra Señora de Candelaria, 38010 Santa Cruz de Tenerife, Spain; dr[email protected] 52 Hospital Universitario Rio Hortega, 47012 Valladolid, Spain; [email protected] 53 Hospital Universitario de Torrejón, 28850 Torrejón de Ardoz, Spain; [email protected] 54 Hospital Universitario Infanta Sofía, 28702 San Sebastián de los Reyes, Spain; [email protected]g 55 Hospital Clínico Universitario de Salamanca, 37007 Salamanca, Spain; [email protected] 56 Hospital Universitario Bellvitge, 08907 L’Hospitalet de Llobregat, Spain; [email protected] 57 Hospital Universitario Gregorio Marañón, 28009 Madrid, Spain; [email protected] 58 Hospital Universitario Vall d’Hebrón, 08035 Barcelona, Spain; vir[email protected] 59 Hospital Universitario Virgen de las Nieves, 18014 Granada, Spain; [email protected] 60 Hospital Universitario San Jorge, 22004 Huesca, Spain; [email protected] 61 Hospital General Universitario Santa Lucía, 30202 Cartagena, Spain; [email protected] 62 Hospital de Sant Joan DespíMoisès Broggi, 08970 Sant Joan Despí, Spain; mnavarr[email protected] 63 Hospital Universitario de Burgos, 09006 Burgos, Spain; [email protected] 64 Hospital General Universitario de Ciudad Real, 13005 Ciudad Real, Spain; [email protected] 65 Janssen, EMEA, 2340 Beerse, Belgium; [email protected] (S.K.); [email protected] (A.P.) 66 Janssen Medical Department, Paseo Doce Estrellas, 5-7, 28042 Madrid, Spain; [email protected] *Correspondence: [email protected] Abstract: Ustekinumab has shown efficacy in Crohn’s Disease (CD) patients. To identify patient profiles of those who benefit the most from this treatment would help to position this drug in the therapeutic paradigm of CD and generate hypotheses for future trials. The objective of this analysis was to determine whether baseline patient characteristics are predictive of remission and the drug durability of ustekinumab, and whether its positioning with respect to prior use of biologics has a significant effect after correcting for disease severity and phenotype at baseline using interpretable machine learning. Patients’ data from SUSTAIN, a retrospective multicenter single-arm cohort study, were used. Disease phenotype, baseline laboratory data, and prior treatment characteristics were documented. Clinical remission was defined as the Harvey Bradshaw Index ≤ 4 and was tracked longitudinally. Drug durability was defined as the time until a patient discontinued treatment. A total of 439 participants from 60 centers were included and a total of 20 baseline covariates considered. Less exposure to previous biologics had a positive effect on remission, even after controlling for baseline disease severity using a non-linear, additive, multivariable model. Additionally, age, body mass index, and fecal calprotectin at baseline were found to be statistically significant as independent negative risk factors for both remission and drug survival, with further risk factors identified for remission. Keywords: Crohn’s Disease; ustekinumab; predictive factors J. Clin. Med. 2022,11, 4518 3 of 17 1. Introduction Crohn’s Disease (CD) is a chronic inflammatory bowel disease that can result in ulceration, thickening of the intestinal wall, and a broad array of symptoms and complications. The immune mechanisms underlying CD are becoming better understood, and therefore more therapeutic options have appeared in recent years, including use of immunomodulators and biologic agents. Ustekinumab is a drug treatment in the latter category, targeting interleukin (IL)-12 and IL-23, that is indicated for the treatment of adult patients with moderately to severely active CD who have had an inadequate response with, lost response to, or were intolerant to either conventional therapy or a TNF α antagonist or have medical contraindications to such therapies [ 1 ]. The positioning of UST in particular versus anti-integrin medications is a topic of active discussion in the clinical literature [2,3]. The use of biomarkers, clinical data, and prediction models holds promise of finding the best drug for the right patient at the right time. Although clinical trials offer the gold standard by way of direct comparison between alternative treatments, they can also have certain limitations based on their generalizability (restrictive eligibility criteria, strict designs, and commonly binary endpoints). On the contrary, real-world studies and pragmatic observational trials can better reflect clinical practice and represent a greater diversity of patients, including ethnic background, biologic constituency, geographic representation, with sometimes longer follow ups and bigger cohort sizes, and hence are an important supplementary source of information to clinical trials. Real-world data from electronic medical records and claims data often come with very large sample sizes but are also constrained by data quality issues. In contrast, observational studies such as SUSTAIN [ 4 ] benefit from high-quality data capture and a pragmatic view of clinical practice. The objective of the current study was to determine whether baseline patient characteristics, including prior treatment history, are predictive of remission and drug durability after controlling for disease severity at baseline. 2. Materials and Methods The SUSTAIN study [ 4 ] is a retrospective multicenter study comprising 463 patients on ustekinumab (UST) with active CD (Harvey-Bradshaw Index (HBI) > 4). The HBI index was recorded by gastroenterologists in the follow-up of patients. HBI measurement was not available for 24 patients, who were excluded from the analysis, leaving a cohort size of 439 from 61 centers. Patients were evaluated at baseline and then during the follow-up for an average of 5.5 visits (median 5), and an average frequency of every 8 weeks, generating a total of 2876 visits in total where HBI was recorded (Figure 1). This creates an opportunity for longitudinal analysis as well as drug durability, which are the focus of this study. J. Clin. Med. 2022, 11, x FOR PEER REVIEW 4 of 18 1. Introduction Crohn’s Disease (CD) is a chronic inflammatory bowel disease that can result in ulceration, thickening of the intestinal wall, and a broad array of symptoms and complications. The immune mechanisms underlying CD are becoming better understood, and therefore more therapeutic options have appeared in recent years, including use of immunomodulators and biologic agents. Ustekinumab is a drug treatment in the latter category, targeting interleukin (IL)-12 and IL-23, that is indicated for the treatment of adult patients with moderately to severely active CD who have had an inadequate response with, lost response to, or were intolerant to either conventional therapy or a TNFα antagonist or have medical contraindications to such therapies [1]. The positioning of UST in particular versus anti-integrin medications is a topic of active discussion in the clinical literature [2,3]. The use of biomarkers, clinical data, and prediction models holds promise of finding the best drug for the right patient at the right time. Although clinical trials offer the gold standard by way of direct comparison between alternative treatments, they can also have certain limitations based on their generalizability (restrictive eligibility criteria, strict designs, and commonly binary endpoints). On the contrary, real-world studies and pragmatic observational trials can better reflect clinical practice and represent a greater diversity of patients, including ethnic background, biologic constituency, geographic representation, with sometimes longer follow ups and bigger cohort sizes, and hence are an important supplementary source of information to clinical trials. Real-world data from electronic medical records and claims data often come with very large sample sizes but are also constrained by data quality issues. In contrast, observational studies such as SUSTAIN [4] benefit from high-quality data capture and a pragmatic view of clinical practice. The objective of the current study was to determine whether baseline patient characteristics, including prior treatment history, are predictive of remission and drug durability after controlling for disease severity at baseline. 2. Materials and Methods The SUSTAIN study [4] is a retrospective multicenter study comprising 463 patients on ustekinumab (UST) with active CD (Harvey-Bradshaw Index (HBI) > 4). The HBI index was recorded by gastroenterologists in the follow-up of patients. HBI measurement was not available for 24 patients, who were excluded from the analysis, leaving a cohort size of 439 from 61 centers. Patients were evaluated at baseline and then during the follow-up for an average of 5.5 visits (median 5), and an average frequency of every 8 weeks, generating a total of 2876 visits in total where HBI was recorded (Figure 1). This creates an opportunity for longitudinal analysis as well as drug durability, which are the focus of this study. Figure 1. Number of unique visits over time ( left ) and number of follow-up visits per patient; red line indicates the mean number of visits (right). J. Clin. Med. 2022,11, 4518 4 of 17 The patient cohort exhibited a median HBI of 8.0 at baseline, with just under 50% with mild disease (5 ≤ HBI < 7), 45% with moderate disease (7 ≤ HBI ≤ 16), and 5% with severe disease (HBI > 16). The cohort was refractory, defined as having failed to respond to or remain on at least one other treatment line in the past, with 96.1% of the cohort being anti-TNF experienced and 23.7% also anti-integrin experienced prior to the onset of treatment with UST. A total of 60% had had at least one surgery; at baseline, 28% were treated with concomitant steroids and 30% with immunosuppressants (65% azathioprine, 28.8% methotrexate, and 8% Mercaptopurine). The analysis was primarily focused on remission as a dichotomous endpoint, achieving or not remission, defined as a HBI value of ≤ 4. Importantly, the analysis was longitudinal and considered all available follow-ups per patient, rather than focusing on a specific crosssection, with temporal dependence handled by adding “days since onset of treatment” as an independent variable. Drop-out bias was handled using the last-value-carried-forward (LVCF) imputation as detailed in Appendix A. As a secondary endpoint, drug durability was considered, defined as the time to dropout, which is right-censored by the last follow-up. In addition to its clinical significance, drug durability served to reinforce findings that replicate across both analyses, given its lack of sensitivity to drop-out bias. A complete table of all features used can be found in Table 1. These included 20 baseline features, time as a longitudinal covariate, and the number of concomitant steroids courses administered during the treatment window. Both models (remission and durability) were able to handle both discrete and numeric variables. For interpretability reasons, all numeric variables, except days since the start of UST treatment, were discretized using a model-driven technique that automatically determines the cut-offs that are optimal for the given multivariable analysis (see Appendix A). For example, the model-driven lower cut-off for the Body-Mass-Index (BMI) is 27, rather than the clinical threshold of 17. The groups employed for discretization are shown in the fourth column of Table 1. Following discretization, a separate model of remission was built on discretized variables only. In particular, the time elapsed between diagnosis and start of treatment, the age at the time of consent, BMI, and laboratory tests (baseline albumin, baseline fecal calprotectin, baseline haemoglobin, and baseline CRP), the number of surgeries prior to onset of treatment, the number of comorbidities, and the number of anti-TNF episodes, and similarly with the number of anti-integrin episodes were used as raw numeric values initially and discretized for interpretability. Treatment episodes refer to the number of times a patient has been treated with the respective drug class. So, for example, a patient with 2 treatment episodes on anti-integrin would have been started and stopped anti-integrin treatment twice, but in both cases, using vedolizumab, which is the only anti-integrin used in this cohort. Alternative codifications of exposure to anti-TNF or anti-integrin were tried out, including total duration of treatment (number of days on drug) and exposure (yes/no). All had a similar effect directionally on the model, but number of treatment episodes offered the best fit overall. J. Clin. Med. 2022,11, 4518 5 of 17 Table 1. Variables employed in their original form and discretized. Variable Original (Non-Discretized) Post-Processed (Discretized) Description of Discretized Version Descriptives Days since the start of UST treatment (longitudinal, not baseline) Numeric Numeric - HBI at baseline, categorized Ordinal (3) - [5–7] (mild), [8–16] (moderate), >16 (severe) Distribution: Mild 50%, Moderate 45%, Severe 5% % missing values: 0% Years between diagnosis and start of UST Numeric Dichotomous ≤18, >18 Median: 10.0 % missing values: 0% Number of surgeries before UST treatment onset Numeric Ordinal (3) 0, 1–2, ≥3 Distribution: 0 (47%), 1–2 (43%), ≥3 (10%) Number of anti-TNF episodes Numeric Dichotomous ≤1, >1 Median: 2.0 96.1% of patients ≥1 anti-TNF Number of anti-integrin episodes Numeric Dichotomous 0, 1 Median: 0.0 23.7% of patients ≥1 anti-integrin Age at the time of signing consent Numeric Dichotomous ≤40, >40 Median: 46.0 % missing values: 0% BMI at baseline Numeric Ordinal (3) <27, 27–38, >38 Median: 23.81 % missing values: 16% Sex Dichotomous Dichotomous Female, Male Distribution: Female 49%, Male 51% % missing values: 0% Number of comorbidities Numeric Ordinal (4) 0–1, 2–3, ≥4 Median: 1.0 Perianal disease Categorical Dichotomous Never/Previous, Current Distribution: Never/Previous 86%, Current 14% % missing values: 0% CD location Categorical Dichotomous Ileocolic/Ileum (L1/L3), Colon (L2) Distribution: Ileocolic/Ileum 88%, Colon 12% % missing values: 0% Patient ever had EIMs Dichotomous Dichotomous Yes/No Distribution: No 59%, Yes 41% Family history of CD Dichotomous Dichotomous Yes/No Distribution: No 13%, Yes 87% % missing values: 16% Upper gastrointenstinal tract (L4) involved Dichotomous Dichotomous Yes/No Distribution: No 92%, Yes 8% % missing values: 0% J. Clin. Med. 2022,11, 4518 6 of 17 Table 1. Cont. Variable Original (Non-Discretized) Post-Processed (Discretized) Description of Discretized Version Descriptives Baseline albumin (g/L) Numeric Dichotomous ≤3.8, >3.8 Median: 4.0 % missing values: 33% Baseline faecal calprotectin (mcg/g) Numeric Ordinal (3) ≤80, 81–650, >650 Median: 667.33 % missing values: 52% Baseline haemoglobin (g/L) Numeric Dichotomous ≤10.9, >10.9 Median: 13.2 % missing values: 10% Baseline CRP (mg/L) Numeric Ordinal (3) ≤0.5, 0.6–6.1, >6.1 Median: 7.2 % missing values: 10% Steroid use at first UST dose Dichotomous Dichotomous Yes/No Distribution: No 72%, Yes 28% Under Immunosuppressants at first UST dose Dichotomous Dichotomous Yes/No Distribution: No 70%, Yes 30% Number of concomitant steroid courses Numeric Dichotomous Yes/No Median: 0.0 BMI, Body Mass Index; CD, Crohn’s Disease; CRP, C-reactive protein; EIM, extraintestinal manifestation; HBI, Harvey Bradshaw Index; TNF, Tumor Necrosis Factor; UST, ustekinumab. J. Clin. Med. 2022,11, 4518 7 of 17 Baseline measurements were used for two reasons: to control for disease severity at the onset of treatment, but also as predictors of remission given information that is available at a clinically actionable point and could in the future inform a decision to treat by the clinician. HBI at baseline was the only numeric variable that was discretized pre-modelling into mild, moderate, and severe disease using the clinical definition described earlier, to facilitate investigation of its interactions with other variables. In particular, a large percentage of patients were under either steroids or immunosuppressants (or both) at the study onset, which might have masked what would have otherwise been a higher HBI value. These variables were considered as having interaction effects with baseline HBI. BMI at baseline was estimated using height and weight measurements where possible and inferred using regression imputation otherwise. Sex was also considered. Given the wide range of different comorbidities, these were not considered individually but rather a single variable was introduced counting the number of comorbidities present. Disease type was described by the following variables: whether the patient had perianal disease,CD location according to the Montreal classification, ileal (L1), colonic (L2), or ileocolonic (L3), whether the upper gastrointestinal tract (L4) was involved, and whether there was ever an extraintestinal manifestation.CD location was collapsed into a dichotomous variable capturing whether the ileum was involved (Yes/No), as it was observed that there was no statistically significant difference between the effect of ileocolonic (L3) vs. ileal only (L1). Additionally, a family history of CD was also encoded as a dichotomous variable. The patient journey was a semi-structured data source that can be summarized in a tabular form in a variety of ways. This analysis considered the number of surgeries before UST treatment onset, discretized into groups for 0, 1–2, and 3+, and treated like an ordinal factor. Previous biologic use was a key object of interest for the study. It was encoded as the number of anti-TNF episodes and, separately, the number anti-integrin episodes. Sensitivity analyses were conducted with alternative encodings including the cumulative length of treatment episodes measured in days, or measured as a proportion of time since diagnosis, and maintained a similar direction and significance levels. The number of concomitant steroid courses that were given concomitantly with UST treatment was also considered in the model as a control, to improve the generalizability of the findings and attempt to extrapolate on the possibility of steroid-free remission. Three separate analyses were run, namely, longitudinal remission non-discretized, longitudinal remission discretized, and drug durability. The primary analysis was a longitudinal remission analysis that measured the effect of different covariates on the probability of remission at different times since onset of treatment. This model was fitted on continuous variables first and then on discretized versions for interpretability reasons. The secondary analysis was a drug durability analysis focused on the time each patient stayed on ustekinumab. In the longitudinal analysis, every patient follow-up constituted a single observation, and the set of covariates included the time since onset of treatment, allowing for follow-ups made at different times to become comparable by controlling for the effect of time, while at the same time extracting an explicit representation of a non-linear trajectory of the endpoint over the duration of the treatment. Such use of time as a predictor notably increases the effective sample size of the analysis. Baseline covariates featured identical values for all observations coming from the same patient. LVCF techniques were used to correct for possible drop-out bias (see Appendix A). Given a dichotomous dependent variable, we were able to assess goodness-of-fit by evaluating R2 (proportion of variance explained) as well as a standard metric of classification accuracy, Area Under the Curve (AUC), using 10-fold cross validation, treating the remission model as a classifier. Drug durability featured instead one observation per patient and hence had less power than the longitudinal analysis but benefitted from a lack of sensitivity to drop-out bias. All three analyses relied on the Generalized Additive Model (GAM) framework, introduced [ 5 ] as an extension of linear and logistic regression that allows for data-driven discovery of non-linear effects between each predictor and the remission flag using splines, J. Clin. Med. 2022,11, 4518 8 of 17 while preserving its interpretable properties such as p-values and confidence intervals on additive effects. Survival analysis is a recent use case for GAMs proposed in [ 6 ] and has certain advantages over classical techniques such as Cox regression and Kaplan–Meier estimates, such as the possibility to describe the hazard function while at the same time consider non-linear effects of additional covariates via the use of splines. A key benefit of using the GAM modeling framework for both longitudinal [ 7 ] and durability analyses is that it allows for direct comparison of the effects of different covariates across the models. GAMs produce statistically accurate p-values if the number of spline knots (which controls the flexibility of each non-linear effect) is set in advance to a fixed value (in our case, this was 8 for main effects and 4 for interaction terms). Additionally, effect sizes can be read off a partial dependence plot which, for the longitudinal remission model, captures the effect on the log-odds of remission (on the y-axis) for changes in the underlying independent variable (x-axis). These plots and respective confidence bands were also used to guide variable discretization. Statistical analysis was performed in Python 3.0 using the PyGAM package (0.8.0) and Matplotlib 3.3.3 for plotting purposes. 3. Results 3.1. Performance of Models Performance of the longitudinal remission model in terms of cross-validated R 2 was modest, around 0.25, though, as discussed, a more reliable metric of predictive performance in this case is AUC, where the model performed well with a value of 0.796 (0.78, 0.8) over 10 folds of cross-validation, compared to that of a baseline model which only considered time as a predictor variable, whose median AUC was 0.62 (0.61–0.64) (Figure 2). The performance is comparable with similar analyses such as [ 8 ], though the latter additionally used laboratory values shortly after onset of treatment, rather than only at baseline as conducted here. J. Clin. Med. 2022, 11, x FOR PEER REVIEW 9 of 18 All three analyses relied on the Generalized Additive Model (GAM) framework, introduced [5] as an extension of linear and logistic regression that allows for data-driven discovery of non-linear effects between each predictor and the remission flag using splines, while preserving its interpretable properties such as p-values and confidence intervals on additive effects. Survival analysis is a recent use case for GAMs proposed in [6] and has certain advantages over classical techniques such as Cox regression and Kaplan–Meier estimates, such as the possibility to describe the hazard function while at the same time consider non-linear effects of additional covariates via the use of splines. A key benefit of using the GAM modeling framework for both longitudinal [7] and durability analyses is that it allows for direct comparison of the effects of different covariates across the models. GAMs produce statistically accurate p-values if the number of spline knots (which controls the flexibility of each non-linear effect) is set in advance to a fixed value (in our case, this was 8 for main effects and 4 for interaction terms). Additionally, effect sizes can be read off a partial dependence plot which, for the longitudinal remission model, captures the effect on the log-odds of remission (on the y-axis) for changes in the underlying independent variable (x-axis). These plots and respective confidence bands were also used to guide variable discretization. Statistical analysis was performed in Python 3.0 using the PyGAM package (0.8.0) and Matplotlib 3.3.3 for plotting purposes. 3. Results 3.1. Performance of Models Performance of the longitudinal remission model in terms of cross-validated R 2 was modest, around 0.25, though, as discussed, a more reliable metric of predictive performance in this case is AUC, where the model performed well with a value of 0.796 (0.78, 0.8) over 10 folds of cross-validation, compared to that of a baseline model which only considered time as a predictor variable, whose median AUC was 0.62 (0.61–0.64) (Figure 2). The performance is comparable with similar analyses such as [8], though the latter additionally used laboratory values shortly after onset of treatment, rather than only at baseline as conducted here. Given that information about the number of concomitant steroid courses is not available at the time of assessing a patient’s risk of not achieving remission, we reran the analysis without this variable in order to assess the sensitivity of the model to this variable, and the AUC dropped by less than 1%, from a median of 0.796 to one of 0.789. Figure 2. Ten iterations of a training-test split following which the longitudinal remission model is used to classify whether a patient will be in remission at a given week that was not observed in the training data. The full model makes use of all 22 variables (including non-significant ones), whereas the baseline model only uses time. Figure 2. Ten iterations of a training-test split following which the longitudinal remission model is used to classify whether a patient will be in remission at a given week that was not observed in the training data. The full model makes use of all 22 variables (including non-significant ones), whereas the baseline model only uses time. Given that information about the number of concomitant steroid courses is not available at the time of assessing a patient’s risk of not achieving remission, we reran the analysis without this variable in order to assess the sensitivity of the model to this variable, and the AUC dropped by less than 1%, from a median of 0.796 to one of 0.789. 3.2. Significant Variables The longitudinal remission model considered 22 variables out of which 20 were baseline features and 2 were not (days since start of UST treatment and number of concomitant steroid courses). Out of the baseline patient characteristics, 16 proved significant in the J. Clin. Med. 2022,11, 4518 9 of 17 non-discretized longitudinal remission model (Table 2). The discretized longitudinal remission model has by construction lower power and was fitted as an interpretability aid; it was performed only on those features that were already found significant from the non-discretized version of the primary analysis model. Most variables preserved their significance, but some lost it, reflecting the importance of discretizing numeric variables such as laboratory measurements. Reassuringly, the shapes, directions, and magnitudes of all reported effect sizes were preserved when discretizing, with the exception of baseline albumin and baseline hemoglobin which in the non-discretized version seem to have nonlinearities not sufficiently captured by the discretized versions, which likely is the reason why these features ceased to be significant in the discretized model. Indeed, several variables lose significance because of discretization, which is expected as discretization loses both information and power. Table 2. Statistical significance (p-value) of the effect of each variable in longitudinal remission and durability. Longitudinal Remission Durability Feature Non-Discretized Discretized * Days since start of UST treatment <0.05 <0.05 <0.05 HBI at the first UST dose, categorized <0.05 <0.05 0.083 Years between diagnosis and start of UST <0.05 0.547 0.077 Number of surgeries before UST treatment onset <0.05 <0.05 0.24 Number of anti-TNF episodes <0.05 0.09 <0.05 Number of anti-integrin episodes <0.05 <0.05 <0.05 Age at the time of signing consent <0.05 <0.05 <0.05 BMI at baseline <0.05 <0.05 <0.05 Sex 0.888 N/A 0.842 Number of comorbidities <0.05 <0.05 0.296 Perianal disease <0.05 <0.05 0.541 Location of CD <0.05 <0.05 0.137 Patient ever had EIMs 0.806 N/A 0.648 Family history of CD 0.958 N/A 0.998 Upper gastrointestinal tract (L4) involved 0.999 N/A 0.828 Baseline albumin <0.05 0.431 0.357 Baseline fecal calprotectin <0.05 0.26 <0.05 Baseline hemoglobin <0.05 0.77 <0.05 Baseline CRP <0.05 <0.05 <0.05 Under Steroids at first UST dose <0.05 <0.05 0.237 Under Immunosuppressants at first UST dose <0.05 <0.05 0.813 Number of concomitant Steroid courses <0.05 <0.05 <0.05 Variables that were significant in longitudinal remission (non-discretized) and durability models are indicated in bold. * All features discretized, except for “days since start of UST treatment”. BMI, Body Mass Index; CD, Crohn’s Disease; CRP, C-reactive protein; EIM, extraintestinal manifestation; HBI, Harvey Bradshaw Index; TNF, Tumor Necrosis Factor; UST, ustekinumab. A number of interaction terms were also included, as shown in Table A1. The drug durability model considered the same baseline variables and given its smaller effective sample size detected seven significant variables, all of which were reassuringly also significant in the remission model. J. Clin. Med. 2022,11, 4518 16 of 17 measurement and were then imputed using the hot-deck K-nearest neighbors imputation as per [14]. Appendix A.2. Drop-Out in Longitudinal Analysis Drop-out bias is a well-known challenge in longitudinal analysis. Analyses that have a non-linear estimation shape are more robust to such bias in the sense that it is unlikely to impact the shape of the remission trajectory earlier on in the treatment window but can still be affected at later parts of the trajectory where the prevalence of drop-out increases. In a longitudinal setting, additional timestamps were generated for drop-out patients using the LVCF imputation, so that each of these patients would not have a number of visits smaller than a value randomly drawn from the distribution of the number of doses received by non-drop-out patients. Appendix A.3. Shape Constraints The user can either manually pre-specify the number of splines (which control the smoothness of the function sj ) or use cross-validation to determine the optimal number instead, in which case, however, the p-values and confidence intervals are no longer valid. For this reason, splines were pre-set to fixed values. All continuous variables were modeled using splines, whereas factors were modeled as in logistic regression as constant effects. Modern implementations of GAM additionally allow the introduction of shape constraints into the smoothing functions, including monotonically decreasing, monotonically increasing, convex and concave shapes. Monotonicity is a natural generalization of linearity. For each variable, a monotonic decreasing and a monotonic increasing constraint were tried out, and the one resulting in the best fitting model was used, presented in Table 3. Appendix A.4. Interaction Effects In all models, a number of interaction effects were introduced in the model only where the main effects were first found to be significant: these were an interaction between hemoglobin and sex, interactions between each baseline lab and the HBI at the first UST dose, and, for longitudinal remission only, an interaction term between the days since the start of treatment and years since diagnosis, as well as an interaction between days since the start of UST treatment and BMI. These improved the fit of the model without modifying the shape or direction of the main effects. Table A1. Interaction effects. Feature Interactions p-Value of Effect Longitudinal Remission Durability Numeric Values Used for Several Features All Features Discretized, Except for “Days since Start of UST Treatment” (Baseline hemoglobin) ×(Sex) 0.626 N/A 0.821 (HBI at the first UST dose, categorized) ×(Baseline fecal calprotectin) <0.05 N/A 0.307 (HBI at the first UST dose, categorized) × (Baseline albumin) <0.05 N/A 0.239 (HBI at the first UST dose, categorized) × (Baseline hemoglobin) <0.05 N/A 0.306 (HBI at the first UST dose, categorized) ×(Baseline CRP) <0.05 N/A <0.05 (Days since the start of UST treatment) ×(Years between diagnosis and start of UST) <0.05 N/A N/A (Days since the start of UST treatment) ×(BMI) <0.05 N/A N/A J. Clin. Med. 2022,11, 4518 17 of 17 References 1. Feagan, B.G.; Sandborn, W.J.; Gasink, C.; Jacobstein, D.; Lang, Y.; Friedman, J.R.; Blank, M.A.; Johanns, J.; Gao, L.L.; Miao, Y. ; et al. Ustekinumab as induction and maintenance therapy for Crohn’s disease. N. Engl. J. Med. 2016 ,375, 1946–1960. [CrossRef] [PubMed] 2. Danese, S.; Bonovan, S.; Peyrin-Biroulet, L. Positioning ustekinumab in Crohn’s disease: From clinical evidence to clinical practice. J. Crohn’s Colitis 2017,11, 1258–1266. [CrossRef] [PubMed] 3. Pacou, M.; Mesana, L.; Gauthier, A.; Naessens, D.; Sloan, S.; Danese, S.; Bonovas, S.; Abrams, K. Indirect treatment comparison of ustekinumab versus other biologics in moderate to severe Crohn’s disease: A 1-year treatment sequence analysis. Value Health 2016,19, 576. [CrossRef] 4. Chaparro, M.; Baston-Rey, I.; Fernández-Salgado, E.; González García, J.; Ramos, L.; Diz-Lois Palomares, M.T.; Argüelles-Arias, F.; Iglesias Flores, E.; Cabello, M.; Rubio Iturria, S.; et al. Long-Term Real-World Effectiveness and Safety of Ustekinumab in Crohn’s Disease Patients: The SUSTAIN Study. Inflamm. Bowel Dis. 2022; online ahead of print. 5. Hastie, T.J.; Tibshirani, R.J. Generalized Additive Models; CRC Press: Boca Raton, FL, USA, 1990; Volume 43. 6. Argyropoulos, C.; Unruh, M.L. Analysis of time to event outcomes in randomized controlled trials by generalized additive models. PLoS ONE 2015,10, e0123784. [CrossRef] 7. Berhane, K.; Tibshirani, R.J. Generalized Additive Models for Longitudinal Data. Can. J. Stat. 1998,26, 517–535. [CrossRef] 8. Waljee, A.K.; Wallace, B.I.; Cohen-Mekelburg, S.; Liu, Y.; Liu, B.; Sauder, K.; Stidham, R.W.; Zhu, J.; Higgins, P.D. Development and validation of machine learning models in prediction of remission in patients with moderate to severe Crohn disease. JAMA Netw. Open 2019,2, e193721. [CrossRef] [PubMed] 9. Biemans, V.B.; van der Meulen-de Jong, A.E.; Van Der Woude, C.J.; Löwenberg, M.; Dijkstra, G.; Oldenburg, B.; De Boer, N.K.; Van Der Marel, S.; Bodelier, A.G.; Jansen, J.M.; et al. Ustekinumab for Crohn’s Disease: Results of the ICC Registry, a Nationwide Prospective Observational Cohort Study. J. Crohn’s Colitis 2020,14, 33–45. [CrossRef] [PubMed] 10. Straatmijer, T.; Biemans, V.B.; Hoentjen, F.; De Boer, N.K.; Bodelier, A.G.; Dijkstra, G.; Van Dop, W.A.; Haans, J.J.; Jansen, J.M.; Maljaars, P.J.; et al. Ustekinuma b for Crohn’s Disease: Two-Year Results of the Initiative on Crohn and Colitis (ICC) Registry, a Nationwide Prospective Observational Cohort Study. J. Crohn’s Colitis 2021,15, 1920–1930. [CrossRef] [PubMed] 11. Iborra, M.; Beltrán, B.; Fernández-Clotet, A.; Iglesias-Flores, E.; Navarro, P.; Rivero, M.; Gutiérrez, A.; Sierra-Ausin, M.; Mesonero, F. ; Ferreiro-Iglesias, R.; et al. Real-world long-term effectiveness of ustekinumab in Crohn’s disease: Results from the ENEIDA registry. Aliment. Pharm. Ther. 2020,52, 1017–1030. 12. Park, J.; Chun, J.; Yoon, H.; Cheon, J.H. Feasibility of a Clinical Decision Support Tool for Ustekinumab to Predict Clinical Remission and Relapse in Patients with Crohn’s Disease: A Multicenter Observational Study. Inflamm. Bowel Dis. 2022 , izac105. [CrossRef] [PubMed] 13. Little, R.J.; Rubin, D.B. Statistical Analysis with Missing Data; John Wiley and Sons: Hoboken, NJ, USA, 2019. 14. Garcia-Laencina, P.J.; Sancho-Gomez, J.-L.; Figueiras-Vidal, A.R.; Verleysen, M. KNearest Neighbours with Mutual Information for Simultaneous Classification and Missing Data Imputation. Neurocomputing 2009,72, 1483–1493. [CrossRef]