Science, Faculty of
Statistics, Department of
DSpace
UBCV
Kempf, Graeme
Master of Science - MSc
University of British Columbia
Multiple sclerosis (MS) is an autoimmune disease that progresses over the course of years. MS is treated with disease-modifying drugs (DMDs). The efficacy of DMDs was established in clinical trials that were short and excluded older individuals and individuals living with comorbidities. This has led to a lack of knowledge of the effects of chronic DMD use and the effects of DMDs on individuals that do not meet the traditional eligibility criteria for clinical trials. This thesis applies multi-state models, which can advance the understanding of a disease beyond that offered by time-to-event models alone, to administrative healthcare data in a retrospective study to explore the long-term, real-world efficacy of DMDs. Chapter 2 applies time-to-event models and weighted cumulative exposure to explore the association between DMDs and mortality as well as to explore reverse causality. Chapters 3 and 4 implement an illness-death model to investigate the main research questions of this thesis: whether exposure to any DMD is associated with fewer hospitalizations, shorter hospitalizations, and/or a reduction in the chance of dying inside or outside the hospital. Specifically, in Chapter 3, Cox-type regressions estimate the hazard ratios (HRs) for exposure to any DMDs. In Chapter 4, pseudo-value regression models examine whether exposure is associated with more time spent alive and out of the hospital or more time hospitalized. Chapters 3 and 4 use interactions between DMD and comorbidity to analyze if the impact of DMDs changed across comorbidity burden. The change of the effects of exposure over time is also investigated in Chapters 2, 3, and 4. The results suggest that DMD use is initially associated with fewer hospitalizations but long-term DMD exposure is associated with more hospitalizations. Long-term DMD exposure is also associated with a greater hazard of dying outside the hospital. Pseudo-value regression suggested that exposure amongst individuals with high comorbidity burden is associated with more time spent alive and out of the hospital than it is for individuals with no comorbidity.
The impact of disease-modifying drugsfor multiple sclerosis onhospitalizations and mortality inBritish ColumbiaA retrospective study using an illness-death multi-statemodelbyGraeme KempfB.B.A., Wilfrid Laurier University, 2020B.Sc., Wilfrid Laurier University, 2022A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Statistics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2024© Graeme Kempf 2024The following individuals certify that they have read, and recommend tothe Faculty of Graduate and Postdoctoral Studies for acceptance, the thesisentitled:The impact of disease-modifying drugs for multiplesclerosis on hospitalizations and mortality in BritishColumbia: a retrospective study using an illness-death multi-state modelsubmitted by Graeme Kempf in partial fulfillment of the requirements forthe degree of Master of Science in StatisticsExamining Committee:John Petkau, Professor Emeritus, Department of Statistics, UBCCo-supervisorYinshan Zhao, Research Associate, Department of Medicine (Neurology),Senior Data Scientist, Population Data BC, UBCCo-supervisorLang Wu, Professor, Department of Statistics, UBCSupervisory Committee MemberiiAbstractMultiple sclerosis (MS) is an autoimmune disease that progresses over thecourse of years. MS is treated with disease-modifying drugs (DMDs). Theefficacy of DMDs was established in clinical trials that were short and ex-cluded older individuals and individuals living with comorbidities. This hasled to a lack of knowledge of the effects of chronic DMD use and the effectsof DMDs on individuals that do not meet the traditional eligibility criteriafor clinical trials. This thesis applies multi-state models, which can advancethe understanding of a disease beyond that offered by time-to-event modelsalone, to administrative healthcare data in a retrospective study to explorethe long-term, real-world efficacy of DMDs.Chapter 2 applies time-to-event models and weighted cumulative exposureto explore the association between DMDs and mortality as well as to explorereverse causality. Chapters 3 and 4 implement an illness-death model to in-vestigate the main research questions of this thesis: whether exposure toany DMD is associated with fewer hospitalizations, shorter hospitalizations,and/or a reduction in the chance of dying inside or outside the hospital.Specifically, in Chapter 3, Cox-type regressions estimate the hazard ratios(HRs) for exposure to any DMDs. In Chapter 4, pseudo-value regressionmodels examine whether exposure is associated with more time spent aliveand out of the hospital or more time hospitalized. Chapters 3 and 4 use in-teractions between DMD and comorbidity to analyze if the impact of DMDschanged across comorbidity burden. The change of the effects of exposureover time is also investigated in Chapters 2, 3, and 4.The results suggest that DMD use is initially associated with fewer hospital-izations but long-term DMD exposure is associated with more hospitaliza-tions. Long-term DMD exposure is also associated with a greater hazard ofdying outside the hospital. Pseudo-value regression suggested that exposureamongst individuals with high comorbidity burden is associated with moretime spent alive and out of the hospital than it is for individuals with nocomorbidity.iiiLay SummaryEven though many outcomes are relevant when studying diseases, most stud-ies on drug efficacy investigate a drug’s effect on only a single outcome.Multi-state modelling techniques take a wider view and allow the drug’seffects on multiple events to be quantified which can elucidate a greater un-derstanding of the impact of the drug.This thesis applies multi-state modelling to healthcare data of British Columbiaresidents living with multiple sclerosis. In doing so, the association betweendrugs that are prescribed to treat multiple sclerosis and three outcomes areexplored: the frequency of hospitalizations, the lengths of hospital visits,and the chances of dying. The impact of these drugs for individuals livingwith comorbidities is studied as well.ivPrefaceThis thesis used as its foundation work already completed as part of the“Disease-modifying Drug Safety and Effectiveness in Multiple Sclerosis”(DRUMS) project by Dr. Helen Tremlett’s Pharmacoepidemiology in Mul-tiple Sclerosis (PiMS) research lab. The work in this thesis extends theanalysis of two papers [1], [2] by Dr. Huah Shin Ng and co-authors by ex-ploring the outcomes together where they were previously only consideredseparately.This thesis is an original, unpublished work which I, Graeme Kempf, au-thored under the supervision of Drs. John Petkau and Yinshan Zhao, whohad the idea of this project (application of multi-state models to the hospitaland death data of individuals living with MS) before I joined as a researchassistant. I formulated the research questions to be addressed, identified thespecific analysis approaches, implemented these approaches into code, andwrote up the results.Access to data provided by the Data Stewards is subject to approval butcan be requested for research projects through the Data Stewards or theirdesignated service providers. The following data sets were used in thisstudy: Medical Services Plan (MSP) Payment Information File, DischargeAbstracts Database, BC Vital Events and Statistics - Deaths, Pharmanet,Central Demographics File (MSP Registration and Premium Billings, ClientRoster and Census Geodata). You can find further information regardingthese data sets by visiting the PopData project webpage at: https://my.popdata.bc.ca/project listings/18-120/collection approval dates. All inferences,opinions, and conclusions drawn in this publication are those of the au-thor(s), and do not reflect the opinions or policies of the Data Steward(s).vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Study cohort . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Exposure to disease-modifying drugs . . . . . . . . . . . . . . 91.2.1 A note about treatment switching . . . . . . . . . . . 91.2.2 Exposure measures . . . . . . . . . . . . . . . . . . . 101.2.3 Rationale for using progressive measures of exposure 121.2.4 Description of exposure . . . . . . . . . . . . . . . . . 142 Time-to-event models . . . . . . . . . . . . . . . . . . . . . . . 172.1 Intensity-based modelling . . . . . . . . . . . . . . . . . . . . 172.1.1 Intensity functions . . . . . . . . . . . . . . . . . . . . 172.1.2 Cox regression . . . . . . . . . . . . . . . . . . . . . . 182.2 Non-parametric analysis of mortality . . . . . . . . . . . . . 232.3 The association between disease-modifying drugs and mortal-ity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.3.1 Limitations due to group sizes . . . . . . . . . . . . . 282.3.2 Model selection . . . . . . . . . . . . . . . . . . . . . 292.3.3 Some comments on interpreting results . . . . . . . . 302.3.4 Results of survival analysis . . . . . . . . . . . . . . . 30viTable of Contents2.3.5 Discussion of survival models . . . . . . . . . . . . . . 312.4 Factors associated with exposure . . . . . . . . . . . . . . . . 332.4.1 Results of the time-to-exposure model . . . . . . . . . 342.5 Exploring reverse causality using weighted cumulative expo-sure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.5.1 Results of weighted cumulative exposure Cox model . 373 An illness-death model to investigate the association be-tween disease-modifying drugs and hospitalizations and mor-tality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.1 Multi-state models . . . . . . . . . . . . . . . . . . . . . . . . 403.1.1 Illness-Death Model . . . . . . . . . . . . . . . . . . . 423.2 Non-parametric estimates . . . . . . . . . . . . . . . . . . . . 433.2.1 Cumulative transition intensity . . . . . . . . . . . . 433.2.2 Transition probabilities . . . . . . . . . . . . . . . . . 433.2.3 State occupancy probabilities . . . . . . . . . . . . . 443.3 Intensity-based multi-state model analysis . . . . . . . . . . 453.3.1 Likelihood . . . . . . . . . . . . . . . . . . . . . . . . 453.4 Non-parametric analysis . . . . . . . . . . . . . . . . . . . . . 513.5 Cox models for the illness-death model . . . . . . . . . . . . 553.5.1 O → H transition . . . . . . . . . . . . . . . . . . . . 563.5.2 O → D transition . . . . . . . . . . . . . . . . . . . . 583.5.3 H → O transition . . . . . . . . . . . . . . . . . . . . 593.5.4 H → D transition . . . . . . . . . . . . . . . . . . . . 613.6 Results of multi-state intensity-based analysis . . . . . . . . 623.7 Discussion of intensity-based analysis . . . . . . . . . . . . . 644 Analyzing the effect of exposure on the mean time spent inand out of the hospital using pseudo-values . . . . . . . . . 684.1 Mean sojourn time . . . . . . . . . . . . . . . . . . . . . . . . 684.2 Pseudo-values . . . . . . . . . . . . . . . . . . . . . . . . . . 694.3 Pseudo-value regression models . . . . . . . . . . . . . . . . . 714.4 Results from the pseudo-value models . . . . . . . . . . . . . 754.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 805.1 Advantages of the multi-state approach . . . . . . . . . . . . 815.2 Limitations of our study . . . . . . . . . . . . . . . . . . . . 825.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 83viiTable of ContentsBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83AppendicesA Appendix for Chapter 1 . . . . . . . . . . . . . . . . . . . . . . 92B Appendix for Chapter 2 . . . . . . . . . . . . . . . . . . . . . . 93B.1 Survival models . . . . . . . . . . . . . . . . . . . . . . . . . 93B.1.1 Candidate models . . . . . . . . . . . . . . . . . . . . 93B.1.2 Final model for primary exposure measure . . . . . . 94B.1.3 Final model for secondary exposure measure . . . . . 98B.2 Cox model results using the cohort from Ng et al. . . . . . . 100B.3 Time-to-exposure model . . . . . . . . . . . . . . . . . . . . 101B.3.1 Functional form . . . . . . . . . . . . . . . . . . . . . 101B.3.2 Model diagnostics . . . . . . . . . . . . . . . . . . . . 103B.3.3 Time-to-exposure final model results . . . . . . . . . 103C Appendix for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . 105C.1 Transition OH . . . . . . . . . . . . . . . . . . . . . . . . . . 105C.1.1 Model with primary measure of exposure . . . . . . . 105C.1.2 Model with secondary measure of exposure . . . . . . 112C.2 Transition OD . . . . . . . . . . . . . . . . . . . . . . . . . . 116C.2.1 Model with primary measure of exposure . . . . . . . 116C.2.2 Model with secondary measure of exposure . . . . . . 122C.3 Transitions HO . . . . . . . . . . . . . . . . . . . . . . . . . 124C.3.1 Model with primary measure of exposure . . . . . . . 124C.3.2 Model with secondary measure of exposure . . . . . . 135C.4 Transition HD . . . . . . . . . . . . . . . . . . . . . . . . . . 139C.4.1 Model with primary measure of exposure . . . . . . . 139C.4.2 Model with secondary measure of exposure . . . . . . 144D Appendix for Chapter 4 . . . . . . . . . . . . . . . . . . . . . . 147D.1 Pseudo-value Model O for Time Period 1 . . . . . . . . . . . 147D.2 Pseudo-value Model O for Time Period 2 . . . . . . . . . . . 150D.3 Pseudo-value Model O for Time Period 3 . . . . . . . . . . . 153D.4 Pseudo-value Model O for Time Period 4 . . . . . . . . . . . 156D.5 Pseudo-value Model H for Time Period 1 . . . . . . . . . . . 158D.6 Pseudo-value Model H for Time Period 2 . . . . . . . . . . . 161D.7 Pseudo-value Model H for Time Period 3 . . . . . . . . . . . 163viiiTable of ContentsD.8 Pseudo-value Model H for Time Period 4 . . . . . . . . . . . 166D.9 Summary of covariate effects . . . . . . . . . . . . . . . . . . 169ixList of Tables1.1 Characteristics of the cohort . . . . . . . . . . . . . . . 61.2 Annual summary of the study . . . . . . . . . . . . . . 81.3 Summary of DMD usage . . . . . . . . . . . . . . . . . 151.4 Summary of number of DMDs used . . . . . . . . . . . 152.1 HRs of Exposure Measures . . . . . . . . . . . . . . . . 314.1 Coefficients for exposure and interaction terms in pseudo-value regression model O . . . . . . . . . . . . . . . . . 764.2 Coefficients for exposure and interaction terms in pseudo-value regression model H . . . . . . . . . . . . . . . . . 78A.1 Time of first prescriptions . . . . . . . . . . . . . . . . . 92B.1 Unstratified survival model . . . . . . . . . . . . . . . . 93B.2 Survival model with fixed stratification by index period 94B.3 Model with time-dependent stratification by calendarperiod . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94B.4 Final model with primary exposure measure . . . . . . 98B.5 Final model with secondary exposure measure . . . . . 100B.6 Cox model using the cohort from Ng et al. [1] . . . . . 101B.7 Initial time-to-exposure model with time-dependent strat-ification by calendar period . . . . . . . . . . . . . . . . 101B.8 Time-to-exposure model with stratification by calendarperiod and time-dependent coefficient for age . . . . . . 104C.1 Unstratified model . . . . . . . . . . . . . . . . . . . . . 105C.2 Model with fixed stratification by index period . . . . . 106C.3 Model with time-dependent stratification by calendarperiod . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106C.4 Final model for the transition from state O to state H 110xList of TablesC.5 Final model using the secondary exposure measure forthe transition from state O to state H . . . . . . . . . . 116C.6 Unstratified model . . . . . . . . . . . . . . . . . . . . . 117C.7 Model with fixed stratification by index period . . . . . 117C.8 Model with time-dependent stratification by calendarperiod . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118C.9 Final model for the transition from state O to state D 120C.10 Final model using the secondary exposure measure forthe transition from state O to state D . . . . . . . . . . 124C.11 Unstratified model . . . . . . . . . . . . . . . . . . . . . 125C.12 Model with fixed stratification by index period . . . . . 125C.13 Model with time-dependent stratification by calendarperiod . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126C.14 First set of coefficients where t∗ < 1 . . . . . . . . . . . 132C.15 Second set of coefficients where 1 ≤ t∗ ≤ 90 . . . . . . . 132C.16 Third set of coefficients where t∗ > 90 . . . . . . . . . . 133C.17 First set of coefficients where t∗ < 1 . . . . . . . . . . . 137C.18 Second set of coefficients where 1 ≤ t∗ ≤ 90 . . . . . . . 138C.19 Third set of coefficients where t∗ > 90 . . . . . . . . . . 139C.20 Unstratified model . . . . . . . . . . . . . . . . . . . . . 140C.21 Model with fixed stratification by index period . . . . . 140C.22 Model with time-dependent stratification by calendarperiod . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141C.23 Final model for transition from state H to state D . . . 143C.24 Final model for transition from state H to state D . . . 146D.1 Summary of Pseudo-value Model O for Time Period 1 . 147D.2 Summary of Pseudo-value Model O for Time Period 2 . 151D.3 Summary of Pseudo-value Model O for Time Period 3 . 154D.4 Summary of Pseudo-value Model O for Time Period 4 . 156D.5 Summary of Pseudo-value Model H for Time Period 1 159D.6 Summary of Pseudo-value Model H for Time Period 2 161D.7 Summary of Pseudo-value Model H for Time Period 3 164D.8 Summary of Pseudo-value Model H for Time Period 4 167xiList of Figures1.1 A histogram of days of follow-up . . . . . . . . . . . . . 61.2 Histogram of number of treatment switches amongstindividuals that used a DMD at least once . . . . . . . 101.3 A histogram of days until exposure for the cohort . . . 162.1 Estimated survival for the cohort over the duration ofthe study. . . . . . . . . . . . . . . . . . . . . . . . . . . 252.2 Estimated survival given comorbidity at the index date:comorbidity of 0 = black, comorbidity of 1 = red, co-morbidity of 2 = blue. . . . . . . . . . . . . . . . . . . . 262.3 Estimated survival based upon exposure at any pointduring follow-up: never exposed = red, exposed at somepoint during follow-up = blue. . . . . . . . . . . . . . . 262.4 Estimated survival stratified by socioeconomic statusquintile: 1 = blue, 2 = red, 3 = black, 4 = green, 5 =purple. . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.5 Estimated survival stratified by pre-index hospitaliza-tions: zero = blue, one or more = red. . . . . . . . . . 272.6 Estimated survival stratified by age at index date: youngerthan 44 years old = blue, 44 or older = red. . . . . . . 282.7 Average pre-index hospitalizations in exposed individ-uals over time . . . . . . . . . . . . . . . . . . . . . . . 342.8 Weighted cumulative exposure function over 10-yearwindow for hazard of death over a 10-year window . . . 382.9 Weighted cumulative exposure function over 10-yearwindow for hazard of death over a 10-year window dif-ferent numbers of knots . . . . . . . . . . . . . . . . . . 383.1 Illness-Death model: O = out of hospital, H = hospi-talized, D = death . . . . . . . . . . . . . . . . . . . . . 423.2 Cumulative intensity of transitions out of state O . . . 523.3 Cumulative intensity of transitions out of state H . . . 52xiiList of Figures3.4 The probability of being hospitalized at time t. Toreduce noise in the AJ estimate and CIs, the plot showsa moving average with a window of 250 days. . . . . . . 533.5 The probability of dying at time t. To reduce noisein the AJ estimate and CIs, the plot shows a movingaverage with a window of 250 days. . . . . . . . . . . . 533.6 The probability of hospitalization by time t stratifiedby baseline comorbidity: comorbidity of 0 = black, co-morbidity of 1 = red, comorbidity of 2 = blue. To re-duce noise in the AJ estimates, the plot shows a movingaverage with a window of 250 days. . . . . . . . . . . . 54B.1 Age at index date modelled with a natural spline versusa linear fit in red. . . . . . . . . . . . . . . . . . . . . . 95B.2 Previous hospitalizations modelled with a natural splineversus a log(x+ 1) transformation in red. . . . . . . . . 95B.3 Scaled Schoenfeld residuals for the coefficient of the pri-mary exposure. Model estimate in red. . . . . . . . . . 96B.4 Deviance residuals for the final model using the primaryexposure measure . . . . . . . . . . . . . . . . . . . . . 97B.5 Scaled Schoenfeld residuals for the coefficient of lowexposure. . . . . . . . . . . . . . . . . . . . . . . . . . . 98B.6 Scaled Schoenfeld residuals for the coefficient of highexposure. . . . . . . . . . . . . . . . . . . . . . . . . . . 99B.7 Previous hospitalizations modelled with a natural splineversus a log(x+ 1) transformation in red. . . . . . . . . 102B.8 Age at index date modelled with a natural spline versusa linear fit in red. . . . . . . . . . . . . . . . . . . . . . 102B.9 Scaled Schoenfeld residuals with time on a log scale forthe coefficient of age. . . . . . . . . . . . . . . . . . . . 103C.1 Previous hospitalizations modelled with a natural spline(5 d.f.) versus log(x+ 1) in red. . . . . . . . . . . . . . 107C.2 Age at index date modelled with a natural spline (3d.f.) versus a linear fit in red. . . . . . . . . . . . . . . 107C.3 Scaled Schoenfeld residuals for the coefficient of exposure.108C.4 Scaled Schoenfeld residuals for the coefficient of pre-index hospitalizations. . . . . . . . . . . . . . . . . . . . 109C.5 Scaled Schoenfeld residuals for the coefficients of co-morbidity. . . . . . . . . . . . . . . . . . . . . . . . . . 109xiiiList of FiguresC.6 Deviance residuals for final O → H model . . . . . . . 111C.7 Scaled Schoenfeld residuals for the coefficient of lowexposure. . . . . . . . . . . . . . . . . . . . . . . . . . . 113C.8 Scaled Schoenfeld residuals for the coefficient of highexposure. . . . . . . . . . . . . . . . . . . . . . . . . . . 113C.9 Scaled Schoenfeld residuals for the coefficient of sex. . . 114C.10 Scaled Schoenfeld residuals for the coefficient of thefirst basis function of the spline for age at the indexdate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114C.11 Scaled Schoenfeld residuals for the coefficient of SES5. 115C.12 Previous hospitalizations modelled with a natural spline(5 d.f.) versus log(x+ 1) in red. . . . . . . . . . . . . . 118C.13 Age at index date modelled with a natural spline (5d.f.) versus a linear fit in red. . . . . . . . . . . . . . . 119C.14 Scaled Schoenfeld residuals for the coefficient of exposure.119C.15 Scaled Schoenfeld residuals for the coefficient of pre-index hospitalizations. . . . . . . . . . . . . . . . . . . . 120C.16 Deviance residuals for the final O → D model . . . . . 121C.17 Scaled Schoenfeld residuals for the coefficient of lowexposure. . . . . . . . . . . . . . . . . . . . . . . . . . . 122C.18 Scaled Schoenfeld residuals for the coefficient of highexposure. . . . . . . . . . . . . . . . . . . . . . . . . . . 123C.19 Scaled Schoenfeld residuals for the coefficient of pre-index hospitalizations. . . . . . . . . . . . . . . . . . . . 123C.20 Previous hospitalizations modelled with a natural spline(5 d.f.) versus log(x+ 1) in red. . . . . . . . . . . . . . 127C.21 Age at admission date modelled with a natural spline(3 d.f.) versus a linear fit in red. . . . . . . . . . . . . . 127C.22 Scaled Schoenfeld residuals for the coefficient of exposure.128C.23 Scaled Schoenfeld residuals for the coefficients of co-morbidity. . . . . . . . . . . . . . . . . . . . . . . . . . 129C.24 Scaled Schoenfeld residuals for the coefficient of pre-index hospitalizations. . . . . . . . . . . . . . . . . . . . 129C.25 Scaled Schoenfeld residuals for the coefficient of SES5. 130C.26 Scaled Schoenfeld residuals for the coefficient of sex. . . 130C.27 Scaled Schoenfeld residuals for the coefficient of thethird basis function for the spline modelling age at ad-mission. . . . . . . . . . . . . . . . . . . . . . . . . . . . 131C.28 Deviance residuals for the final H → O model . . . . . 134xivList of FiguresC.29 Scaled Schoenfeld residuals for the coefficient of lowexposure. . . . . . . . . . . . . . . . . . . . . . . . . . . 135C.30 Scaled Schoenfeld residuals for the coefficient of highexposure . . . . . . . . . . . . . . . . . . . . . . . . . . 136C.31 Scaled Schoenfeld residuals for the coefficient of thethird basis function for the spline modelling age at ad-mission. . . . . . . . . . . . . . . . . . . . . . . . . . . . 136C.32 Previous hospitalizations modelled with a natural spline(5 d.f.) versus log(x+ 1) in red. . . . . . . . . . . . . . 141C.33 Age at admission modelled with a natural spline (5 d.f.)versus a linear fit in red. . . . . . . . . . . . . . . . . . 142C.34 Scaled Schoenfeld residuals for the coefficient of exposure.142C.35 Deviance residuals for the final H → D model . . . . . 144C.36 Scaled Schoenfeld residuals for the coefficient of lowexposure . . . . . . . . . . . . . . . . . . . . . . . . . . 145C.37 Scaled Schoenfeld residuals for the coefficient of highexposure . . . . . . . . . . . . . . . . . . . . . . . . . . 145D.1 Age at index date modelled with a natural spline . . . 148D.2 Residual plot for Model O for time period 1 . . . . . . 149D.3 Residual plot for Model O for time period 1, smallest1% of fitted values removed . . . . . . . . . . . . . . . . 149D.4 Histogram of fitted pseudo-values for Model O for timeperiod 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 150D.5 Residual plots for covariates in Model O for Time Pe-riod 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150D.6 Age at index date modelled with a natural spline forModel O for Time Period 2 . . . . . . . . . . . . . . . . 152D.7 Residual plot for pseudo-value Model O for Time Pe-riod 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152D.8 Residual plot for pseudo-value Model O for Time Pe-riod 2 with smallest 1% of fitted values removed . . . . 153D.9 Residual plots for covariates in Model O for Time Pe-riod 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153D.10 Age at index date for Model O for Time Period 3 . . . 154D.11 Residual plot for pseudo-value Model O for Time Pe-riod 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155D.12 Residual plot for pseudo-value Model O for Time Pe-riod 3 with smallest 1% of fitted values removed . . . . 155xvList of FiguresD.13 Residual plots for covariates in Model O for Time Pe-riod 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156D.14 Age at index date for Model O for Time Period 4 . . . 157D.15 Residual plot for pseudo-value Model O for Time Pe-riod 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158D.16 Residual plots for covariates in Model O for Time Pe-riod 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158D.17 Age at index date for Model H for Time Period 1 . . . 160D.18 Residual plot for pseudo-value Model H for Time Pe-riod 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160D.19 Residual plots for covariates in Model H for Time Pe-riod 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161D.20 Age at index date for Model H for Time Period 2 . . . 162D.21 Residual plot for pseudo-value Model H for Time Pe-riod 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163D.22 Residual plots for covariates in Model H for Time Pe-riod 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163D.23 Age at index date for Model H for Time Period 3 . . . 165D.24 Residual plot for pseudo-value Model H for Time Pe-riod 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165D.25 Residuals against age for pseudo-value Model H forTime Period 3 . . . . . . . . . . . . . . . . . . . . . . . 166D.26 Residual plots for pre-index hospitalization in Model Hfor Time Period 3 . . . . . . . . . . . . . . . . . . . . . 166D.27 Age at index date for Model H for Time Period 4 . . . 168D.28 Residual plot for pseudo-value Model H for Time Pe-riod 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168D.29 Residual plots for covariates in Model H for Time Pe-riod 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169xviAcknowledgementsUBC’s Department of Statistics has been an incredible place to learn for thepast two years. I would like to thank all of my professors for their enthu-siasm and eagerness to share their knowledge. I have learned a great dealhere. Specifically, I want to thank Dr. Lang Wu, who was a great teacherin two courses, for his help as an external reviewer for this thesis.But most of all, I would like to thank my two supervisors on this thesis, Drs.John Petkau and Yinshan Zhao, for their generosity. Before this thesis, Iknew little about survival analysis, and their patience during my researchassistantship, even as I proved during our meetings that there is in fact sucha thing as a stupid question, was crucial. Their expertise and feedback wasinvaluable to not just my research, which benefitted greatly from their di-rection, but also my development as a statistician. Under their guidance, Ilearned how to think critically about statistics, approach a research questionproperly, interrogate the data, and communicate my results. This thesis isas much a result of their efforts as it is mine. But of course, any errors,should they exist, are mine.This thesis would also not have been possible without the Pharmacoepidemi-ology in Multiple Sclerosis research lab. Not only does my research extendtheir work but support from the PiMS team was vital, from their early ad-vice on beginning my project to answering any questions I had throughoutthe process.xviiChapter 1IntroductionMultiple sclerosis (MS) is a chronic, demyelinating, autoimmune disease[3], [4]. It is chronic because the course of the disease is slow and maytake decades to progress. Demyelinating refers to what causes the symp-toms. Myelin covers and protects nerves in the human body, and diseasesin which myelin is degraded, such as MS, are demyelinating diseases. De-graded myelin leads to nerve signals being disrupted which manifests as,among other symptoms, the loss of vision and motor function (e.g., lack ofcoordination, weakness), cognitive decline, mental illness, such as depres-sion, and fatigue. MS is an autoimmune disease because the body’s ownimmune system attacks and degrades the myelin in the central nervous sys-tem.MS is treated with disease-modifying drugs (DMDs) [4], [5]. DMDs havedifferent mechanisms and are administered differently (e.g., orally, intra-venously), but they aim to decrease the immune response that is thought tocause demyelination. The efficacy of DMDs was established in clinical trialsthat were of short duration and excluded older individuals and individualsliving with comorbidities [6]–[8]. This is common for clinical trials but hasled to a lack of knowledge of the effects of DMDs over a long time—which isespecially problematic for chronic diseases—and the effects of DMDs on in-dividuals that do not meet the traditional eligibility criteria for clinical trials.The motivation of this thesis is to use health administrative data collectedfrom individuals residing in British Columbia over 22 years to remedy thisgap by investigating the following research questions:1. Are DMDs associated with a decrease in the frequency of hospitaliza-tion or shorter hospital visits?2. Are DMDs associated with a decrease in mortality?3. Does an individual’s comorbidity burden change either of the precedingassociations?11.1. Study cohort4. How, if at all, do the effects of DMDs change over time?This thesis extends and builds upon previous work in the “Disease-modifyingDrug Safety and Effectiveness in Multiple Sclerosis” (DRUMS) study. Twopapers by Ng et al. investigated, among other outcomes, the associationbetween DMDs and mortality [1] and the association between DMDs andhospitalization [2]. The studies by Ng et al. investigated mortality and hos-pitalizations separately, whereas this thesis aims to study hospitalizationrates, length of hospitalizations, and mortality simultaneously using multi-state modelling.The remainder of Chapter 1 describes the cohort in this retrospective studyand discusses exposure to disease-modifying drugs. Chapter 2 uses time-to-event models to investigate the association between DMDs and mortalityand what factors dictate when and if someone becomes exposed. Chapter2 also lays the foundations for the introduction of multi-state modelling inChapter 3. In addition to outlining the theoretical foundation of multi-statemodels, Chapter 3 applies multi-state methods to the health administrativedata of individuals living with MS to answer the research questions. Theeffect of disease-modifying drugs on the time spent in and out of the hos-pital is estimated in Chapter 4 using pseudo-value regression. This thesisconcludes with Chapter 5 which summarizes and discusses the findings, anylimitations, and potential directions for future work.1.1 Study cohortThe study cohort was composed of 9,124 individuals whose data was linkedfrom several provincial datasets from Population Data BC (PopData): Phar-maNet recorded prescription information; the Discharge Abstract Databasecontained hospitalizations coded using the International Classification ofDiseases (ICD) system; Medical Services Plan (MSP) Payment InformationFile captured physician claims; BC Vital Events and Statistics - Deathsrecorded the date of death for individuals in the study cohort; informationon BC residency and covariates such as sex, date of birth, and socioeco-nomic status was pulled from the Central Demographics File. Data fromthe Discharge Abstract Database, Medical Services Plan, and PharmaNetwere utilized to identify individuals with MS and identify comorbidities.To be included in this study, individuals must:21.1. Study cohort• be living with multiple sclerosis, as identified by a validated algorithm:≥ 3 MS diagnostic codes (ICD-9/10 340/G35) or ≥ 1 dispensation ofa DMD [9].• experience their first demyelinating or MS-related event on or afterJanuary 1, 1996. The date of this event is referred to henceforth asthe individual’s index date. The index date is assumed to be the timethat multiple sclerosis clinically manifested itself. One corollary isthat it is assumed no individuals have received treatment for multiplesclerosis (such as disease-modifying drugs) before their index date.• be at least 18 years of age at the index date.• be a resident of BC (and thus registered in the BC healthcare database)for more than five years (> 365×5 days) prior to the index date. Thisreduces the chances—but is no guarantee—that the first demyelinatingevent in the database was not the first such event that the personexperienced.Follow-up began on the index date and ended on the earliest of death, emi-gration from BC, dispensation of ocrelizumab or daclizumab, or December31, 2017, the administrative censoring date. Subjects were censored uponocrelizumab and daclizumab dispensation since fewer than six individualswere ever prescribed either drug, and daclizumab ceased to be prescribedsoon after the end of the study due to safety reasons.Five covariates were considered for the analysis:• Age at index date (unless otherwise stated)• Sex• Socioeconomic status (SES) at index dateSES was categorized by income quintile where one was the lowest quin-tile and five the highest. Income was estimated from neighborhood-level income via postal codes. SES during the year of the index datewas not always available so, when required, the closest value within3 years to the index date was used. Twenty-six people satisfied theinclusion criteria but had missing income data and therefore were ex-cluded.31.1. Study cohort• ComorbidityConsistent with how comorbidity was determined in the two previousDRUMS studies [1], [2], comorbidity is determined by the CharlsonComorbidity Index (CCI), excluding hemiplegia and paraplegia, usinghealth records in the previous year (365 person-days). A comorbidityof zero corresponds to a CCI of zero (no comorbidities), a comorbidityof one to a CCI of one, and a comorbidity of two to a CCI ≥ 2. Indi-viduals in this final category face complicated MS treatment decisions.They either have one serious comorbidity that itself has a CCI of atleast 2, such as AIDS, cirrhosis, or cancer, or multimorbidity (two ormore comorbidities) [8], [10].In Chapter 2, comorbidity is measured at the index date and fixed.In Chapter 3, comorbidity is treated as a time-dependent covariatethat is updated annually on an individual’s birthday, if an individualis not hospitalized. If an individual is hospitalized on their birthday,their comorbidity is not updated until they are discharged. This wasdone so that comorbidity is constant during hospitalization. Once dis-charged, comorbidity is immediately updated to the value that wascalculated on the person’s most recent birthday. For example, if anindividual’s birthday is on January 10th but they are hospitalized fromJanuary 9th, 2010 until January 15th, 2010, their comorbidity duringthat hospitalization was last updated on January 10th, 2009. Upondischarge on January 15th, their comorbidity was updated to the valuecalculated on January 10th, 2010. In summary, after the baseline mea-surement, it is assumed that an individual’s comorbidity when theyare x years old depends on events that occurred when they were x− 1years old.In Chapter 4, pseudo-value regression necessitates only fixed covari-ates, and the effect of comorbidity is assumed to only depend on itsvalue at certain landmark times. This is elaborated upon in Chapter4.• The number of hospitalizations in the five years prior to the index dateAdjusting for the number of previous hospitalizations may capturemore information about a person’s health that is not reflected in co-41.1. Study cohortmorbidity.Table 1.1 shows a summary of the characteristics of the entire cohort (n =9124). More than ten percent (10.3%) of individuals that satisfied the in-clusion criteria died during follow-up, and approximately one-third (31.8%)of study participants had zero hospitalizations post-index.This study collects data from 1996 to 2017. There was 94,578.5 total person-years. The maximum follow-up time was 22.0 years while the median lengthof follow-up was 9.9 years. Figure 1.1 shows a histogram of the number ofdays of follow-up for the cohort.As stated above, clinical trials for disease-modifying drugs are compara-tively short (typically ≤ 3 years) and often exclude individuals over 50 andindividuals with comorbidities [6]–[8]. In our study, 1,658 (18.2%) and 847(9.3%) of subjects had a baseline comorbidity of one and two respectively.When comorbidity is used as a time-dependent covariate, there are 15,029person-years of follow-up on individuals while they have a comorbidity ofone and 8,446 person-years of follow-up on individuals while they have acomorbidity of two. The mean age at the index date was 44.8; 33.1% of thecohort were 50 years or older (27,453.7 person-years of follow-up) and 13.4%were 60 years or older (9,304.3 person-years of follow-up).51.1. Study cohortTable 1.1: Characteristics of the cohortCharacteristicNumber of individuals 9,124Age at index date — mean (s.d.) — years 44.75 (13.43)Female sex — # (%) 6,668 (73.08)SES at index date — # (%)1 1,635 (17.92)2 1,792 (19.31)3 1,909 (20.92)4 1,997 (21.89)5 1,821 (19.96)Comorbidity at index date — # (%)0 6,608 (72.42)1 1,658 (18.17)2 858 (9.40)Number of previous hospitalizations — mean (s.d.) 1.39 (2.45)Died — # (%) 965 (10.28)Hospitalized at least once during follow-up — # (%) 6,225 (68.23)Ever prescribed a DMD — # (%) 2,424 (26.57)Figure 1.1: A histogram of days of follow-upExploratory analysis of the hospitalization patterns revealed two potentialissues. Of the 25,792 hospital visits experienced by the cohort where the pa-tient was discharged alive, 13,825 (53.6%) were shorter than one day. These61.1. Study cohortsame day discharges may be fundamentally different than visits longer thanone day and require consideration in the multi-state model in Chapter 3.Another issue arises for the opposite reason: some hospitalizations are un-expectedly long. The median length of a hospital visit was half a day or,excluding same day discharges, four days. But the longest visit in the datasetwas 2,394 days and, for visits that did not result in a same day discharge,7.3% were greater than 30 days. Lengthy hospital visits in the databasemay not be hospitalizations at all but some kind of long-term care that mayhave an outsized influence on model estimates.Table 1.2 provides a yearly breakdown of the study. The ratio of hospitalvisits shorter than one day to longer hospitalizations appeared to increaseafter 2002 compared to the first seven years of the study. 277 people werecensored prior to the end of the study, and 7,882 people were censored onthe last day of the study.71.1.StudycohortTable 1.2: Annual summary of the studyYear Entered Died CensoredAt risk onDec 31Same daydischargesOtherhospitalizationsRatio ofsame day to otherhospitalizations1996 420 < 6 < 6 419 51 116 0.441997 388 < 6 < 6 802 74 187 0.401998 425 8 < 6 1,215 110 295 0.371999 425 13 7 1,620 157 370 0.422000 446 9 8 2,049 212 493 0.432001 430 21 9 2,449 274 573 0.482002 468 21 11 2,885 301 613 0.492003 431 31 8 3,277 397 733 0.542004 397 29 9 3,636 422 809 0.522005 431 43 13 4,011 498 928 0.542006 460 50 13 4,408 543 1044 0.522007 406 50 17 4,747 636 1172 0.542008 440 44 10 5,133 709 1301 0.542009 423 71 9 5,476 775 1420 0.552010 453 48 21 5,860 850 1524 0.562011 411 60 16 6,195 975 1728 0.562012 435 60 25 6,545 1021 1881 0.542013 392 63 20 6,854 1097 1996 0.552014 456 85 18 7,207 1210 2148 0.562015 418 91 16 7,518 1302 2405 0.542016 333 110 29 7,712 1341 2455 0.552017 236 54 12* 7,882 870 1601 0.54*individuals administratively censored on December 31, 2017 are not included81.2. Exposure to disease-modifying drugs1.2 Exposure to disease-modifying drugsThis section discusses the measures of exposure upon which conclusions ofthis thesis are based. An active prescription for a drug does not mean anindividual is exposed to it. Ideally, an individual is “exposed” to a treatmentif they experience the effects from that treatment and “unexposed” if theydo not, but determining when an individual is affected by a treatment is notsimple. Prescription data on residents in BC is collected via PharmaNetwhich tracks when and how much of a prescription drug was dispensed incommunity-based pharmacies. The measures of exposure outlined in thissection are determined solely from PharmaNet data so other useful informa-tion on prescription drug use, such as adherence, which may be collected indoctor visits is unavailable. Disease-modifying drugs often take months ofcontinuous use to have a clinical effect which may continue for years afterdiscontinuation.General issues with defining a measure of exposure in observational studiessuch as treatment switching and more specific problems that are prevalent inthis study such as possible reverse causality are described below. Exposureis defined in Section 1.2.2 and discussed in the subsequent sections. Anexploratory analysis at the end of this section summarizes exposure withinthe cohort.1.2.1 A note about treatment switchingTreatment switching occurred if treatment was discontinued or initiated.Treatment was discontinued if there were no dispensations for the samedrug within 30 days of the supply running out or if a different drug wasprescribed (a new treatment was initiated) before the supply of the old drugwas finished. Gaps in supply ≤ 30 days are ignored. Treatment was initiatedwhen an individual filled a prescription after not using any drug for morethan 30 days.Most of the 2,424 DMD users changed treatment. A histogram of the num-ber of changes in treatment amongst DMD users in the cohort is shown inFigure 1.2. The maximum number of changes in treatment was 47 but themedian number of changes amongst DMD users was three. By definition, anindividual that uses a disease-modifying drug must have at least one treat-ment switch since no individual was prescribed a DMD at their index date.After being prescribed a DMD, 318 people did not stop using DMDs, and of91.2. Exposure to disease-modifying drugsthis subset, 261 people stayed on the same DMD from their first prescrip-tion until the end of their follow-up. The frequency of treatment switchingis also evident in Table 1.3 and Table 1.4.Figure 1.2: Histogram of number of treatment switches amongst individualsthat used a DMD at least onceDeciding upon the measures of exposure when there is treatment switchingis difficult since study subjects may be prescribed several drugs over thecourse of their follow-up. This is further complicated by the variability inthe number of times individuals changed treatment—the exposure measuremust be applicable to the individual who changed treatment 47 times aswell as the individual who changed only once.1.2.2 Exposure measuresDuring the study, nine disease-modifying drugs were prescribed to individ-uals living with MS in BC. They became available at different times (thedates that they were first prescribed are listed in Appendix A). As the studywent on, there were more treatment options for individuals living with MS.In this thesis, days supplied and days prescribed are synonymous. Fewerthan six individuals were ever prescribed either daclizumab (which was with-drawn from the market shortly after the end of the study) or ocrelizumab.The first prescriptions of daclizumab and ocrelizumab occurred within the101.2. Exposure to disease-modifying drugsfinal year of the study. Individuals were censored on the first day of theirdaclizumab or ocrelizumab prescription.The primary measure of exposure is based on the minimum time thoughtto produce a clinical effect [11]. It was assumed that individuals adheredcompletely to drugs that were prescribed and took only one DMD at a time.If a drug prescription was renewed before the old prescription ended, it wasassumed that individuals stopped using the old supply and began using thenew supply.An individual is considered as exposed for the remainder of follow-up oncethey fulfill any one of the following three criteria:i. Beta-interferon, glatiramer acetate: ≥ 180 days of use during follow-upii. Dimethyl fumarate, fingolimod, natalizumab, teriflunomide: ≥ 90days of use during follow-upiii. Alemtuzumab: 90 days after prescription dateExposure is time-dependent in this study. To illustrate this, consider a hy-pothetical individual that entered the study on January 1, 2010 and wascensored on December 31, 2017. Like all subjects, this individual is “un-exposed” at study entry. Suppose that this person was prescribed beta-interferon on February 1, 2010 to last for two months until April 1, 2010.This does not meet the threshold of (i.) above and this individual is still“unexposed”. Next, suppose that this same individual was prescribed glati-ramer acetate for 10 months beginning on June 1, 2010. This individual’sexposure would be “unexposed” until November 26, 2010 (179 days of useincluding June 1, 2010). As of November 27, the 180th day of use, thisindividual is exposed until the end of follow-up.A secondary measure of exposure was also considered based upon the num-ber of days supplied of each DMD. Individuals were categorized into oneof three exposure levels: no exposure if individuals did not satisfy any ofthe exposure definitions listed above, low exposure if individuals satisfied atleast one of the exposure definitions above but did not have high exposure,and high exposure if they satisfied at least one of the following definitions:i. Beta-interferon: ≥ 4.5 years of useii. Glatiramer acetate: ≥ 3.5 years of use111.2. Exposure to disease-modifying drugsiii. Dimethyl fumarate: ≥ 1.5 years of useiv. Fingolimod: ≥ 2 years of usev. Natalizumab: ≥ 2.5 years of usevi. Teriflunomide: ≥ 1.25 years of usevii. Alemtuzumab: 2 yearly courses or ≥ 2 years of useCategories were defined based upon the distribution of DMD exposure inthe cohort. The thresholds were rounded for convenience, and the person-years of follow-up for low exposure was approximately three times that ofthe follow-up for high exposure. The one exception to this is alemtuzumabwhere the ratio of low to high exposure follow-up was 2.3:1. Alemtuzumabis prescribed for a yearly course, and the first prescription for alemtuzumabin the cohort was in 2014. Individuals exposed to alemtuzumab receivedeither 1, 2, or 3 courses, which restricted the definition of high exposure—athreshold of only 1 year included too many people and any threshold greaterthan 2.3 years included too few. Most but not all individuals exposed toalemtuzumab were also exposed to another DMD beforehand.To illustrate the secondary measure of exposure, consider the same hypo-thetical individual as above but suppose that this individual was prescribedglatiramer acetate for 5 years continuously beginning on June 1, 2010. Thisindividual would be “unexposed” from the index date until November 27,2010, have “low” exposure beginning on November 28, 2010 until November28th, 2013. From November 29th, 2013 (3.5 × 365 days of use includingJune 1, 2010) until December 31, 2017, this individual’s exposure is “high”.1.2.3 Rationale for using progressive measures of exposureAn exposure is progressive if an individual cannot return to a lower levelof exposure once they have satisfied the definition of a higher level. Forthe primary measure, an exposed person cannot become unexposed, andaccording to the secondary measure, a highly exposed person cannot returnto a state of low exposure and a person with low exposure cannot becomeunexposed. A progressive measure of exposure has three important benefits.First, once a DMD has had a clinical effect on an individual, the courseof their disease may be fundamentally different from that of a person whois unaffected. Recent reviews on studies of DMD discontinuation reported121.2. Exposure to disease-modifying drugsthat there was a limited or negligible risk of relapse or disability progressionafter discontinuing DMD treatment. One of these review articles [12] re-ported this result in patients that discontinued DMDs after two years of usebut another systematic review and meta-analysis [13] demonstrated that thelonger an individual was treated prior to discontinuation, the more likely itwas that their disease would not worsen once they stopped treatment. Thissuggests that an individual can still experience clinical benefits after discon-tinuing DMD treatment but a non-progressive definition may label such aperson as unexposed. Non-progressive definitions also require a decision onhow long after discontinuation someone is still considered exposed. Such achoice is made under assumptions that may not hold true.Second, short-lived treatment switches in response to the health of the in-dividual are not reflected to the same extent in the exposure variable. Pro-gressive measures reduce the “noise” of treatment switching. An individual’sexposure depends upon the sustained treatments that they have undertaken,as defined by clinical guidelines.Third, the exposure measure may be defined in such a way that it reducesthe impact reverse causality has on the estimates, if it exists. In lieu of arandomized controlled trial to test the efficacy of some drug interventionand establish a causal relationship, an observational study seeks to establishan association between the drug and some outcome. Confounders are con-trolled for in order to strengthen the findings of the study and reduce thechances that some factor other than the drug is responsible for the associ-ation. Reverse causality occurs when, instead of the intervention affectingthe outcome, the association exists because the outcome has an effect onwhether the intervention is prescribed or not. This is another potential is-sue that may arise in observational studies and weaken the legitimacy oftheir findings.In the context of a study on the association between exposure to disease-modifying drugs and mortality, reverse causality occurs if the associationbetween DMD exposure and mortality is not because DMDs reduce (or in-crease) the chance of death but because individuals at a higher risk of deathare less likely (or more likely) to be prescribed DMDs. Older individuals aremore likely to have life-threatening comorbidities, experience lower DMDefficacy, and stop DMD treatment [14]. In a systematic review of 151 stud-ies on DMDs, Giovannoni, Southam, and Waubant reported that adverseeffects were the most common reason for stopping treatment with DMDs,131.2. Exposure to disease-modifying drugswhich may contribute to reverse causality [15]. Thus, no DMD use may beassociated with old age, frailty, and comorbidities which is in turn associatedwith a greater chance of dying.Three hypothetical examples of this in a study on mortality are: (1) if anindividual experiences a downturn in their health, they may not be ableto tolerate DMDs anymore and cease using them; (2) a clinician may stoptreating a patient with a DMD if there are signs that they may die soon;(3) an older individual with MS, who is at a higher risk of dying, may ex-perience fewer benefits from DMDs and stop treatment. These situationswould also introduce reverse causality into an analysis on all-cause hospital-izations since an older individual with worsening health is more likely to behospitalized and, in these hypotheticals, less likely to be exposed to DMDs.A common strategy in pharmacoepidemiology studies to avoid bias from re-verse causality is to use a lag time [16]. This decouples exposure from recentchanges in treatment which could have been made in response to some un-derlying sign of failing health. The exposure measures for this study are notlagged but nonetheless are defined to be resistant to any vagaries of treat-ment choice. It is expected that this will similarly reduce reverse causalitybias.The presence of reverse causality in this study is explored in Section 2.5.1.2.4 Description of exposureJust over one-quarter of the 9,124 people in the cohort (2,424) had a disease-modifying drug dispensed to them during the study and 2,134 met at leastone definition of exposure. Table 1.3 breaks down exposure by individualDMD and Table 1.4 shows the number of individuals prescribed and ex-posed to between one and five DMDs. More individuals were exposed toonly one DMD than were dispensed only one DMD. This is because theformer includes individuals prescribed several DMDs but who only reachedthe definition of exposure for one of them.There was 18,017 person-years of follow-up while exposed (19% of totalfollow-up), and among individuals that were exposed to at least one DMD,the median number of years exposed was 8.0 years. Individuals that had acomorbidity of one were exposed for 2,213 person-years and individuals thathad a comorbidity of two were exposed for 959 person-years.141.2. Exposure to disease-modifying drugsTable 1.3: Summary of DMD usageTreatment Dispensed Exposed Highly exposedAny disease-modifying drugs 2,428 2,134 1,185Beta-interferon 1,511 1,236 517Glatiramer acetate 965 638 261Alemtuzumab 121 116 37Dimethyl fumarate 474 387 220Fingolimod 222 204 105Natlizumab 151 133 61Terifluomide 334 277 152The counts in the dispensed, exposed, and highly exposed columns are re-spectively the number of individuals that: had ≥ 1 days of supply, satis-fied the primary definition of exposure at some point during follow-up, andreached the threshold of high exposure at some point during follow-up.Table 1.4: Summary of number of DMDs usedNumber of DMDs Dispensed Exposed1 1,426 1,4522 726 5723 210 1204 55 <255 11 <25The counts in the dispensed and exposed columns are respectively the num-ber of individuals that had ≥ 1 days of supply but were not necessarilyexposed and that satisfied the primary definition of exposure at some pointduring follow-up.Figure 1.3 shows the distribution of time-to-exposure (measured since theindex date) for the 2,134 individuals who were exposed during follow-up.Users of disease-modifying drugs were treated with and exposed to anywherefrom one to five different drugs over follow-up.151.2. Exposure to disease-modifying drugsFigure 1.3: A histogram of days until exposure for the cohort16Chapter 2Time-to-event modelsTime-to-event models study how covariates may impact the time until in-dividuals experience some event. Section 2.1 discusses methods used inintensity-based time-to-event models. Section 2.1.1 introduces intensityfunctions which are the foundation of many time-to-event models. Cox re-gression is a popular time-to-event model that is discussed in Section 2.1.2along with inference using the model, useful extensions, and residuals formodel assessment.Section 2.2 presents Kaplan-Meier curves to illustrate the broad trends ofmortality in the cohort. The next sections are three different applications ofCox models. Because Cox models are most often applied to time-to-deathdata, we first implement a Cox model to examine the association betweenDMD exposure and mortality in Section 2.3. Then Section 2.4 applies thesame technique but in a time-to-exposure context to determine the factorsassociated with DMD exposure. The most complex Cox model is left untilthe last section. Section 2.5 extends the Cox model using a technique calledweighted cumulative exposure to explore potential reverse causality.2.1 Intensity-based modelling2.1.1 Intensity functionsThree processes are required to introduce intensity-based models [17]: theevent, covariate, and censoring processes. For the event of interest, let N(t)be the number of events that have occurred up to and including time t.Then, {N(t), t ≥ 0} is the event process. Let X(t) denote a vector ofpossibly time-dependent covariates at time t. Then, {X(t), t ≥ 0} is thecovariate process. Let Y (t) denote the status at time t where Y (t) = 1 if anindividual is not censored by time t and 0 otherwise. Then {Y (t), t ≥ 0} isthe censoring process.Intensities are conditioned on the past. The history up to and including time172.1. Intensity-based modellingt is denoted as H(t) = {N(s), X(s), Y (s), 0 ≤ s ≤ t}. For small enough∆t, it is assumed that at most one event can occur in an interval [t, t+∆t).The intensity function at time t, which is denoted λ(t), is defined in termsof ∆N(t), the number of events that occurred in the interval [t, t+∆t), as:λ(t) = lim∆t→0P (∆N(t) = 1 | H(t))∆t. (2.1)Assume for simplicity that λ(t) = 0 when Y (t) = 0. The intensity may beinterpreted as the instantaneous risk of the event occurring at time t. Theprobability of an event occurring in the interval [t, t+∆t) is:P (∆N(t) = 1 | H(t)) = λ(t)∆t+ o(∆t) . (2.2)2.1.2 Cox regressionTime-to-event models are specified with intensity functions [17], and theCox regression model is one example. Cox regression is a semi-parametrictime-to-event model that takes the form:λ(t) = λ0(t) exp(x(t)′β), (2.3)where λ0(t) is the baseline intensity function, β is a vector of coefficients,and x(t) is a vector of possibly time-dependent covariates.Cox regression assumes the intensity depends on H(t) only through x(t),which is derived from {X(s), 0 ≤ s ≤ t}. Covariates may increase or de-crease the intensity function. No assumptions are made about the shape ofthe baseline intensity function.The Cox model specified by (2.3) assumes every person shares the samebaseline hazard function. In many contexts, this assumption may not berealistic. A stratified Cox model relaxes this assumption. Individuals maybe in only one stratum at any given time. Individuals in the same stratahave the same baseline hazard function but this hazard is allowed to differbetween strata. If there are K strata, the simplest version of a stratifiedCox model is specified as follows:λk(t) = λk0(t) exp(x(t)′β), where k = 1, ..., K. (2.4)Equation (2.4) assumes the effect of every covariate is constant across strata.In some cases, this may not be the case, and stratum-specific coefficients βk182.1. Intensity-based modellingmay be specified.The DRUMS study is over two decades long. Medical practices, treatmentoptions, and technology changed markedly from 1996 to 2017 with respectto diagnosis and treatment of MS but also generally. Patterns of survivallikely changed because of these advances, and it may be unrealistic to as-sume that individuals followed early in the study have the same baselinehazard function as individuals followed later. It may be necessary to strat-ify based upon when an individual entered the study so that the baselinehazard could depend on when MS first manifested itself clinically.The stratified Cox model may be extended even further to enable individualsto move between strata during follow-up (but only belonging to one strata ata time). The strata are thus time-dependent. Time-dependent strata wouldbe more appropriate in this study if the form of the baseline hazard does notdepend on when the disease was diagnosed but if it changed for everyone asmedical practices change. For instance, a person diagnosed with MS in 1996and a person diagnosed with MS in 2006 may have similar baseline hazardfunctions in 2006 if they both benefit from the same advances in medicine;when their disease began may matter less. Stratifying by calendar periodmay be appropriate.Two disadvantages of the stratified model are that no estimate is providedfor the impact of the stratification variable on the hazard and that there isa loss of power [18]. The effect of strata on the hazard may not be directlyestimable but strata-specific coefficients may offer insight indirectly intothe effects of the stratification variable. The loss of power is greater as thenumber of strata increases but does not become problematic unless there aremany strata relative to the sample size or there are only a few individuals insome strata. This may be attenuated by collapsing strata that are similarand have a small sample size. Moreover, the reduction in power is negligiblewhen many subjects are event-free at the end of follow-up, as is true in thisstudy where only 10.6% of individuals died.Partial LikelihoodTime-to-event data is the realization of several random processes: the eventprocess, the censoring process, and the covariate process [17]. A completelikelihood for the Cox model would incorporate each of these processes. Thisis typically not done.192.1. Intensity-based modellingThe goal of time-to-event analysis is to understand the effect of covariateson the event process. These effects are the parameters of interest. The co-variate and censoring processes are said to be non-informative if they do notcontribute additional information about the parameters of interest beyondwhat may be gleaned from the event process alone. Under the assumptionthat they are non-informative, the covariate and censoring processes couldof course provide information on the forces that impact the censoring andcovariate distributions, but this is rarely the objective. For the purposes ofmost studies, this extra information is not pertinent.A full likelihood may sometimes be factored into separate terms that con-tain only the parameters of interest and the nuisance parameters. If thiscan be done, the term that contains the parameters of interest is referred toas the partial likelihood [19]. In the context of time-to-event data, the termof the full likelihood that only incorporates the event process is the partiallikelihood. Standard maximum likelihood theory still applies to the partiallikelihood, and it has many of the useful properties of full likelihoods [20].Let ti be the observed event time of individual i and let R(t) be the riskset at time t, i.e., the individuals who have not experienced the event be-fore time t and are not censored. Assume that individuals are observedin continuous time and there are no tied event times. In practice, time isrecorded discretely, e.g., in days, so adjustments for tied event times aremade. These are ignored in this development. Let Ni(t) = 1 if individuali observably experienced the event by time t and 0 otherwise. Then, let∆Ni(t) = Ni(t+∆t)−Ni(t−) and let dNi(t) = lim∆t→0∆Ni(t). The partiallikelihood of the event process is constructed as follows [20].First, note that at ti, by (2.2), the probability that any individual experi-ences the event is∑j∈R(ti) λ(ti|xj(ti))∆t + o(∆t), if the individuals are in-dependent of one another, and the probability that individual i experiencesthe event is λ(ti|xi(ti))∆t + o(∆t). Then, the probability that individual iexperiences the event at ti given that an event has occurred at ti is:P (∆Ni(ti) = 1 | ∆N(ti) = 1,H(t)) = λ(ti|xi(ti))∆t+ o(∆t)∑j∈R(ti) λ(ti|xj(ti))∆t+ o(∆t).202.1. Intensity-based modellingReplacing the intensity function with the Cox model specification yields:=λ0(ti) exp(xi(ti)′β)∆t+ o(∆t)∑j∈R(ti) λ0(ti) exp(xj(ti)′β)∆t+ o(∆t).The baseline hazard function and ∆t may be factored out of the numeratorand the denominator, leaving=exp(xi(ti)′β) + o(1)∑j∈R(ti) exp(xj(ti)′β) + o(1).Thus, the baseline hazard function does not need to be estimated in orderto estimate the coefficients. Furthermore, o(1)→ 0 as ∆t→ 0, so:P (dNi(ti) = 1 | dN(ti) = 1,H(t)) = exp(xi(ti)′β)∑j∈R(ti) exp(xj(ti)′β). (2.5)Finally, let D be the set of individuals that are observed to experience theevent. If we assume that the event at time t depends only on H(t), thepartial likelihood of the Cox model is∏i∈DP (dNi(ti) = 1 | dN(ti) = 1,H(t)) =∏i∈Dexp(xi(ti)′β)∑j∈R(ti) exp(xj(ti)′β).Estimates for the Cox model are derived by maximizing this partial like-lihood. Typically, along with the maximum partial likelihood estimate,confidence intervals are created using the Wald statistic and p-values arereported from the Wald test [18].Hazard RatiosResults from Cox regression are rarely interpreted in terms of the coefficientsbut as a hazard ratio (HR). If the effect of some drug intervention on thehazard is of interest, then let X = 1 if an individual received treatment andX = 0 if they did not. Then, the hazard ratio of treatment is:λ(t | X = 1)λ(t | X = 0) =λ0(t) exp(1 ∗ β)λ0(t) exp(0 ∗ β) = exp(β),where β is the coefficient for treatment. When the covariate is continuous,the HR for X = x1 relative to X = x2 is exp(β(x1 − x2)). However, unlessotherwise specified, the HR of a continuous variable is reported as the HR212.1. Intensity-based modellingcorresponding to a change of one unit in X, namely exp(β). If the Coxmodel has many variables, the HR reported for any variable is interpretedas the effect of a one unit change in that variable on the hazard, holding theother variables constant.Proportional hazards assumptionThe form of the Cox model in (2.3) assumes that the impact of a variableon the hazard function is constant, i.e., the HR does not depend on time.This is the proportional hazards assumption. More complicated specifica-tions of the Cox model that do not assume proportional hazards are possiblethrough the inclusion of time-dependent coefficients. Time-dependent coef-ficients, denoted β(t) for the value of the coefficient at time t, allow the HRto change over time. Time-dependent coefficients may exist in many forms.For example, β(t) may be discontinuous, i.e., a step function over intervalsof t, or it may be continuous of the form β(t) = β∗f(t) where f is a specifiedfunction of time. Common choices of f are f(t) = log(t) or simply f(t) = t.The time-dependent component of the HRs in this thesis are reported inyears, e.g., for β(t) = β ∗ t, it is reported as the change in the HR in oneyear. The effect size in terms of days would be so small as to not be insightful.The HR of a time-dependent coefficient may be interpreted as a functionof time, i.e., the HR at time t is HR(t) = exp(β(t)). Time-dependent co-efficients complicate the interpretation of the hazard ratio but if a variablehas a time-dependent effect and it is ignored in the model specification, theresulting estimate of β is estimating a weighted average of β(t) over time.This may obfuscate the true impact of a covariate on the hazard.Time-dependent ExposureExposure in this study is time-dependent. To accommodate this, an indi-vidual’s follow-up could be divided into segments of time with each newsegment of follow-up time beginning when an individual’s exposure statuschanged. Each segment was represented by a separate row in the dataset.This is the “start/stop” data format [18].ResidualsDfbeta residuals are the changes in a coefficient estimate if a segment offollow-up time were to be removed from the dataset. For every segment of222.2. Non-parametric analysis of mortalitytime, dfbeta residuals are calculated for each covariate and are commonlyscaled by the standard error of that coefficient. As individuals may haveseveral rows in the dataset, residuals for each of their segments of time aresummed. Therneau and Grambsch state that the sum of scaled dfbeta resid-uals less than 1/3 are considered small [18]. If an individual’s sum of scaleddfbeta residuals is greater than 1/3, they are removed from the dataset andthe model is refit to determine how much the estimates change.Another important residual for the Cox model is the Schoenfeld residualwhich is calculated for each coefficient at each event time tk. Grambschand Therneau [21] showed that Schoenfeld residuals may be scaled by thevariance of the covariate vector at time tk to investigate the proportional haz-ards assumption. Grambsch and Therneau demonstrated if skj is the scaledSchoenfeld residual for βˆj at time tk then E(skj)+ βˆj is approximately equalto βj(tk), the jth time-varying coefficient at time tk. If the proportionalhazards assumption holds, then the mean of the scaled Schoenfeld residu-als will be approximately constant across time. If it is not, this suggeststhat a time-dependent coefficient should be used. Potentially time-varyingeffects of the covariates were investigated using the cox.zph function fromthe survival package. cox.zph models skj using a natural spline with fourd.f. This was done for each coefficient in the Cox model. If the spline isroughly constant over time, the coefficient is assumed to not vary over time.Deviance residuals compare the observed event time to the event time ex-pected under the model. Negative residuals indicate that the model expecteda small event time but a large one was observed or that the individual wascensored at the end of follow-up. Positive residuals indicate the opposite—the model expected long survival but the individual died quickly. Devianceresiduals can show if many individuals are poorly fit by the model and helpidentify which individuals the model does not explain well.2.2 Non-parametric analysis of mortalityIf T is the time of death, the survival function at time t is S(t) = P (T > t),the proportion of subjects still alive at t. When individuals may be censored,this cannot be estimated by (# alive at t / sample size) since this would beunbiased only if every censored individual had died.In the presence of censoring, the survival function at time t is commonly232.2. Non-parametric analysis of mortalityestimated with the Kaplan-Meier estimator. Let D be the set of timeswhere at least one person died during (0, t]. The Kaplan-Meier estimator isŜ(t) =∏u∈D(1− Number of individuals that died at uNumber of individuals at risk at u),where the denominator is the number of individuals in the risk set R(u) asit was defined in Section 2.1.2.The probability of survival given some categorical covariate X is denotedS(t|X = x) and may be estimated by calculating the Kaplan-Meier estimateusing only individuals who have X = x.Figure 2.1 shows the Kaplan-Meier estimate across the duration of the studyfor the cohort. At the end of the study, the estimated survival is 77.4%. Thedecrease in survival appears steady throughout the study.Figure 2.2 shows estimated survival given comorbidity at the index date.This figure is as one would expect: individuals with a lower comorbidityburden had better survival throughout the study. By the end of the study,the Kaplan-Meier estimate of survival was 82.2% among individuals withno comorbidity at baseline, 70.1% for individuals with a comorbidity of one,and 45.8% for individuals with a comorbidity of two. The clear separation ofthe curves indicates that even baseline comorbidity (which almost certainlybecomes less prognostic as the study goes on) is an important predictor.Figure 2.3 plots the estimated survival curves for individuals ever exposedversus never exposed. This shows clearly that individuals exposed to DMDsat some point during follow-up had a better prognosis but this is unadjustedfor any other covariates. The Cox models in the next section should providea better picture as to how DMDs may impact survival and will account forthe time-varying definition of exposure which was ignored here.Figure 2.4 shows survival curves stratified by socioeconomic status withoutpointwise confidence intervals. Individuals whose SES differs by only onequintile seem to have relatively similar survival but differences in survivalare clearer for individuals whose SES differs by two quintiles, and there isa marked separation between the curves for individuals in the lowest andhighest quintiles.242.2. Non-parametric analysis of mortalityFigure 2.5 shows Kaplan-Meier estimates of survival for individuals withno pre-index hospitalizations (n = 4,184) and individuals with at least onepre-index hospitalization (n = 4,940). There was a clear survival advantageamongst individuals with no pre-index hospitalizations.Figure 2.6 shows the Kaplan-Meier estimates of survival stratified into twogroups by age. The median age at index date is 44 years old and individualsare stratified based upon whether they are younger than 44 or not at theindex date. Of the six figures in this section, this figure shows the largestseparation between survival curves, indicating the prognostic power of age.Figure 2.1: Estimated survival for the cohort over the duration of the study.252.2. Non-parametric analysis of mortalityFigure 2.2: Estimated survival given comorbidity at the index date: comor-bidity of 0 = black, comorbidity of 1 = red, comorbidity of 2 = blue.Figure 2.3: Estimated survival based upon exposure at any point duringfollow-up: never exposed = red, exposed at some point during follow-up =blue.262.2. Non-parametric analysis of mortalityFigure 2.4: Estimated survival stratified by socioeconomic status quintile:1 = blue, 2 = red, 3 = black, 4 = green, 5 = purple.Figure 2.5: Estimated survival stratified by pre-index hospitalizations: zero= blue, one or more = red.272.3. The association between disease-modifying drugs and mortalityFigure 2.6: Estimated survival stratified by age at index date: younger than44 years old = blue, 44 or older = red.2.3 The association between disease-modifyingdrugs and mortality2.3.1 Limitations due to group sizesOne of the goals of this study is to investigate whether the impact of DMDson mortality changed for individuals with comorbidities. However, only 16individuals exposed to DMDs died with a baseline comorbidity of one andsix subjects exposed to DMDs died with a baseline comorbidity of two. Theremaining subjects that died had no comorbidities at baseline. Any compar-ison of the effect of DMDs on individuals with comorbidities would be basedon estimates calculated using a very small sample size. When comparing lowversus high exposure, the sizes of these groups are even smaller (e.g., fewerthan six individuals with a comorbidity of two died while exposed to DMDsat a high level). The accuracy of estimates are limited by the relativelyfew events experienced by these groups. However, in total, individuals inthis study experienced far more hospitalizations (25,792) than deaths (965),so it is more feasible to determine whether the association of DMDs andhospitalization was consistent across comorbidity levels. This goal is thusdeferred to Chapter 3.282.3. The association between disease-modifying drugs and mortality2.3.2 Model selectionThe association between DMD exposure and mortality in the cohort wasmodelled using Cox regression. The model was adjusted for socioeconomicstatus, sex, age at baseline, number of hospitalizations in the five years priorto the index date, and comorbidity at baseline.Three candidate models were considered: an unstratified model, a modelwith fixed stratification by index period, and a model with time-dependentstratification by calendar period. The stratified models had five strata: 1996- 1999, 2000 - 2004, 2005 - 2009, 2010 - 2014, 2015 - 2017. Time-dependentstratification was handled as a time-dependent variable in the dataset, i.e., anew segment of follow-up time if an individual changed strata. The outputof these models may be found under Section 1.1 of Appendix B.During the comparison of the three models, all were adjusted for the samevariables, and the variables were assumed to have a linear functional form.The time-dependent stratified model had the lowest AIC and was selectedfor the final model. After selection, the functional form of the continuouscovariates was investigated for the final model: there was no evidence that alinear form was unsuitable for age but hospitalizations in the five years beforeindex date appeared to be modelled better using a log(x+1) transformationas shown in Appendix B.1.2. Coefficients were assumed to be constant acrossstrata. The final form of the primary exposure model is:λk(t) = λk0(t) exp(β1Exposure + β2Sex +β3Age + β4SES2 + β5SES3 + (2.6)β6SES4 + β7SES5 + β8Comorbidity1 +β9Comorbidity2 + β10log(Hospitalizations + 1)),where k represents the time-dependent strata by calendar period. Schoen-feld residuals did not show any meaningful time-dependent effect for anyof the coefficients nor did scaled dfbeta residuals identify any influential in-dividuals. Residual plots and model assessment for the models using theprimary measures of exposure may be found in Appendix B.1.2.The model with the secondary definition of exposure was specified similarly.The functional forms of the continuous covariates were not changed norwas the stratification. Schoenfeld residuals were used to investigate time-dependent effects. Low DMD exposure exhibited a time-dependent effect292.3. The association between disease-modifying drugs and mortalityafter 4000 days post-index as shown in Figure B.5. Residual plots andmodel assessment for the models using the secondary measure of exposuremay be found in Appendix B.1.3. The final form of the model was:λk(t) = λk0(t) exp([β1 + β2 ∗ I(t ≥ 4000) ∗ (t− 4000)] Low exposure +β3High exposure + β4Sex + β5Age + β6SES2 +β7SES3 + β8SES4 + β9SES5 + β10Comorbidity1 +β11Comorbidity2 + β12log(Hospitalizations + 1)).2.3.3 Some comments on interpreting resultsTo interpret the models, we focus on results showing strong evidence of aneffect. Broadly speaking, this is where p < 0.05 and the confidence intervaldoes not include zero. Results where this is not the case are often reportedusing the shorthand that a covariate was “not associated” with the outcome.This should not be read as concluding that there is no effect. That state-ment would be too absolute. Instead, it merely suggests that that the datais consistent with no covariate effect. However, the data is also consistentwith effects contained in the confidence interval.In spite of this, interpretations are often worded as if such effects are zero,i.e., a HR of 1, even though the effect of a covariate will never be exactlyzero. All covariates have some effect; it just might be exceedingly small.If the evidence of an effect is weak or non-existent, or if the effect size issmall, then the covariate can be interpreted as having no meaningful clinicaleffect on the outcome. This is especially true in the next chapter: if theeffect of exposure does not meaningfully change across comorbidity levels,it is often concluded that “exposure does not appear to change for differentcomorbidity burdens”.2.3.4 Results of survival analysisThe estimated hazard ratios of the primary and secondary exposure mea-sures are shown in Table 2.1. Model output from R and a brief discussionof the other covariates from the models may be found in Appendix B.1.2(for the model using the primary measure of exposure) and Appendix B.1.3(for the model using the secondary measure of exposure). Exposure to anyDMD does not appear associated with mortality [HR = 1.02, CI = (0.80,302.3. The association between disease-modifying drugs and mortality1.29) p-value = 0.88]. Despite estimates further away from 1, this does notappear to change for individuals with low exposure before t = 4000 [HR =0.82, CI = (0.57, 1.16), p-value = 0.26] or high exposure [HR = 0.87, CI= (0.59, 1.31), p-value = 0.51]. There was a time-dependent effect of lowexposure. After 4000 days post-index, the fitted model indicated that theHR for low exposure increased 18% annually [HR = 1.18, CI = (1.07, 1.29),p-value < 0.001].Table 2.1: HRs of Exposure MeasuresCovariate HR CI p-valuePrimary exposure modelExposure 1.02 (0.80, 1.29) 0.88Secondary exposure modelLow exposure before 4000 days post-index 0.82 (0.57, 1.16) 0.26Time-dependent effect of low exposureafter 4000 days post-index (in years) 1.18 (1.07, 1.29) < 0.001High exposure 0.87 (0.59, 1.31) 0.512.3.5 Discussion of survival modelsThese results are inconsistent with previous work of the DRUMS projectreported by Ng et al. [1], where it was concluded that exposure to anyDMD was associated with lower mortality. However, there are importantdifferences in the design of the two studies:• Ng et al. used provincial databases from BC, Manitoba, Saskatchewan,and Nova Scotia.• The cohort used in this study was a subset of the BC cohort of Ng etal.’s study. That study included the following individuals that wereexcluded from this study: individuals that were not in the MSP regis-tration and Premium Billing (R&PB) files for at least five years priorto their index date; individuals that had missing socioeconomic statusdata; individuals that were identified by the algorithm as living withMS prior to the beginning of the study (assigned an index date of Jan-uary 1, 1996); individuals identified before their 18th birthday (whoentered the study on their 18th birthday).312.3. The association between disease-modifying drugs and mortalityThe cohort of the present study had 131,431 fewer person-years offollow-up. Removed individuals indexed on January 1, 1996 had longfollow-up (about 93,995 person-years) since, in theory, they could befollowed for nearly 22 years.• The Cox model in Ng et al. was stratified based upon index date.Refitting the primary model (2.6) with the cohort of this study butstratifying on index period instead of calendar period did not changethe estimated HR of exposure by more than 0.002.The results are probably not due to differences in the modelling. Ng et al.reported combined results from four provinces, so to compare the resultsof the primary model with Ng et al.’s BC cohort directly, model (2.6) withcalendar period stratification was re-estimated using the BC cohort fromthat paper. The number of pre-index hospitalizations was not included inthe model since this data is not available for all of the individuals in Nget al.’s study. In any case, pre-index hospitalization was not influential inour model. Excluding it and refitting the primary model changes the HR ofexposure estimate by less than 0.01 [HR = 1.01, CI = (0.79, 1.28), p-value= 0.89], so this variable is likely not the reason that the conclusions differed.The output for the model using Ng et al.’s cohort may be found in Section 2of Appendix B. The estimated HR using the larger cohort was lower [HR =0.84, CI = (0.75, 0.95), p-value = 0.004] so, unlike with the models above,there appeared to be an association between DMD and reduced mortality.A plausible explanation is that the individuals excluded from this studywere older and further along in their disease course—their first demyelinat-ing event was before the beginning of the study—compared to the includedindividuals which could make them likelier to die. These individuals mayalso have lower rates of exposure, so the unexposed group in the Ng et al.study would be less healthy than the unexposed group in this study, and asa result, DMD exposure seemed more beneficial.Another paper investigating the link between mortality and DMDs by King-well et al. [22] concluded that more than 6 months of cumulative use of beta-interferon was associated with a 32% reduction in the risk of death but morethan 6 months of cumulative use of glatiramer acetate was not. Kingwell etal. analyzed older individuals than this study (mean age was over 50 versus45.8) but individuals included in the Kingwell et al. cohort died at a younger322.4. Factors associated with exposureage (mean age at death was 61 versus 67.8); the distributions of covariatesbetween the two cohorts were similar other than age. Age may be one reasonfor the different results but there are two other differences that may be moreimportant. Kingwell et al. matched cases and controls by disability level atstudy entry which was not available for our study. This is a potentially im-portant covariate that could not be adjusted for here. Further, Kingwell etal. investigated the effect of only two DMDs, beta-interferon and glatirmiracetate, separately whereas this study considers any disease-modifying drug.Kingwell et al. also found that high exposure (defined as ≥ 3 years of beta-interferon use) was the driving force of this risk reduction. In the currentstudy, the Cox model estimate of high exposure had a lower estimated HRthan for low exposure but the CI was quite wide, so there is no evidencethat high exposure has the same effect in this study. This may be a resultof the relatively few individuals (32) that had high exposure to DMDs anddied; low p-values and small CIs from such a small group are unlikely evenif an association is present.2.4 Factors associated with exposureIt is not obvious what determines exposure. For example, Figure 2.7 showsthe average number of pre-index hospitalizations amongst exposed individ-uals over time up until 1000 days post-index, after which the graph changedvery little. At 270 days post-index, the average number of hospitalizationsamong the exposed group is 1.38 but at the end of follow-up that number is1.06. There seems to be a disproportionate number of individuals with highpre-index hospitalizations that were exposed early. A table comparing thecharacteristics of individuals that were not exposed during follow-up to in-dividuals that were exposed during follow-up would miss this. Determiningwhich covariates are associated with earlier exposure is complicated whenthere are time-dependent covariates such as comorbidity.A time-to-exposure Cox model may be useful to investigate this. To createthe final model, only a model with stratification by calendar period was con-sidered. The remaining steps (i.e., investigating functional form and propor-tional hazards, calculating deviance and dfbeta residuals) of model buildingare the same as for the survival models above. Appendix B.3 provides anoverview of the time-to-exposure model.332.4. Factors associated with exposureFigure 2.7: Average pre-index hospitalizations in exposed individuals overtimeLet λE(t) be the intensity of exposure. The final form of the Cox model isλEk (t) = λEk0(t) exp(β1Sex + [β2 + β3 ∗ log(t)](Age) +β3SES2 + β4SES3 + β5SES4 + β6SES5 +β7(Comorbidity1) + β8(Comorbidity2) +β9log(Hospitalizations + 1)),where k represents the time-dependent strata by calendar period. The in-teraction of age at the index date with time was modelled using the logfunction. This is a common choice when the event times are right-skewedand works well for this model too.2.4.1 Results of the time-to-exposure modelA full summary of the final model may be found in Appendix B.3.3. Morepre-index hospitalizations, modelled using a log(x+1) transformation, wereassociated with a lower hazard of being exposed [HR = 0.88, CI = (0.82,0.94), p-value < 0.001)], contrary to what seemed to be suggested in Figure2.1. Older individuals had a lower hazard of exposure. Older individualswere exposed less. The HR at the index date for age was 0.95 [HR = 0.95,CI = (0.95, 0.96), p-value < 0.001] and the hazard ratio of the interaction342.5. Exploring reverse causality using weighted cumulative exposurebetween age and time decreased slightly over each year [HR = 0.99, CI =(0.98, 0.99), p-value < 0.001]. A comorbidity of 1 was associated with a17.5% decrease in the intensity of exposure and a comorbidity of 2 was asso-ciated with nearly a 30% decrease. The other covariates were not associatedwith the intensity of exposure.2.5 Exploring reverse causality using weightedcumulative exposureDanieli et al. used weighted cumulative exposure (WCE) to investigate re-verse causality [23]. Weighted cumulative exposure attempts to summarizethe cumulative effect of an exposure that changes in status or dosage overtime. This is difficult to capture in a single categorical variable, and it isnot always clear how the effects of such exposures accumulate over time [24].The WCE at time t (where t is the time since the index date) is a weightedsum of past exposure history. Let u, which is specified beforehand, repre-sent the number of days before time t for which exposure is relevant to thehazard at time t. Then, each day between time t− u and time t is assigneda weight based on how exposure on that day alters the hazard at time trelative to an unexposed individual. If exposure on a given day is associatedwith an increase in the hazard, that day is assigned a positive weight. Dayswhen the exposure is expected to decrease the hazard are assigned negativeweights. Any exposure before time t−u is assumed to have no effect on thehazard at time t. Cox models that use WCE as a covariate may be fit in Rusing the WCE package [25].In the paper that introduced WCE modelling, exposure history for an in-dividual at time t is defined as occurring between time 0 and time t, andthe application section of the paper uses a cohort that was assumed to haveno exposure prior to the beginning of the study. That the WCE packageassumes exposure cannot occur before time 0 may be problematic in othersettings but is consistent with the cohort in this study. Thus, any intervalof the time window [t− u, t] that falls before the index date is treated as noexposure.Weights are assigned via the weight function w. Let w(j) be the weightcorresponding to the jth day before time t. To facilitate estimating of theweight function from the data, the weight function w is represented as a352.5. Exploring reverse causality using weighted cumulative exposureB-spline with m internal knots of the form: w(j) =∑mk=1 θkhk(j), wherehk is the kth B-spline basis function and θk is its coefficient. B-splinesare a computationally efficient method to estimate a spline function [26].A vignette included in the WCE package documentation indicates B-splinesthat have 3 internal knots or less are sufficient for all but the most complexWCE functions. Let xe(j) represent the exposure of an individual at timej, j = 1, . . . , t. The WCE at time t is calculated as:WCE(t) =t∑j=t−uw(t− j)xe(j) (2.7)A Cox model with only the WCE as a covariate takes the form:λ(t) = λ0(t) exp[β WCE(t)](2.8)The term in the exponent can be expressed as:β WCE(t) = βt∑j=t−uw(t− j)xe(j)= βt∑j=t−um∑k=1θkhk(t− j)xe(j) =m∑k=1βθkt∑j=t−uhk(t− j)xe(j)Sylvestre and Abrahamowicz [24] simplify this expression by letting let γk =βθk and Dk(t− u, t) =∑tj=t−u hk(t− j)xe(j) to re-express (2.8) asλ(t) = λ0(t) exp[ m∑k=1γkDk(t− u, t)]. (2.9)Once the Dk variables are calculated, the coefficients are estimated usingcoxph from the survival package [24], [25]. The γks are the coefficientsreturned by WCE function.It is not possible to unambiguously show reverse causality is present in anobservational study since the complete information on how treatment de-cisions are made is not available. At best, there may only be a suggestionthat it exists. If the weight function of a WCE suggests an effect of exposurewhen it is clinically implausible for one to exist (e.g., exposure to DMD afterone day is associated with a marked decrease in the hazard) then reversecausality may be present.362.5. Exploring reverse causality using weighted cumulative exposureFor example, Danieli et al. investigated the association between radiationfrom cardiac procedures and cancer incidence [23]. It was not known whethersuch radiation increased cancer incidence or whether individuals exhibitingsymptoms of cancer underwent more procedures. The weight function ofa WCE Cox model indicated that the hazard of cancer incidence was as-sociated with cardiac procedures done in the previous two years. Cancerdue to radiation of the type being investigated takes more than two yearsto manifest, and Danieli et al. concluded that reverse causality was likelypresent: doctors ordered more procedures on individuals that had cancer inresponse to the symptoms of the cancer.2.5.1 Results of weighted cumulative exposure Cox modelA WCE in a Cox model that was also adjusted for the covariates outlined inSection 1.1 was fit to the mortality data of the cohort using the WCE pack-age [25]. The weight function was restricted to go to zero at the end of a10-year time window. Splines with 0, 1, 2, or 3 internal knots placed evenlywere tested and the model using 0 internal knots in the weight function wasselected using BIC. Figure 2.8 shows the estimated 10-year weight functionfrom this model and Figure 2.9 shows all four weight functions that wereconsidered.The fitted weight function suggests that more recent exposure to DMDs isprotective against mortality but earlier exposure, more than 1200 days be-fore, seems to have the opposite effect. One would expect to see the weightfunction go towards zero but not become positive, i.e., after the benefits ofDMD exposure wear off, exposure in the past would be expected to have noeffect on the hazard of death. The weight function in Figure 2.8 suggeststhat the HR is lowest at the onset of exposure and appears to be greatest(HR exceeds one) about 3.5 years after exposure.372.5. Exploring reverse causality using weighted cumulative exposureFigure 2.8: Weighted cumulative exposure function over 10-year window forhazard of death over a 10-year window(a) Weight functions with 0, 1, 2, and 3internal knots on the same graph(b) A closer look at Figure 2.9aFigure 2.9: Weighted cumulative exposure function over 10-year window forhazard of death over a 10-year window different numbers of knotsThis shift in the weights from negative to positive around day 1200 implieseither short-term (≤ 1200 days) benefit but long-term (> 1200) harm as-sociated with the exposure or withdrawal effects from the exposure thatarise after day 1200 [27]. Harm associated with long-term use of DMDswas investigated by Luchetta et al., who conducted a systematic review ofextension studies using participants from randomized controlled trials [28].This review pointed out issues and a lack of transparency in many of thestudies that precluded comparison and argued that no conclusions aboutwhether the safety of DMDs changed or stayed the same over time could bedrawn. Previous work on the association between adverse events and DMDs382.5. Exploring reverse causality using weighted cumulative exposurein the DRUMS database [29] broadly agreed with results from clinical trials,when they were possible to compare. Clinical trials had a much shorter du-ration suggesting similar safety between short-term and long-term DMD use.If long-term DMD use is not associated with worse safety, then the shapeof the weight function may be due to withdrawal symptoms. A drug whereusers experience withdrawal symptoms after discontinuation would have aweight function that looked like the estimated function in Figure 2.8 [27].The negative weights before day 1200 would represent the benefits offered byexposure to the drug and the positive weights afterwards indicate the higherhazard due to harmful withdrawal symptoms. However, reviews on DMDdiscontinuation make no mention of withdrawal symptoms associated withDMDs [12], [13] and withdrawal symptoms emerging only after 1200 daysseems improbable. Instead, this shape may not be a result of withdrawalsymptoms but reverse causality. For reasons stated earlier, an individualmay stop treatment as their health worsens and they are at a higher risk ofdeath. This association appears just as withdrawal symptoms would exceptthat the outcome affects the exposure and not vice-versa.39Chapter 3An illness-death model toinvestigate the associationbetween disease-modifyingdrugs and hospitalizationsand mortalityIn this chapter, the limits of time-to-event models are briefly discussed beforeintroducing multi-state models and how they remedy this (Section 3.1). Sec-tions 3.2 and 3.3 introduce non-parametric and semi-parametric techniquesfor multi-state models. The former discusses methods of non-parametricestimation of features of interest in multi-state models (i.e., cumulative in-tensities, transition probabilities, and state occupancy probabilities) andthe latter expands the intensity-based techniques of the previous chapter tomulti-state models. The remainder of the chapter applies these techniquesto the DRUMS cohort in order to investigate how disease-modifying drugsare associated with hospitalization and mortality as well as if this variesacross comorbidity levels.3.1 Multi-state modelsConventional time-to-event models take a narrow view of diseases. Theyinvestigate the effect of covariates on a single outcome, often death, but thisis rarely enough to understand a disease. When individuals may experiencethe event multiple times, recurrent event models can be used. This hasthe same drawback of time-to-event models since only one type of event isconsidered. Competing risks models add a little more flexibility. In com-peting risks models, individuals may experience different types of events. Acommon example is when individuals may die from different causes. Thecompeting risk scenario assumes that individuals experience only one type403.1. Multi-state modelsof event.More commonly, physicians and researchers are interested in many eventssimultaneously during the course of a disease, some of which may occur morethan once. Time-to-event, recurrent events, or competing risks models areunable to accomodate these complexities together. For example, patientswith acute myeloid leukemia may experience remission, relapse, bone mar-row transplants, and death, and a multi-state model is able to investigatethese outcomes together [30].A multi-state model is specified through a finite number K of states and thetransitions between those states. The collection of states is called the statespace. A state represents the status of an individual at a particular time intheir follow-up, and a transition describes how an individual can move fromone state to another [31].States from which an individual can leave are transient states and states fromwhich individual cannot leave are absorbing states; an example of the latteris death. Transitions between states are dictated by intensity functions. Let{Nkl(t), t ≥ 0} be the number of transitions from state k to state l for allk ̸= l observed up to and including time t. Let ∆Nkl(t) be the number of kto l transitions observed in the interval [t, t+∆t). For the multi-state model,let N(t) = {Nkl(t), for all k ̸= l}. As in Section 2.1, let {Y (t), t ≥ 0} and{X(t), t ≥ 0} be the censoring and covariate processes respectively. Then,the intensity function of observed transitions from state k to state l, denotedλkl(t), is defined as follows:λkl(t) = lim∆t→0P (∆Nkl(t) = 1 | H(t))∆t, (3.1)where H(t) = {N(s), X(s), Y (s), 0 ≤ s ≤ t}. Let λkl(t) = 0 for anyone notat risk of a k to l transition at time t.The simplest multi-state model is a survival model: there are two states,alive and dead, and there is one transition from alive to dead. Multi-statemodels are almost always more complex than this. There are usually morethan two states such as in the case of patients with acute myeloid leukemia.Le-Rademacher et al. considered five states: initial treatment, complete re-mission, relapse, transplant, and death [30]. In practice, the number of statesin a multi-state model is limited by interpretability and sample size. Enoughpatients need to have experienced a transition for it to be adequately mod-413.1. Multi-state modelselled and too many states may make the model impractical. For example,in their study of diabetes-related disabilities, Lovblom et al. used a 9-statemodel but stated that a 27-state model would have been more realistic yetharder to interpret [32].In this study, multi-state models are applied to the DRUMS health ad-ministration data to determine the association between disease-modifyingdrugs and two outcomes: hospitalization and death. This analysis focuseson HRs for transitions between states, probabilities of state occupancy, andthe mean time spent in a state.3.1.1 Illness-Death ModelThe illness-death model is a common multi-state model with three states:healthy, illness, and death [17], [33]. Individuals usually begin in the healthystate and may transition from healthy to illness or from either to death.Some illness-death models allow individuals to recover, transitioning fromthe illness state back to the healthy state, and re-enter the illness state morethan once. This is the version of the illness-death model used in our analysis.The three states of the illness-death model in this study are out of hospital(O), hospitalized (H), and death (D). Every individual begins the study outof hospital and possible transitions are O → H, H → O, O → D, andH → D, as illustrated in Figure 3.1.Figure 3.1: Illness-Death model: O = out of hospital, H = hospitalized, D= death423.2. Non-parametric estimates3.2 Non-parametric estimatesNon-parametric estimates generally ignore covariates and describe the wholecohort. In some cases, one or two categorical covariates may be consideredto describe a subset of study subjects.3.2.1 Cumulative transition intensityThe cumulative intensity function for transitions from state k to state l isΛkl(t) =∫ t0 λkl(u)du and may be estimated by the Nelson-Aalen (NA) esti-mator.Let Yik(t) be an indicator that is 1 if individual i occupies state k and isnot censored at time t and 0 otherwise, and let Y·k(t) =∑ni=1 Yik(t) be thenumber of people at risk in state k at time t.Let {Nikl(t), t ≥ 0} be the counting process of the number of k → l transi-tions experienced by individual i up to and including time t. The subscripti is made explicit in this section whereas it was suppressed for Nkl(t) in Sec-tion 3.1. That section concerned a single, generic individual whereas thissection uses information from every individual and requires a distinction tobe made between the individuals. Let N·kl(t) =∑ni=1Nikl(t) be the numberof k to l transitions experienced by all the individuals up to t. The NAestimate of the cumulative transition intensity is [17]:Λ̂kl(t) =∫ t0dΛ̂kl(u) =∫ t0dN·kl(u)Y·k(u), (3.2)a step function that changes by dΛ̂kl(u) at every event time u. dΛ̂kl(u) isthe number of people who experienced the transition at u divided by thenumber of people at risk of the transition at u.3.2.2 Transition probabilitiesLet Z(t) be the state occupied at time t. Then, for t > s, Pkl(s, t) =P (Z(t) = l | Z(s) = k) is the probability of occupying state l at time t giventhat state k was occupied at time s. Let P (s, t) be the K × K matrix oftransition probabilities where element k, l is Pkl(s, t). The Aalen-Johansen(AJ) non-parametric estimate of the matrix of transition probabilities is aK × K matrix denoted P̂ (s, t) where entry k, l is an estimate of Pkl(s, t)433.2. Non-parametric estimatesdenoted P̂kl(s, t).Let D be the set of all transition times in (s, t] for every individual. LetM̂(u) be a K ×K matrix where u ∈ D. At event time u:• the non-diagonal entries of M̂(u) are dΛ̂kl(u), the non-parametric es-timate of the intensity of transitioning from k to l at time u.• the diagonal entry of row k of M̂(u) is 1 −∑l ̸=k dΛ̂kl(u), i.e., 1 - thesum of the non-diagonal entries in row k.Then, the AJ transition probability estimates are:P̂ (s, t) =∏u∈DM̂(u), (3.3)the product of M̂(u) across every event time after s and up to and includingtime t.AJ estimates were originally derived for Markov processes but are validtransition probability estimates for non-Markov processes if the censoringand multi-state processes are independent [17]. This is assumed to be thecase in the DRUMS study. Administrative censoring is independent of themulti-state process; there is little reason to assume that censoring due toemigration has an effect on the multi-state process; so few individuals werecensored because of daclizumab or ocrelizumab dispensation that it likelydoes not make a difference. The potential effect of censoring is discussedfurther in Section 3.7.3.2.3 State occupancy probabilitiesThe probability of occupying state l at time t is Ql(t) = P (Z(t) = l).Then for any time t > s, Ql(t) =∑Kk=1 Pkl(s, t)P (Z(s) = k). If Y··(s) =∑Kk=1 Y·k(s) is the total number of individuals at risk at time s then P (Z(s) =k), the state k occupancy probability at time s, is estimated by Y·k(s)/Y··(s).Pkl(s, t) may be estimated by (3.3), so the estimated probability of occupyingstate l at time t isQ̂l(t) =K∑k=1P̂kl(s, t)Y·k(s)Y··(s).443.3. Intensity-based multi-state model analysisIf every individual is in state k at time s, then Ql(t) = Pkl(s, t) andQ̂l(t) = P̂kl(s, t). For the remainder of this thesis, unless otherwise stated,s = 0 and Q̂l(t) = P̂Ol(0, t) since every individual begins in state O.State occupancy probabilities may be conditional on covariates. The prob-ability of occupying state l at time t given some categorical covariate X = 1is denoted Ql|X=1(t) = P (Z(t) = l | X = 1). This can be estimated by cal-culating Q̂l(t) using only individuals with X = 1, denoted Q̂l|X=1(t). Theprobability of occupying state l given an individual is exposed (or not ex-posed) is of interest. However, when exposure is time-dependent, this cannotbe done. Instead, state occupancy probabilities may only be conditioned onexposure at a certain time. Let E(s) be the time-dependent exposure attime s, where E(s) = 1 is exposed at time s and E(s) = 0 is unexposed attime s. One suggestion by Cook and Lawless [34] is to compare Q̂l|E(s)=1(t)to Q̂l|E(s)=0(t) for different values of s and t (where s < t) to reveal the re-lationship between the probability of occupying state l and time-dependentexposure.3.3 Intensity-based multi-state model analysis3.3.1 LikelihoodThe full likelihood of a multi-state model is rarely used to estimate param-eters, and similar to the simpler survival context, the partial likelihood ismaximized instead. A more detailed discussion of the likelihood for multi-state models may be found in [17] but the main points are summarizedbelow. The partial likelihood of a multi-state model appears complicatedinitially but, with a few assumptions, simplifies dramatically. The usualassumptions are as follows:• The exact event times are known.• Two or more events cannot occur in the same vanishingly small intervalof time (i.e., simultaneously).• There are no shared covariate effects between transitions. This as-sumption is sometimes eschewed when there is a biological reason orto reduce the number of parameters to be estimated [35].• The censoring process and the multi-state process are conditionallyindependent given the history.453.3. Intensity-based multi-state model analysisThe censoring process is not modelled, and yet, the conclusions of anintensity-based analysis are ideally about the entire population, whichincludes censored individuals [36]. If censoring is independent, thenthe censoring process does not depend at all on the multi-state pro-cess. This is a strict assumption. A more reasonable one is that themulti-state process and censoring process are conditionally indepen-dent given the history, i.e., whether an individual is uncensored at timet may depend only on the states that have been previously visited andcovariates that are being adjusted for in the model.• The individuals are independent of one another.The likelihood of the multi-state model is then the product of thelikelihood contribution from each of the n individuals.To begin developing the likelihood, consider a generic individual from thecohort and define three counting processes that describe how this individualmoves through the states. The subscript i is once again suppressed and asdefined in Section 3.1, {Nkl(t), t ≥ 0} is the number of transitions fromstate k to state l experienced by this individual by time t. Let Nk(t) be avector of length K − 1 with elements Nkl(t) where l ̸= k = 1, . . . ,K, and letN(t) = (N1(t), N2(t), . . . , NK(t)).As defined in Section 2.1, let X(t) be the covariate process and let Y (t) bethe censoring process. Furthermore, let Yk(t) be defined as Yik(t) was inSection 3.2 (that is, it equals 1 if the individual is uncensored in state K attime t and 0 otherwise) except the subscript i is now suppressed. So, if theindividual is censored from the study at t then Y (t) = 0 and Yk(t) = 0 fork = 1, . . . ,K.Assume that the entire length of follow-up [0, C] is partitioned into Rintervals, [ur−1, ur), r = 1, ..., R, where it is assumed that the intervalsare small enough that at most one event may occur in any interval. Let∆ur = ur − ur−1 be the length of the rth interval [ur−1, ur).Let ∆Nkl(ur) = Nkl(u−r ) − Nkl(u−r−1) be the number of times the indi-vidual transitions from k to l in the rth interval [ur−1, ur). Similar tothe definition of Nk(t) and N(t) above, let ∆Nk(ur) be a vector of lengthK − 1 with elements ∆Nkl(ur) where l ̸= k = 1, . . . ,K, and let ∆N(ur) =(∆N1(ur), . . . ,∆NK(ur)).463.3. Intensity-based multi-state model analysisLet H(ur) be the history of the individual at ur. This includes informationfrom the multi-state process from time 0 to time u−r and information fromthe covariate and censoring processes from time 0 to time ur. The proba-bility of a transition occurring during [ur−1, ur) depends on H(ur−1).Writing out the partial likelihood of an intensity-based multi-state modelcan be broken down into five steps:i. Factor the full likelihood into the partial likelihood.ii. Represent the likelihood contribution of an individual from a singleinterval by the transition probabilities, which enables the next step.iii. Express these probabilities in terms of intensity functions.iv. Convert the likelihood into continuous time because the subjects wereobserved in continuous time.v. Factor the partial likelihood into components for each transition. Asshown below, this makes estimation of the HRs easier.Cook and Lawless [17] show that the full likelihood of the multi-state, cen-soring, and covariate processes can be factored into three terms such that theprobability of the realization of the multi-state process—the partial likeli-hood as stated in (3.4)—is separated from that of the covariate and censoringprocesses. This factorization relies on the assumption that the censoring andcovariate processes are non-informative. If the censoring and/or covariateprocesses are informative, they must be jointly modelled with the multi-stateprocess to produce unbiased estimates, so the full likelihood cannot then befactored into separate terms. The partial likelihood across all R intervals isR∏r=1P (∆N(ur) | H(ur−1)). (3.4)Second, we evaluate the likelihood contribution for a single interval usingprobabilities of transitions between state k and state l. From (3.4), thecontribution for a single interval is clearlyP (∆N(ur) | H(ur−1)) . (3.5)∆N(ur) captures every possible transition in the rth interval [ur−1, ur) viathe vectors ∆Nk(ur), k = 1, . . . ,K. The goal in this step is to express (3.5)473.3. Intensity-based multi-state model analysisin terms of probabilities related to the elements of the vectors ∆Nkl(ur),k ̸= l = 1, . . . ,K.This depends crucially on assuming that there may only be one transitionin any interval. We only need to consider whether the individual made atransition from the state occupied at the beginning of the interval to anotherstate and do not have to incorporate the possibilities that arise when severaltransitions are possible in the interval (e.g., the individual may transitionfrom k to l then from l to any other state then from that other state toanother and so on).As a result, only one element of ∆N(ur)—the element corresponding to thestate occupied at the beginning of the interval—contributes to the likeli-hood in an interval. This allows (3.5) to be expressed as the product ofprobabilities corresponding to each element of ∆N(ur):K∏k=1P (∆Nk(ur) | H(ur−1))Yk(ur−1). (3.6)Let ∆Nk·(ur) =∑l ̸=k∆Nkl(ur) be the number of transitions out of state kexperienced by this individual during the rth interval, which may be 0 or 1.Then (3.6) may be expressed asK∏k=1K∏l ̸=k=1P (∆Nkl(ur) | H(ur−1))Yk(ur−1) (3.7)=K∏k=1K∏l ̸=k=1P (∆Nkl(ur) = 1 | H(ur−1))∆Nkl(ur)×P (∆Nk·(ur) = 0 | H(ur−1))Yk(ur−1)(1−∆Nk·(ur)).(3.8)When a k to l transition occurs in the interval, ∆Nkl(ur) = 1 and 1 −∆Nk·(ur) = 0 so the term outside the parentheses reduces to one and theterm inside the parentheses is the contribution to the likelihood. When notransition out of state k occurs, ∆Nkl(ur) = 0 for l ̸= k and 1−∆Nk·(ur) =1, so the term inside the parenthesis is one and the term outside the paren-thesis is the contribution. Neither term contributes to the likelihood if theindividual does not occupy state k at the beginning of the rth interval. Forthe term outside the parentheses, this is easy to see since Yk(ur−1) = 0.For the term inside the parentheses, note that ∆Nkl(ur) must be zero if483.3. Intensity-based multi-state model analysisYk(ur−1) = 0.Third, the partial likelihood of (3.8) may be described in terms of intensityfunctions. By (3.1), the probability in the first term of (3.8) can be expressedin terms of the k to l transition intensity and the probability of no transitionoccurring in the second term can be expressed as 1−λk·(ur | H(ur−1))∆ur+o(∆ur) where λk·(ur | H(ur−1) =∑Kl ̸=k=1 λkl(ur | H(ur−1)). So, (3.8) maybe rewritten in terms of intensity functions asR∏r=1K∏k=1K∏l ̸=k=1[λkl(ur | H(ur−1))∆ur + o(∆ur)]∆Nkl(ur)× [1− λk·(ur | H(ur−1))∆ur + o(∆ur)]Yk(ur−1)(1−∆Nk·(ur)) .(3.9)The expression is simplified slightly by noting that by a Taylor series expan-sion:1− λk·(ur | H(ur−1))∆ur + o(∆ur) ≈ exp (−λk·(ur | H(ur−1))∆ur)= exp K∑l ̸=k=1−λkl(ur | H(ur−1))∆ur = K∏l ̸=k=1exp (−λkl(ur | H(ur−1))∆ur) .Using this, and then dividing (3.9) by∏Rr=1∏Kk=1∏Kl ̸=k=1∆u∆Nkl(ur)r givesthe following expression for the partial likelihood:∝R∏r=1K∏k=1K∏l ̸=k=1[λkl(ur | H(ur−1)) + o(1)]∆Nkl(ur)× exp(−λkl(ur | H(ur−1))∆ur)1−∆Nk·(ur).(3.10)Yk(ur−1) is no longer necessary in the exponent on the second term as theintensity is defined as being 0 when Yk(ur−1) = 0, so the second term reducesto 1 when the individual does not occupy state k.Fourth, let R → ∞ since the individuals are observed in continuous time,not discrete intervals. As the number of intervals in [0, C] increases, as-sume that they all become smaller, and as that number approaches infinity,maxr(∆ur)→ 0.493.3. Intensity-based multi-state model analysisThe likelihood can now be written in terms of observed transition timesinstead of intervals where transitions occurred. Let Dkl be the set of tran-sition times from state k to state l experienced by the individual. Whent /∈ Dkl, ∆Nkl(ur) = 0 and the first term in (3.10) is one. So, the likelihoodcontribution when t /∈ Dkl islimR→∞R∏r=1K∏k=1K∏l ̸=k=1exp(−λkl(ur | H(ur−1))∆ur)= limR→∞K∏k=1K∏l ̸=k=1exp(−R∑r=1λkl(ur | H(ur−1))∆ur). (3.11)Noting that as R → ∞ with max(∆ur) → 0,∑Rr=1 g(ur)∆ur →∫ C0 g(u)dufor any integratable function g, (3.11) bcomesK∏k=1K∏l ̸=k=1exp(−∫ C0λkl(s | H(s))ds). (3.12)On the other hand, when t ∈ Dkl—after the limit R→∞ with max(∆ur)→0 is taken in (3.10)—the contribution to the likelihood is λkl(t | H(t)) since∆Nk·(ur) = 1. Thus, the contribution to the likelihood of the individualwhen treating time as continuous is:=K∏k=1K∏l ̸=k=1∏t∈Dklλkl(t | H(t)) exp(−∫ C0λkl(s | H(s))ds). (3.13)Fifth, the likelihood may be separated into distinct components that repre-sent each type of transition. Let Lkl be the contribution for transitions fromk to l as described in (3.13). The contribution to the partial likelihood bythe individual may be written as:K∏k=1K∏l ̸=k=1Lkl . (3.14)The partial likelihood for all individuals is obtained by taking the productof the contributions to the partial likelihood of each individual.503.4. Non-parametric analysisIf there are no shared effects between the transitions, the partial likelihoodcorresponding to each type of transition may be maximized separately. Tran-sitions between states may be modelled with separate time-to-event modelswhere the outcome of interest is the time-to-transition. For example, themulti-state model in Figure 3.1 may be modelled with four Cox models in-vestigating the following events: the time-to-hospitalization of individualsalive and out of hospital (OH transition); the time-to-discharge of individu-als currently hospitalized (HO transition); the time-to-death of individualsalive and out of the hospital (OD transition); the time-to-death of indi-viduals alive and hospitalized (HD transition). The likelihood L of theillness-death model in Figure 3.1 is L = LOHLODLHOLHD and, under theassumption of no shared parameters, the four likelihood components can bemaximized separately.3.4 Non-parametric analysisFigure 3.2 shows the NA estimates of the cumulative intensity of transi-tioning to state H and state D from state O. The cumulative intensity ofO → H transitions was much greater than that of O → D transitions. Bythe end of follow-up, the former was 44 times larger than the latter. Bothplots show the cumulative hazard increasing at a consistent rate.Figure 3.3 shows the NA estimates of cumulative intensity of discharge anddeath from stateH. These are on a different time scale (days since admissionto the hospital) than Figure 3.2 (years since index date). The NA estimatesfor transitions out of state H show that the intensity rose quickly early onand then increased more slowly. The H → O cumulative intensity is onlyabout 5.3 times larger than the H → D cumulative intensity at 2000 days.After this, the intensities did not change.513.4. Non-parametric analysis(a) O → H (b) O → DFigure 3.2: Cumulative intensity of transitions out of state O(a) H → O (b) H → DFigure 3.3: Cumulative intensity of transitions out of state HThe probability of being hospitalized at time t is shown in Figure 3.4. Hospi-talization was highest within 500 days of the index date (perhaps in responseto disease onset) but afterwards remained fairly consistent around 0.42%.The confidence intervals become wider after 6000 days since relatively fewerpeople (1,895) were followed for that long.The AJ-based estimate of the probability of dying by time t is shown inFigure 3.5. When time was divided into 1000 day intervals—e.g., [0, 1000),523.4. Non-parametric analysis[1000, 2000), and so on—the probability of dying was greater in each sub-sequent interval. The probability of dying was lowest (0.0240) in the first1000 days and highest (0.0324) in the final 1000 days. The probability ofdying by time t increased by an average of 0.0280 every 1000 days for theduration of follow-up.Figure 3.4: The probability of being hospitalized at time t. To reduce noisein the AJ estimate and CIs, the plot shows a moving average with a windowof 250 days.Figure 3.5: The probability of dying at time t. To reduce noise in the AJestimate and CIs, the plot shows a moving average with a window of 250days.533.4. Non-parametric analysisA non-parametric analysis for more than one or two covariates would takeup too much space and require too many plots. Instead, we focus only onthe probability of hospitalization at time t stratified by baseline comorbidityas shown in Figure 3.6. Confidence intervals made the figures too busy andhard to read so they were not included.There is a noticeable difference in the hospitalization levels for individualswith different baseline comorbidities. Individuals with no comorbidities atthe index date generally had the lowest hospitalization levels during follow-up. Individuals with a comorbidity of one had a higher probability of beinghospitalized during the study, and both groups had lower chances of beinghospitalized than individuals with a baseline comorbidity of two. The curvescrossed occasionally (e.g., the red and blue curves cross at about 3000 dayspost-index) but not often. The most noticeable feature of Figure 3.6 is thatindividuals with a comorbidity of two had exceptionally high hospitalizationlevels in the first 1000 days of follow-up relative to the other two groups. Thisappears to the driving force of the elevated probability of hospitalizationnear the start of follow-up in Figure 3.4. This feature abated but there werefairly sharp increases and decreases at different points throughout the study.Figure 3.6: The probability of hospitalization by time t stratified by baselinecomorbidity: comorbidity of 0 = black, comorbidity of 1 = red, comorbidityof 2 = blue. To reduce noise in the AJ estimates, the plot shows a movingaverage with a window of 250 days.543.5. Cox models for the illness-death model3.5 Cox models for the illness-death modelThe transitions between the states are represented through separate Coxmodels; the model selection procedure for the Cox model in the survivalsection was repeated here for each transition. First, three candidate modelswere considered: an unstratified model, a model with fixed stratification byindex period, and a model with time-dependent stratification by calendarperiod. The stratified models had five strata: 1996 - 1999, 2000 - 2004, 2005- 2009, 2010 - 2014, 2015 - 2017. Time-dependent stratification was handledas a time-dependent variable in the dataset, i.e., a new segment of follow-uptime if an individual changed strata.During comparison, all models were adjusted with the same covariates, andthe quantitative covariates were assumed to have a linear functional form.The time-dependent stratified model had the lowest AIC and was selectedas the initial working model in all four transitions. After selection, the func-tional form of the continuous covariates was investigated. Hospitalizations inthe five years before index date was more suitably modelled using a log(x+1)transformation in every transition. Age was modelled using a spline in theO → H and H → O transitions. Next, time-dependent coefficients were ex-plored using scaled Schoenfeld residuals for all covariates. If the Schoenfeldresiduals indicated a time-dependent effect, an interaction between time andthe covariate was included in the final model. Time-varying covariate effectsare specified in days when applicable. Deviance residuals and scaled dfbetaresiduals were calculated for the final model to check individuals whose eventtimes were poorly explained by the model and to determine any influentialindividuals respectively. Since individuals may have several events in theO → H and H → O models, each deviance residual corresponds to onesojourn.The model selection procedure was simplified for the secondary measure ofexposure. The functional forms of the covariates and the stratification re-mained the same as in the final model for the primary measure of exposure.The proportional hazards assumption was checked once more for all covari-ates and time-dependent effects were added when necessary.The main text has descriptions of the models (i.e., the variables included,the functional form of the continuous variables, any interactions with time)for the primary and secondary measures of exposure in Section 3.5, the HRsof exposure are reported (including interaction terms) in Section 3.6, and a553.5. Cox models for the illness-death modeldiscussion of the results is in Section 3.7. Supporting information for themodels may be found in Appendix C.Appendix C has separate sections for the OH, OD, HO, and HD Coxmodels, each of which are structured similarly. For example, Section 1 ofAppendix C titled ‘Transition OH’ contains detailed information on theCox models representing the OH transition. Section 1.1 of Appendix C hasdetails on the OH model using the primary measure of exposure. Tablesof the R output summarizing the candidate models are provided; AIC andlog-partial likelihood are summarized above every model in the appendix.The functional forms of the continuous covariates and the proportional haz-ards assumption are then checked. The scaled Schoenfeld residual plots forthe coefficients of exposure as well as for any covariate with an effect thatappeared to vary over time are included. Plots that did not suggest a time-varying effect are not shown. A summary of the final model is providedin ‘Model results’ before deviance and scaled dfbeta residuals are checkedin ‘Model diagnostics’. Finally, a brief discussion of covariates (i.e., ageat index date, sex, SES, pre-index hospitalizations, and comorbidity) canbe found under ‘Model interpretation’. Section 1.2 of Appendix C has theSchoenfeld residual plots for the secondary measure of exposure that justifythe inclusion of time-dependent coefficients, followed by model results andinterpretation.As discussed in Section 2.2.1, the number of individuals who died is insuffi-cient to investigate whether the effect of exposure on mortality varied overcomorbidity levels. Hence, interactions between comorbidity and exposurewere not considered for the O → D and H → D models. Interactions be-tween exposure and comorbidity were considered for the O → H andH → Omodels since many individuals had hospitalizations.3.5.1 O → H transitionPrimary measure of exposureThere was a violation of the proportional hazards assumption in three of thecovariates (as shown in Figures C.3, C.4, C.5): exposure to any DMD, num-ber of pre-index hospitalizations, and comorbidity. The spline fits to thescaled Schoenfeld residuals suggested exposure and pre-index hospitaliza-tions had relatively linear changes in effect across time, so time-dependentcoefficients for exposure and hospitalizations were added to the models563.5. Cox models for the illness-death modelthrough linear interactions with time.The change of both comorbidity variables over follow-up was not as straight-forward. Comorbidity changed during the first 2000 days and then plateauedwhich is more difficult to model. To accommodate this, comorbidity was in-corporated as a second stratifying variable. The final OH Cox model has15 strata, and covariate–exposure interaction terms were thus shared acrossstrata that had the same comorbidity but different calendar period.Age at the index date was modelled by a natural spline with three degreesof freedom, which is represented by three basis functions [26], here denotedh1(·), h2(·), and h3(·). Then, the Cox model representing the OH transitionis specified as:λOHk (t) = λOHk0 (t) exp([β1 + β2 ∗ t]Exposure + β3Sex + β4h1(Age) +β5h2(Age) + β6h3(Age) + β7SES2 + β8SES3 +β9SES4 + β10SES5 + [β11 + β12 ∗ t]log(Hospitalizations + 1) +β13Exposure ∗ Comorbidity1 + β14Exposure ∗ Comorbidity2),where k represents the time-dependent strata by comorbidity and calendarperiod, and β13 and β14 are shared across certain strata.Since individuals could make the O → H transition more than once, confi-dence intervals and p-values from Wald tests were calculated with a robuststandard error to account for within-individual dependence [18].Secondary measure of exposureThe proportional hazards assumption appeared to be violated for low andhigh exposure to any DMDs, the spline for age at the index date (the co-efficient for the first basis function appears to increase and then plateauafter t = 2000), sex, and pre-index hospitalizations. The coefficients forthese variables were incorporated in time-dependent forms suggested by thespline fits to the Schoenfeld residuals. The proportional hazards assumptionfor SES5 also appeared violated but the changes were relatively small andnot easy to handle using an interaction with time. This time-dependent ef-fect is noted but not modelled explicitly. The Cox model using the secondary573.5. Cox models for the illness-death modelmeasure of exposure representing the OH transition is specified as:λOHk (t) = λOHk0 (t) exp([β1 + β2 ∗ I(t ≤ 4000) ∗ t]Low exposure +[β3 + β4 ∗ I(t > 4000) ∗ (t− 4000)]High exposure +[β5 + β6 ∗ I(t > 4000) ∗ (t− 4000)]Sex +I(t ≤ 2000) ∗ [ β7h1(Age) + β8h2(Age) + β9h3(Age)] +I(t > 2000) ∗ [β10h1(Age) + β11 ∗ h2(Age) + β12h3(Age)] +β13SES2 + β14SES3 + β15SES4 +β16SES5 + [β17 + β18 ∗ t]log(Hospitalizations + 1) +β19Low exposure ∗ Comorbidity1 + β20High exposure ∗ Comorbidity2 +β21Low exposure ∗ Comorbidity1 + β22High exposure ∗ Comorbidity2).The models are specified using t in days, e.g., β2 is zero after 4000 days post-index, but the results reported in Section 3.6 for the O → H and H → Omodels are in years so that the effect is easier to interpret.3.5.2 O → D transitionPrimary measure of exposureExposure appeared to have a time-dependent effect, beginning to increaseafter 4000 days post-index (see Figure C.14). Individuals were at risk ofdeath in the survival model of the previous chapter for the entire durationof follow-up but are only at risk of the O → D transition when they arenot hospitalized. This is why the results from this model and the model ofChapter 2 differ. The Cox model for the O → D transition isλODk (t) = λODk0 (t) exp([β1 + β2 ∗ I(t > 4000) ∗ (t− 4000)]Exposure + β3Sex +β4Age + β5SES2 + β6SES3 + β7SES4 + β8SES5 +β9Comorbidity1 + β10Comorbidity2 + β11log(Hospitalizations + 1)),where k represents the time-dependent strata by calendar period.Secondary measure of exposureThe time-dependent effect of the primary measure of exposure appears tobe driven by the effect of low exposure to any DMDs. In addition, when the583.5. Cox models for the illness-death modelsecondary measure of exposure was used, pre-index hospitalizations had aneffect that decreased over time. The Cox model for the O → D transitionusing the secondary measure of exposure isλODk (t) = λODk0 (t) exp([β1 + β2 ∗ I(t > 4000) ∗ (t− 4000)]Low exposure +β3High exposure + β4Sex + β5Age + β6SES2 +β7SES3 + β8SES4 + β9SES5 + β10Comorbidity1 +β11Comorbidity2 + [β12 + β13 ∗ t]log(Hospitalizations + 1)),where k represents the time-dependent strata by calendar period.3.5.3 H → O transitionModel with primary measure of exposureThe transitions out of H are modelled on a semi-Markov time scale. In aMarkov Cox model, the baseline hazard function is a specified function of thetime since entering the study (the index date for this study), but the baselinehazard functions in the HO and HD Cox models are a function of the timesince entering the H state. Some covariates are also changed. Exposure,sex, SES, and pre-index hospitalizations are the same but comorbidity is nolonger time-varying. It is fixed at the value at hospital admission. Similarly,age at the date of admission was used instead of age at the index date sinceit is more useful and consistent with the time scale. Robust SEs were alsoused for this transition.Comparison of the three candidate models and examination of functionalforms for age at admission and the number of pre-index hospitalizationsproceeded as for the OH and OD transitions.More than half (13,825) of the 25,792 hospital visits where the patient wasdischarged alive were shorter than one day. These same day discharges maybe fundamentally different than visits longer than one day, e.g., same dayprocedures versus inpatient care. They require distinct consideration in themulti-state model since the effect of the covariates on same-day dischargesand longer hospital visits may be quite different. In addition, some hospi-talizations are unusually long. The longest hospitalization is 2,223 days and259 hospitalizations are longer than 90 days. Hospitalizations such as thesemay be more akin to long-term care, which falls outside of this study’s in-terest. After the stratification and functional form of continuous covariates593.5. Cox models for the illness-death modelhas been decided, we consider separate coefficients for sojourns of < 1 day,between 1 and 90 days, and > 90 days to estimate the separate effects.As a result, investigating violations of the proportional hazards assumptionproceeded differently for the OH models. Time-varying effects were inves-tigated for the second set of coefficients, reflecting covariate effects between1 and 90 days, only. Time-varying effects were not considered for the othertwo sets of coefficients. One day is too short for a time-dependent effect toexist. Time-varying effects longer than 90 days are unlikely to be modelledaccurately and are not of direct relevance. Detailed reasoning is given inSection 3.1 of Appendix C.If t∗ is the time since the most recent admission to the hospital, i.e., thecurrent sojourn/visit time, then the HO Cox model is of the following form:λHOk (t∗) = λHOk0 (t∗) exp(x(t∗)′β(t∗)),where k represents the time-dependent strata by calendar period.Exposure, sex, SES, comorbidity at admission, pre-index hospitalizations,and the interaction between exposure and comorbidity have three sets ofcoefficients. Age at hospital admission is modelled with a 3 d.f. spline. Thescaled Schoenfeld residual plots that indicate a violation of the proportionalhazards assumption may be found in Section 3.1 of Appendix C beginningwith Figure C.22. The coefficient of the third basis function in the splinefor age appeared to have an effect that increased over time when 1 ≤ t∗ ≤90. The effect, indicated in the Schoenfeld residual plot in Figure C.27,appears linear but was modelled with an indicator function to avoid evenmore complexity. The basis functions for the age spline thus have four setsof coefficients (meaning there are four different splines) corresponding tothe following intervals of sojourn time: t∗ < 1, 1 ≤ t∗ ≤ 10, 10 < t∗ ≤ 90,t∗ > 90. When 1 ≤ t∗ ≤ 90, SES5 and pre-index hospitalizations both hadan interaction with log time. There are 41 coefficients in the final model so,for brevity, the form of the final model is not stated here. The final modelfit may be found in Table C.14.Model with secondary measure of exposureThe scaled Schoenfeld residual plots changed very little between the modelwith the primary or secondary measure of exposure. Plots for low exposure,603.5. Cox models for the illness-death modelhigh exposure, and admission age may be found in Section 3.2 of AppendixC. The effect of neither low nor high exposure changed over time but thecoefficient of the third basis function in the spline for age changed over time.The same solution as in the model with the primary measure of exposure—four different splines—was also implemented for this model.3.5.4 H → D transitionModel with primary measure of exposureThe covariates and time scale are the same in the HO and HD Cox mod-els. Separate effects for different sojourn times are not considered in theHD model. The distinction between short-, medium-, and long-term hos-pitalizations does not apply to deaths. For example, a patient’s death onthe first day of their hospitalization does not imply that hospitalization wasshort-term, medium-term or long-term.None of the Schoenfeld residual plots indicated a violation of the propor-tional hazards assumption but the plots for exposure (Figure C.34) is pro-vided nevertheless.The final HD Cox model is specified asλHDk (t∗) = λHDk0 (t∗) exp(β1Exposure + β2Sex + β3(Age at admission) +β4SES2 + β5SES3 + β6SES4 + β7SES5 +β8(Comorbidity1 at admission) + β9(Comorbidity2 at admission) +β10log(Hospitalizations + 1)),where k represents the time-dependent strata by calendar period.Model with secondary measure of exposureThe exposure covariate is the only difference between the model with theprimary measure of exposure and the model with the secondary measure ofexposure. The other covariates remain the same and none exhibited a time-dependent effect in the scaled Schoenfeld residuals. The scaled Schoenfeldresidual plots for low and high exposure are shown in Figures C.36 and C.37.613.6. Results of multi-state intensity-based analysis3.6 Results of multi-state intensity-based analysisO → H transitionThere was evidence that exposure to any DMD was associated with thechances of hospitalization and that this effect varied over time. No indi-vidual was exposed at the index date—the earliest exposure was 96 dayspost-index date—so the initial estimate was extrapolated using the linearinteraction with time. It should be interpreted alongside the time-dependenteffect of exposure.At the index date, exposure for individuals with no comorbidity was associ-ated with a 25% decrease in the hazard of hospitalization [HR = 0.75, CI =(0.66, 0.85), p-value < 0.001]. There was no evidence that the effect of expo-sure to any DMD changed amongst individuals with different comorbidities.For simplicity, we proceed as if there is no change at all across comorbidityunless there was evidence to support the presence of interactions, as dis-cussed in Section 2.3.3. The HR of exposure changed over time and wasassociated with a 2.83% annual increase in the hazard [HR = 1.03, CI =(1.02, 1.04), p-value < 0.001]. At one year post-index, the HR of exposurewas 0.77, and at five years post-index, it was 0.86. The HR of exposure is 1at 10.4 years after the index date. After this point, exposure to DMDs wasassociated with a higher hazard of hospitalization. By 15 years post-index,exposure is associated with a 14% increase in the hazard of hospitalization.Both the protective effect of exposure to DMDs as well as the increase in theHR of hospitalization over time associated with exposure are largely drivenby the effect of high exposure to DMDs. For the first 4000 days, the HR ofhospitalization for high exposure among individuals with no comorbidity was0.79 [HR = 0.79, CI = (0.69, 0.90), p-value < 0.001]. After about 4000 dayspost-index, this HR increased by a factor of 1.06 every year [HR = 1.06, CI= (1.03, 1.09), p-value < 0.001]. In addition, there was some evidence that,under high exposure, there is a greater reduction in the hazard among thosewith a comorbidity of two compared to individuals with no comorbidities[HR = 0.74, CI = (0.57, 0.98), p-value = 0.03]. There was insufficientevidence to suggest that low exposure was associated with hospitalizations[HR = 0.91, CI = (0.81, 1.03), p-value = 0.14], that this association changedacross comorbidity levels, or, despite indications from Figure C.7, that therewas a time-dependent effect over the first 4000 days [HR = 1.01, CI = (1.00,1.03), p-value = 0.14].623.6. Results of multi-state intensity-based analysisO → D transitionAs indicated by Figure C.14, proportional hazards did not hold for expo-sure; around 4000 days post-index, the HR of exposure appeared to begin toincrease in a roughly linear manner. There was no evidence that exposureto any DMD had an association with mortality in the first 4000 days post-index [HR = 0.71, CI = (0.43, 1.17), p-value = 0.18]. After 4000 days, theHR increased by about 5% every subsequent year [HR = 1.05, CI = (1.01,1.09), p-value = 0.02].Low exposure to any DMD shared this pattern of association: there was noevidence that there is an association in the first 4000 days post-index [HR= 0.70, CI = (0.42, 1.16), p-value = 0.17] and then an annual increase inthe hazard ratio thereafter [HR = 1.25, CI = (1.11, 1.40), p-value < 0.001].High exposure to any DMD was not associated with mortality outside of thehospital [HR = 0.99, CI = (0.60, 1.63), p-value = 0.96].An interaction between comorbidity and exposure was not included in thismodel but comorbidity was prognostic of death outside of the hospital. Forthe primary measure of exposure, a comorbidity of one was associated witha HR of 1.88 [CI = (1.52, 2.33), p-value < 0.001] and a comorbidity of twowas associated with a HR of 3.10 [CI = (2.52, 3.81), p-value < 0.001] relativeto no comorbidities. The results for the model with the secondary measureof exposure are similar and can be found in Table C.10. This is a markedincrease over the effect of comorbidity shown for the time-to-death modelof Chapter 2 but that model used baseline comorbidity whereas this Coxmodel uses time-dependent comorbidity.H → O transitionThere was weak evidence that the effect of exposure changed across comor-bidity levels between the first and 90th day of hospitalization. Comparedto exposed individuals with no comorbidities, exposed individuals with acomorbidity of one had a nearly 16% greater HR of discharge [HR = 1.16,CI = (1.01, 1.33), p-value = 0.04]. This was also apparent for high exposure[HR = 1.24, CI = (1.02, 1.52), p-value = 0.035] but not low exposure. Noeffect was evident in the other sets of coefficients nor was there any evidencefor an effect for exposed individuals with a comorbidity of two.Exposure to DMDs was associated with a lower hazard of discharge in the633.7. Discussion of intensity-based analysisfirst day of hospitalizations although the effect was not large [HR = 0.94,CI = (0.89, 0.99), p-value = 0.02]. This also seemed to be the case for lowexposure to DMDs [HR = 0.91, CI = (0.86, 0.97), p-value < 0.001] but nothigh exposure [HR = 0.99, CI = (0.92, 1.08), p-value = 0.89].H → D transitionThere was no evidence that exposure to any DMD had an association withdeath during hospitalization [HR = 1.13, CI = (0.78, 1.62), p-value = 0.52].There was also no evidence that low exposure [HR = 1.17, CI = (0.79,1.73), p-value = 0.43] or high exposure [HR = 1.03, CI = (0.54, 1.99), p-value = 0.92] to any DMD was associated with death during hospitalization.Comorbidity was once again strongly associated with mortality, and, relativeto an individual with no comorbidities, a comorbidity of one [HR = 1.54, CI= (1.15, 2.04), p-value < 0.001] had a lower HR than a comorbidity of two[HR = 2.56, CI = (1.99, 3.30), p-value < 0.001]. The results were similarwhen using a secondary measure of exposure. Once more, this was a muchstronger effect than was shown in the time-to-death model of Chapter 2.3.7 Discussion of intensity-based analysisExposure had a time-dependent effect in three transitions.1. In the O → H transition, exposure to any disease-modifying drug(primary measure) was initially associated with fewer hospitalizations.This effect changed over time and the HR of exposure increased. At10.4 years post-index, the HR was one. After 10.4 years, exposure wasassociated with more hospitalizations.High exposure was initially associated with a 20% reduction in thehazard but after 4000 days of follow-up, the HR increased.2. In the O → D transition, the effect of exposure on mortality outsidethe hospital exhibited a time-dependent effect as well. In the precedingchapter, there appeared to be a slight increase in mortality after 6000days post-index. When deaths outside the hospital are consideredseparately, this increase became sharper and seems to occur earlier ataround 4000 days post-index.643.7. Discussion of intensity-based analysis3. In the H → O transition, exposure was associated with a 6% decreasein the hazard of discharge during the first day. For the remainder ofthe hospital visit, exposure is not associated with discharge.It is hard to come up with a plausible explanation of the time-dependenteffects in the O → H and O → D transitions. Given that follow-up could be22 years long, it is not surprising that many of the covariate effects changedover time but it is still worthwhile to discuss why this may be the case.The hazard ratio has a “built-in selection bias” [37]. The HR at time t isestimated only using individuals who were at risk at time t. For the O → Htransition, this is individuals who have not died, were not hospitalized, andwere not censored at time t. The types of individuals becoming exposedchanges as well. Two very general explanations: exposure to DMDs re-ally is associated with more hospitalizations and a greater chance of deaththan being unexposed after about 10 years post-index; the “built-in selec-tion bias” causes an increase in the hazard ratio but the effect of exposurehas not changed.Herna´n calls individuals more likely to experience the event of interest (orevents in this case) “susceptible” [37]. This may be due to an unmeasuredconfounder that causes individuals to be susceptible or because of generalfrailty not related to any single variable. If the “built-in selection bias” isresponsible for the time-dependent effect, then later in the study, censoring,mortality, hospitalization, or exposure itself (and perhaps, some combinationof several) cause a greater proportion of susceptible individuals to remainamong the exposed group relative to the proportion of susceptible individ-uals that are unexposed. More exposed individuals are hospitalized or diebut this is due to susceptibility, which is hidden, and not from exposure.It seems highly unlikely that censoring would be responsible for the increasein the HR. Censoring occurred at the end of the study (December 31, 2017),when an individual emigrated from the province, or, for very few individu-als (< 6), when daclizumab or ocrelizumab was prescribed. Administrativecensoring and the majority of emigration can be assumed to be independentof hospitalization and death. Daclizumab or ocrelizumab prescription maynot be, but this applies to so few individuals and occurred only in the finalyear of the study. Thus, even if daclizumab and ocrelizumab prescriptionis associated with hospitalization and death, this probably does not mean-ingfully affect the HRs. Similarly, some emigration may be associated with653.7. Discussion of intensity-based analysishealth status but this is a relatively small proportion of censored individualsand emigration for health reasons is likely not very common.The hazard of hospitalization between exposed and unexposed individualsis equal at around the same time (about 10.5 to 11 years post-index) thatthe hazard of mortality for exposed individuals begins to increase relative tounexposed individuals. That individuals are hospitalized and thus, not ableto experience O → H or O → D transitions, is likely not responsible forthe time-varying effect of exposure since the fraction of individuals in thehospital at any one time is low (about 0.42% according to Figure 3.4), andhospitalizations tended to be short so hospitalized individuals were soon atrisk of an O → H or O → D transition once again.The competing risk of death may be responsible for the time-dependent na-ture of the hazard ratios [38]. For example, near the beginning of the study,perhaps a lower proportion of individuals susceptible to hospitalization diedin the exposure group rather than the unexposed group. This may be dueto exposure or some other unknown reason. Then, later on in the study, theexposed group would have a higher proportion of susceptible individuals sorelatively more exposed individuals would be hospitalized compared to theunexposed group. A similar argument could be made for the O → D transi-tion about individuals likely to die. There is no evidence for this argumentbut it is not possible to disprove either.The exposure group did change over time but this does not seem to explainthe time-dependent HRs either. The time-to-exposure model indicated thatyounger individuals, individuals with no comorbidities, and individuals withfewer pre-index hospitalizations were exposed sooner (see Table B.8) butthe Cox models adjust for these covariates. It is possible that the exposuregroup changed over time with respect to some unmeasured covariate thatcauses “susceptibility” and that this was responsible for a higher hazard ofhospitalization, but this does not explain the changes in the HR that well.Of the individuals that were exposed during their follow-up, over 80% wereexposed by 5 years post-index and over 90% were exposed by 7 years post-index. If the characteristics of the exposed group were responsible for thetime-dependent effect, one would expect to see the changes occur as thegroups were changing, i.e., an increase in the HR of mortality earlier andthe time-dependent effects to stabilize at some point.Ng et al. [2] reported that current exposure (defined as a prescription in the663.7. Discussion of intensity-based analysisprevious 90 days) was associated with a 24% lower hazard of hospitalizationand no time-dependent effect on hospitalization of exposure to any DMD[HR = 0.76, CI = (0.71, 0.82)] which is consistent with the HR for exposureestimated in this study near the beginning of follow-up. This may be fora few reasons. First, as stated in the discussion of the preceding chapter,this study uses a subset of the cohort in Ng et al’s study. Second, pre-indexhospitalizations was not a covariate in Ng et al.’s model. Finally, the twostudies used different models to analyze the hazard of hospitalizations—Nget al. used a proportional means model and the covariates were adjusteddifferently.DMDs have been shown to reduce the frequency and severity of relapseswhich Ng et al. suggest may be the reason that DMD exposure was asso-ciated with fewer hospitalizations [2]. This is not necessarily the case forall stages of MS, and the efficacy of DMDs decreases during more advancedstages which may explain why the HR changed over time [12]. Discontinu-ation of DMDs may also be associated with a rebound of the disease [12].Although this varies from person to person, it could contribute to the in-crease in the HR of hospitalization and mortality outside of the hospital.If DMDs are not discontinued, negative reactions and adverse events are aparticular issue in older individuals due to age-related diseases and greaterdisability [39]. Ng et al. reported several studies with similar conclusionsbut noted that there were few comparable studies and did not report anythat investigated hospitalizations over the long-term. A systematic reviewpublished in 2020 [28] found that many studies on the long-term outcomesassociated with DMD use are biased. This emphasizes the need for the long-term observational studies such as this one but also indicates that there arefew studies that allow direct comparison with our results.67Chapter 4Analyzing the effect ofexposure on the mean timespent in and out of thehospital using pseudo-values4.1 Mean sojourn timeThe mean sojourn time in state k from time s to time t denoted ψk(s, t) isan population average, ignoring covariates.Consider an individual who was at risk in any state at time s and letI(Z(u) = k) be an indicator function that equals one when this individ-ual is in state k at time u and zero otherwise. The sojourn time between sand t in state k for this individual is∫ tsI(Z(u) = k) du.The mean sojourn time in state k during [s, t) among all individuals isψk(s, t) = E[∫ tsI(Z(u) = k) du]=∫ tsE[I(Z(u) = k)]du =∫ tsQk(u) du ,where Qk(u) = P (Z(u) = k) as defined in Section 3.2.3.Estimates of the mean sojourn time are based upon Q̂k(t). A curve of esti-mated state occupancy probabilities may be constructed by treating Q̂k(t)as a function of time t. This is a step function since estimated probabilitiesonly change at transition times. ψk(s, t) may be estimated by the area underthe curve (AUC) of Q̂k treated as a function of time from s to t [17], [33].This estimate is denoted ψ̂k(s, t).684.2. Pseudo-values4.2 Pseudo-valuesPseudo-value regression may be used to estimate the effects of the covari-ates on the mean time in a state (e.g., ψH(s, t), the mean time spent in thehospital between time s and t) [17]. The benefit of pseudo-value regressionover comparing the AUCs directly is that the effects of exposure and comor-bidity on ψ̂k(s, t) are easily adjusted for other covariates measured at time s.Let Sik be the sojourn time in state k for individual i when there is nocensoring. Then the sojourn times for every individual would be observed,and the mean sojourn time in state k denoted S¯k could be easily determinedby taking the average of the sojourn times for all n individuals, where n isthe number of individuals alive and uncensored at time s. If S¯(−i)k is themean sojourn time for all individuals except for individual i in the case ofno censoring, note that Sik may be rewritten asSik = nS¯k − (n− 1)S¯(−i)k .When individuals are censored, Sik is not known for some individuals, andthe mean sojourn time in state k from s to t may be estimated with ψ̂k(s, t),which is defined in the presence of censoring. The foregoing equation maybe used to estimate Sik. Denote the estimate ψ̂k(s, t) when individual i isexcluded from the dataset as ψ̂(−i)k (s, t). Then, the pseudo-value of ψ̂k(s, t)for individual i is:ψ˜ik(s, t) = n ∗ ψ̂k(s, t)− (n− 1) ∗ ψ̂(−i)k (s, t).This pseudo-value is then utilized as the response variable in a regressionmodel. Time-dependent covariates may be handled with landmarking: theeffect of the time-dependent covariate is assumed to depend on its value atlandmark time s [40]. This approach allows an investigation of the effect oftime-dependent exposure on the expected time spent in a state but there aretwo major limitations. First, as is the case with state occupancy probabili-ties, in the presence of a time-dependent exposure, pseudo-value regressionmay only explore associations between sojourn times and exposure at a cer-tain time. One solution is to adapt for pseudo-values the strategy articulatedby Cook and Lawless for state occupancy probabilities and time-dependentexposure [34]—that is, to calculate pseudo-values for different intervals [s, t)to explore the relationship between state sojourn times during [s, t) and ex-posure at different times s.694.2. Pseudo-valuesSecond, individuals exposed after s but before t are grouped with individualsunexposed between [s, t). It is unavoidable that the true effect of exposure,if any exists, will be slightly obscured, but judicious choices of s and t mayminimize this. The estimated effect of exposure is likely to be understatedsince, if anything, some of the effect of exposure will be experienced by the‘unexposed’ group. As a result, a large estimated effect of exposure may besafely interpreted as such, and a small effect of exposure, or cases where itis hard to tell, may be cautiously interpreted as an effect of exposure.For our own cohort, time spent hospitalized was small compared to timespent out of hospital. The average time spent out of the hospital in thestudy was 6,518 days and the average time spent hospitalized was 31 days.This was true even among individuals with at least one pre-index hospital-ization: 6,479 days versus 37 days. It may be the case that changes in theaverage sojourn times in state H associated with covariates may be too smallto be identified by this method. Grand and Putter [41] used landmarkingand pseudo-values to investigate the effect of a time-dependent covariate onthe mean sojourn time in a reversible illness-death model. Follow-up wasalso long-term in that study (often more than a decade) but one major dif-ference compared to this study was that years could be spent in the illnessstate, and for some individuals, more time was spent in the illness state thanthe healthy state.Individual pseudo-values do not have an intuitive interpretation. Andersenand Perme [42] showed the possible values that a pseudo-value may take intwo examples from the basic survival context. When censoring was higherthan 50%, pseudo-values for the one-year restricted mean lifetime were neg-ative for individuals that experienced the event before one year and aboveone for individuals that were censored at one year. Both are impossiblevalues. A similar phenemenon arises in pseudo-values for the Kaplan-Meiersurvival estimator. Even though Kaplan-Meier estimates are between 0 and1, pseudo-values for individuals that are still alive all have the same valueand are above 1. Pseudo-values for censored individuals may be above 1 orbetween 0 and 1, depending on when they were censored, and a pseudo-valuebecomes negative after an individual dies.Pseudo-values for the mean sojourn time in a multi-state model are morecomplicated but the general idea that pseudo-values should not be inter-preted individually is still true. Pseudo-values for the mean time spent instate O were frequently negative, and individuals who died early on and/or704.3. Pseudo-value regression modelshad frequent and long hospitalizations had low pseudo-values. Individualswith more time spent in the hospital had larger pseudo-values for the meansojourn time in state H. Pseudo-values for mean time spent in state H areoften unusually high, e.g., over 4000 in the case of the largest residual. Thesepseudo-values tended to be calculated from individuals that experienced oneor several of the 259 hospitalizations longer than 90 days. If few people werein the hospital during this period, such as near the end of follow-up, thenthe pseudo-values were even larger.That pseudo-values can take on impossible values for the mean sojourntime may seem problematic but is not actually an issue. Regression modelsthat have pseudo-values as the response variable provide interpretable andapproximately unbiased coefficient estimates. As Andersen and Ravn [43]note generally for AJ-based pseudo-values and Overgaard et al. [44] notefor the mean sojourn time specifically, under the assumption of independentcensoringE[ψ˜ik(s, t) | Xi(s)] ≈ E[ψik(s, t) | Xi(s)],where Xi(s) is a vector of possibly time-dependent covariates at time s. Ifcensoring is not independent, then adjustments to ψ˜ik(s, t) are necessary.As specified in the discussion of the previous chapter, censoring is assumedindependent here (or dependent censoring is rare) so no adjustments aremade.4.3 Pseudo-value regression modelsThe pseudo-value for the mean sojourn time is the response in a linearregression model even though some of the assumptions of linear regressionare violated: the pseudo-values are not independent and the errors seemto have unequal variances. The dependence between the pseudo-values isnumerically small enough that it is ignored in practice, and that is how wewill proceed for the pseudo-value regression models in this section.Heteroscedascity-consistent standard errorsUnequal variances still persisted when the pseudo-values were transformedusing log or Box-Cox transformations, so our inferences for these regres-sion models use the heteroscedascity-consistent (HC) variance estimator, anasymptotically consistent estimator of the variance in the presence of het-eroscedascity [45]. Let R be an n × n diagonal matrix where the entry in714.3. Pseudo-value regression modelsrow/colomn i is the square of the ith residual. The HC estimator is(XTX)−1XTRX(XTX)−1.Since this estimator is based upon asymptotic theory, using it can be prob-lematic when the sample size is small. This should not be a concern heresince the pseudo-value regression models in this study have a sample size≥ 4000 individuals. However, the HC variance estimate may still be bi-ased in large sample sizes, especially if there are influential points [45], [46].Flaichaire [47] verified with simulations that issues with the HC estimatorcan exist even when n = 1000.An alternative approach is SEs obtained via the wild bootstrap [45], [47].The wild bootstrap is an extension of the residual bootstrap, a standardbootstrap technique. The most common type of bootstrap, the pairs boot-strap, samples n observations from the dataset with replacement. The resid-ual bootstrap on the other hand modifies every observation from the originaldataset for each iteration and then uses the modified sample to calculatethe SE. The modified response for an observation is obtained by sampling aresidual with replacement (from the set of n residuals of the model fit), andthen adding that residual to the fitted value for the observation.The wild bootstrap modifies the original dataset using a method similarto the residual bootstrap [47], [48]. The wild bootstrap procedure usingpseudo-values is as follows. Let ψ˜∗ik(s, t) be the fitted pseudo-value frompseudo-value regression model k specified below. Let ri be the residual fromthe ith observation, and let ai be a random variable with a distribution suchthat E[ai] = 0 and E[a2i ] = 1. The ith response for a wild bootstrap sampleis yBi = ψ˜∗ik(s, t) + riai. The original model is refit with the same X matrixbut using yBi , i = 1, . . . , n as the response.There are a few popular choices for the distribution of ai. We use thefollowing:fa ={−(√5− 1)/2 with probability (√5 + 1)/(2√5)(√5 + 1)/2 with probability 1− (√5 + 1)/(2√5) .fa is the most popular choice for the wild bootstrap and handles errors withunequal variances that are asymmetrically distributed, although it may notbe the best choice generally [47].724.3. Pseudo-value regression modelsThe HC standard errors were checked against the wild bootstrap standarderrors using 5,000 bootstrap samples. If they were similar, we proceededunder the assumption that there was little to no bias in the HC estimates.Thus, we address the possibility of unequal variances with the HC standarderrors (or wild bootstrap).Choices of s and tFor this cohort, it may be difficult to use pseudo-values when s < 2 yearssince many individuals are becoming exposed during the first two years af-ter the index date. Of the individuals ever exposed, 23% met the definitionof exposed in the first year after the index date, and a further 34% wereexposed between the first and second year post-index. However, the rateof individuals becoming exposed slows down considerably after 2 years andcontinues to slow as time goes on. After two years post-index, 57% of in-dividuals ever exposed are exposed, and this rises to 83% in the next threeyears (five years post-index). By the end of year seven, this has risen to91% and by year ten, it is 95%. Four intervals of time may thus be useful:[2 years, 4 years), [4 years, 7 years), [7 years, 10 years), [10 years, 22 years).The association between a time-dependent exposure and the number of daysspent out of the hospital and alive between time s and time t may be esti-mated by the coefficients in a pseudo-value regression that uses ψ˜iO(s, t) asthe response variable. A pseudo-value regression model for four time periodscan investigate these effects over the study:• Model O for time period 1 uses ψ˜iO(2, 4).• Model O for time period 2 uses ψ˜iO(4, 7).• Model O for time period 3 uses ψ˜iO(7, 10).• Model O for time period 4 uses ψ˜iO(10, 22).Similarly, the impact of time-dependent exposure and comorbidity on thenumber of days spent in the hospital between time s and time t may beestimated by the coefficients in a pseudo-value regression that uses ψ˜iH(s, t)as the response variable. A pseudo-value regression model for four timeperiods can also investigate these effects over the study:• Model H for time period 1 uses ψ˜iH(2, 4).734.3. Pseudo-value regression models• Model H for time period 2 uses ψ˜iH(4, 7).• Model H for time period 3 uses ψ˜iH(7, 10).• Model H for time period 4 uses ψ˜iH(10, 22).Model specificationAge at the index date was modelled by a natural spline with three degreesof freedom, which is represented by three basis functions [26]. Let thesebasis functions be denoted h1(·), h2(·), and h3(·). This appeared to explainthe relationship between age and the pseudo-values sufficiently—more d.f.and there was often too much noise in the estimate while fewer d.f. did notappear to capture the relationship well. The number of pre-index hospital-izations were modelled with a log(x+ 1) transformation.Model k for the time period of [s, t) has the form:ψ˜ik(s, t) = β0 + β1Exposurei(s) + β2Sexi +β3h1(Agei) + β4h2(Agei) + β5h3(Agei) +β6SES2i + β7SES3i + β8SES4i + β9SES5i +β10Comorbidity1i(s) + β11Comorbidity2i(s) + β12 log(Hospitalizationsi + 1) +β13Exposurei(s) ∗ Comorbidity1i(s) + β14Exposurei(s) ∗ Comorbidity2i(s) + ϵi,where ϵi is an error term.In Model H, the coefficients may be interpreted in terms of the number ofdays spent and alive hospitalized. For example, if β2 > 0, on average menspend β2 more days in the hospital than women. Similarly, in Model O,the coefficients may be interpreted in terms of the number of days spentalive and out of the hospital. In each of the models, β1, β13, and β14 are ofprimary interest for answering the research questions and will be discussedin the main body.Model assessmentFor each model and time period, there is a separate appendix section (Sec-tions 1 to 8 of Appendix D) that includes:• A table summarizing the model and its coefficients. Confidence inter-vals and hypothesis tests of no effect are based upon the HC standarderrors.744.4. Results from the pseudo-value models• Within the table is a comparison of HC standard errors and wild boot-strap standard errors. The two methods produced very similar vari-ances for all terms in both pseudo-value regression models in all of thetime periods.• A plot showing the shape of the spline for age with confidence intervalsbased upon the HC standard errors.• Residual plots.Residual plots against fitted values (for example, Figure D.2 in the appendix)indicate that the data are heteroscedastic. In the presence of unequal vari-ances, the expected value of the error terms is still zero if the structuralpart of the regression model is correctly specified. Heteroscedasticity makesassessing the fit of the models using residual plots difficult, and scatter-plot smoothers are necessary to discern any patterns [42]. For the residualplots in this study, a natural spline with five d.f. was fit to the residuals.If the residuals of a pseudo-value regression model are centred at zero, thesmoother applied to a plot of residuals versus fitted values should be at ap-proximately zero. This suggests the structural component of the regressionmodel is approximately correctly specified.Deviations from zero indicate that there may be an issue with the model.In all eight pseudo-value regression models described below, there are afew fitted pseudo-values that are very small or large, and the scatterplotsmoothers can be overly sensitive to these extreme fitted pseudo-values.The smoother may drift away from zero even if there is no serious issue withthe fit of the pseudo-value regression model [42], [43]. Plots of residualsversus age and residuals versus pre-index hospitalizations may be used toprovide an additional check on the structural component and the functionalform specified for age and pre-index hospitalizations. Other more complexgoodness-of-fit methods for pseudo-value regression models exist but are notwidely used [49].4.4 Results from the pseudo-value modelsThe estimated effects of exposure on the time spent in State O and StateH may be found in Tables 4.1 and 4.2. The coefficients are in terms ofdays. For comparison, ignoring censoring, the mean sojourn time in stateO in time periods 1 to 4 was 721 days, 1076 days, 1074 days, and 4048 daysrespectively. The mean sojourn time in state H (ignoring censoring) in time754.4. Results from the pseudo-value modelsPseudo-value Regression for Mean Sojourn Time in State OCoefficient CI p-valueTime Period 1, 2 to 4 yearsExposure 0.54 (-1.59, 2.67) 0.62Exposure * Comorbidity1 4.39 (-2.07, 10.84) 0.18Exposure * Comorbidity2 16.85 (7.84, 25.85) < 0.001Time Period 2, 4 to 7 yearsExposure -0.68 (-4.06, 2.70) 0.69Exposure * Comorbidity1 -3.23 (-21.31, 14.84) 0.73Exposure * Comorbidity2 54.26 (33.98, 74.53) < 0.001Time Period 3, 7 to 10 yearsExposure -3.64 (-8.49, 1.22) 0.14Exposure * Comorbidity1 22.15 (4.41, 39.89) 0.014Exposure * Comorbidity2 28.22 (-4.27, 60.72) 0.089Time Period 4, 10 years to end of studyExposure -71.08 (-144.16, 1.99) 0.057Exposure * Comorbidity1 82.91 (-99.39, 265.21) 0.37Exposure * Comorbidity2 261.38 (-0.98, 523.74) 0.051Table 4.1: Coefficients for exposure and interaction terms in pseudo-valueregression model Operiods 1 to 4 was 3.0 days, 4.7 days, 4.5 days, and 21.0 days respectively.For the first ten years of follow-up, exposure does not appear to be associ-ated with the mean time in state O for individuals with no comorbidities.After the first decade, exposure seemed to be associated with over 70 fewerdays spent in State O [-71.08, CI = (-144.16, 1.99), p-value = 0.057].Effects of exposure emerged primarily for individuals with a comorbiditylevel of two. Compared to exposed individuals with no comorbidities, ex-posure increased the mean sojourn time in state O for individuals with acomorbidity of two in the following way: between two and four years post-index, exposure was associated with 17 more days spent in state O [16.85,CI = (7.84, 25.85), p-value < 0.001]; between four and seven years post-index, it was associated with 54 extra days [54.26, CI = (33.98, 74.53),764.4. Results from the pseudo-value modelsp-value < 0.001]; from 10 years post-index to the end of the study, it ap-peared to be associated with 261 more days in state O [261.38, CI = (-0.98,523.74), p-value = 0.051].Between seven and 10 years post-index, there was less evidence of an effectfor this comparison [28.22, CI = (-4.27, 60.72), p-value = 0.089]. Given thatthe estimated effect is still quite large and there was a consistent effect in thefirst, second, and fourth time periods, this is a case when we are cautiouslyinterpreting the result as an effect. In this same time period that, the modeldid suggest that, relative to exposed individuals with no comorbidity, expo-sure was associated with more time spent in state O for individuals with acomorbidity of one [22.15, CI = (4.41, 39.89), p-value = 0.014].In Model O, the positive interaction effects between exposure and comor-bidity offset some of the negative effects of comorbidity.The effect of exposure on the mean sojourn time in State H is less obvious.Relatively little follow-up time was spent in state H so the effect sizes aremuch smaller. Nevertheless, in Model H for Time Period 2, exposure wasassociated with reducing the number of days spent in hospital by 6.5 forindividuals with a comorbidity of two compared to exposed individuals withno comorbidities.In Time Period 4, exposure was associated with 6.5 more days spent in thehospital. The interaction terms between exposure and comorbidity were notsignificant. However, the effect sizes were quite large which may indicatethat some association exists but was not fully captured by the pseudo-valuemodel.It is clear from Table 1 and Table 2 that the pattern of association withexposure in Model O is not commensurate with that in Model H. Whenthere was an increase in the days spent in state O, this did not necessarilycorrespond to a large decrease in time spent in State H which suggests thatthere was fewer days spent in state D.Section 9 of Appendix D has a short summary of the effects of other covari-ates.774.4. Results from the pseudo-value modelsPseudo-value Regression for Mean Sojourn Time in State HCoefficient CI p-valueTime Period 1, 2 to 4 yearsExposure 0.16 (-0.66, 0.98) 0.70Exposure * Comorbidity1 -0.60 (-3.64, 2.45) 0.70Exposure * Comorbidity2 -1.86 (-4.17, 0.46) 0.12Time Period 2, 4 to 7 yearsExposure 0.39 (-0.62, 1.39) 0.45Exposure * Comorbidity1 -1.82 (-5.78, 2.13) 0.37Exposure * Comorbidity2 -6.51 (-10.67, -2.34) 0.002Time Period 3, 7 to 10 yearsExposure 1.00 (-0.17, 2.17) 0.094Exposure * Comorbidity1 -2.28 (-6.24, 1.67) 0.26Exposure * Comorbidity2 -0.97 (-7.45, 5.51) 0.77Time Period 4, 10 years to end of studyExposure 6.54 (0.86, 12.21) 0.024Exposure * Comorbidity1 6.94 (-8.82, 22.69) 0.39Exposure * Comorbidity2 -9.22 (-26.76, 8.31) 0.30Table 4.2: Coefficients for exposure and interaction terms in pseudo-valueregression model H784.5. Discussion4.5 DiscussionComorbidity complicates treatment decisions about MS [50]. Individualswith comorbidities face worse outcomes. A greater comorbidity burdentended to be associated with fewer days spent in state O (as shown in Sec-tions 1 to 4 of Appendix D) and often was associated with more days spentin state H (Sections 5 to 8 of Appendix D). This is consistent with theintensity-based models in the previous chapter that investigated hospital-izations and death where comorbidity was prognostic and associated withhospitalization (in models where it was not used as a stratifying variable).This chapter demonstrated that exposure to any disease-modifying therapywas consistently associated with a benefit (more time in the “healthy” state)to individuals with a comorbidity of two relative to individuals with no co-morbidity. It is not obvious why this is the case. The effect of exposure wassometimes large enough to overcome the burden of comorbidities. Disease-modifying drugs have been shown to induce new comorbidities but it is alsopossible that they reduce the impact of existing comorbidities or preventnew comorbidities from developing [51].Compared to their unexposed counterparts, individuals with a comorbidityof two that are exposed may have unmeasured qualities that made theirphysicians more likely to prescribe any DMD to them in the first place orqualities that made it more likely they adhered to treatment. Such qualitiesmay also increase the likelihood of staying in state O, so any change in meansojourn times may not be a result of exposure.Comparison of these results to other studies is difficult. The discussion of Nget al. [2] includes a brief review of previous work on the connection betweendisease-modifying drugs and hospitalization. Several studies were outlinedbut focused on the frequency of hospitalization; none looked at time spentin or out of the hospital. Multi-state models have been applied to multiplesclerosis data before but these studies investigated disability, time spent atdifferent levels of disability, and the time to disability progression [52]–[54].79Chapter 5ConclusionThis thesis applied multi-state modelling to British Columbia health ad-ministrative data to answer four questions about the efficacy of disease-modifying drugs for multiple sclerosis. The primary measure of exposurewas defined based upon the minimum time thought to produce a clinicaleffect.Does the effect of DMDs change over time?The effects of DMD exposure often varied with time, which is not surprisinggiven the study was 22 years long.Are disease-modifying drugs associated with a decrease in the frequency ofhospitalization or shorter hospital visits?The intensity-based multi-state analysis showed that within the first tenyears, exposure to DMDs was associated with fewer hospitalizations. As thestudy went on, exposure became associated with more hospitalizations.In the intensity-based analysis, exposure did not appear to shorten hospi-talizations, but the pseudo-value regression models indicated that exposurewas associated with more time spent in the hospital later in the study.Are disease-modifying drugs associated with a decrease in mortality?Results from the survival analysis or intensity-based multi-state analysissuggested that exposure was either not associated with mortality or per-haps associated with a higher chance of dying. As time went on, the hazardratio of death for exposure appeared to increase.Does an individual’s comorbidity burden change either of the preceding as-sociations?805.1. Advantages of the multi-state approachThe pseudo-value regression models indicated that exposure was associatedwith more time spent alive and out of the hospital for individuals with ahigher comorbidity burden throughout follow-up. Exposed individuals witha CCI of two or more consistently experienced greater benefits associatedwith exposure than individuals with a smaller CCI level.Except for an increase in the chance of being discharged between 1 and 90days after admission to the hospital (the second set of coefficients in themodel for the H → O transition), this relationship was undetected in theintensity-based models.5.1 Advantages of the multi-state approachThe illness-death model delineates deaths inside and outside the hospitalwhereas a survival analysis considers them together. The results of theO → D and H → D Cox models of Chapter 3 had subtle differences fromthe survival models of Chapter 2, but there was a marked time-varying effectof DMD exposure in the O → D model that was not as pronounced in thetime-to-death model.The illness-death model revealed, in contrast to what could be known aftera recurrent event analysis of hospitalizations, differences between the effectsof exposure on the O → H and H → O transition. This suggested that theeffect of exposure on the frequency of hospitalizations does not translate toan effect on the length of hospitalizations.The risk of an event is often of greater interest than the hazard ratio, butit is reported infrequently in observational studies [37], probably becauseit is more difficult to calculate. The first obstacle is the competing risk ofdeath. For example, in the DRUMS data, the risk of hospitalization cannotbe properly estimated without considering the competing risk of death [55].This additional complexity is incorporated into multi-state models duringspecification of the state space. By including a state for death (in additionto a state for hospitalization), the risk of hospitalization may be estimatedusing AJ-based estimates of state occupancy probability. Two other prob-lems that arise when reporting the effect of exposure on risk is adjusting theestimate for other covariates and handling censoring. Both are managed inmulti-state models with pseudo-value regression models. As demonstratedin Chapter 4, they also provide intuitive interpretations of covariate effects.815.2. Limitations of our studyWithout an analysis of absolute risks, the apparent benefits of DMD expo-sure for individuals with high comorbidity burdens would have been over-looked. On the other hand, the HRs indicate why an individual makescertain transitions between states, and so, the time-varying nature of DMDexposure on hospitalizations was clearest in the O → H Cox model. Multi-state models provide a framework of incorporating both risk and hazardratio estimates to gain insights into multiple sclerosis.5.2 Limitations of our studyThe administrative health records do not capture some important clinicalinformation about the individuals. For example, clinical studies on MS mayrecord the level of disability that subjects have, as measured by the Ex-panded Disability Status Scale (EDSS), or their specific MS subtypes. Brainatrophy takes place in different locations, severity, and frequency dependingon MS subtype [56], and the extent of MS progression is commonly measuredusing EDSS [57]. This analysis instead could only adjust for MS progressionvia the time since the index date, and it is possible that including EDSS andMS subtype as time-dependent covariates could produce different results forthis study. Moreover, relying on pre-index hospitalizations to capture thehealth of an individual—while useful given its predictive power in the mod-els of this study—does not fully describe an individual’s well-being.Access to healthcare can be affected by many possible variables that are in-dividual, structural, or systemic [58] that we do not consider. These factorsmay also affect an individual’s desire to utilize healthcare services, even ifaccess is not an issue. The primary outcome of this thesis is rates of hospi-talization but there are complex dynamics not considered in this paper thatmay affect the number of hospitalizations or the duration of hospital visits.Notwithstanding socioeconomic status, most of the covariates are probablynot relevant to structural or systemic barriers to healthcare.Other observational studies may have access to information about the ad-herence to treatment among subjects in the study and define measures ofexposure accordingly. This study had only drug dispensation informationand assumed that drugs dispensed were adhered to fully. Studies using in-formation collected at clinics would have a smaller sample size and be lesspowerful but may better adjust the results of Cox and pseudo-value regres-825.3. Future worksion. As a result, the findings of this paper should be interpreted withcaution.5.3 Future workMulti-state models are a flexible tool to answer questions about multiplesclerosis but previous work that has applied multi-state modelling to MSfocused on disability and had state spaces defined by EDSS score [52]–[54].Future studies may apply multi-state models to a wider variety of problems.For instance, this thesis investigated the efficacy of DMDs in the presence ofcomorbidities but disease-modifying drugs and comorbidity have a ‘bidirec-tional relationship’ [51] which could be explored further through multi-statemodels. Two examples of multi-state models that could investigate thisrelationship follow. Comorbidity burden could be used to define a statespace (e.g., three states for the three levels of comorbidity in this study)and the impact of DMD exposure on comorbidity burden could be inves-tigated through intensity-based or pseudo-value analysis. One could alsodefine the states using comorbidity burden and prescription status (i.e.,prescribed DMDs vs. not prescribed), and the bidirectional relationship be-tween DMD use and comorbidity could be investigated by exploring howindividuals move through these states.83Bibliography[1] H. S. Ng, F. Zhu, E. Kingwell, et al., “Disease-modifying drugs formultiple sclerosis and association with survival,” Neurology Neuroim-munology & Neuroinflammation, vol. 9, no. 5, e200005, 2022. doi:10.1212/NXI.0000000000200005. eprint: https://www.neurology.org/doi/pdf/10.1212/NXI.0000000000200005. [Online]. Avail-able: https : / / www . neurology . org / doi / abs / 10 . 1212 / NXI .0000000000200005.[2] H. S. Ng, F. Zhu, E. Kingwell, et al., “Disease-modifying drugs formultiple sclerosis and subsequent health service use,” Multiple Scle-rosis Journal, vol. 28, no. 4, pp. 583–596, 2022, PMID: 34949130.doi: 10.1177/13524585211063403. eprint: https://doi.org/10.1177/13524585211063403. [Online]. Available: https://doi.org/10.1177/13524585211063403.[3] R. Marcus, “What Is Multiple Sclerosis?” JAMA, vol. 328, no. 20,pp. 2078–2078, Nov. 2022, issn: 0098-7484. doi: 10 . 1001 / jama .2022.14236. eprint: https://jamanetwork.com/journals/jama/articlepdf/2798738/jama\_marcus\_2022\_pg\_220027\_1668792283.21183.pdf. [Online]. Available: https://doi.org/10.1001/jama.2022.14236.[4] B. S. Travers, B. K.-T. Tsang, and J. L. Barton, “Multiple sclero-sis: Diagnosis, disease-modifying therapy and prognosis,” AustralianJournal for General Practitioners, vol. 51, no. 4, pp. 199–206, Apr.2022. [Online]. Available: https://www1.racgp.org.au/ajgp/2022/april/multiple-sclerosis-diagnosis-therapy-and-prognosis.[5] A. Finkelsztejn, “Multiple sclerosis: Overview of disease-modifyingagents,” Perspectives in Medicinal Chemistry, vol. 6, PMC.S13213,2014, PMID: 25336899. doi: 10.4137/PMC.S13213. eprint: https://doi.org/10.4137/PMC.S13213. [Online]. Available: https://doi.org/10.4137/PMC.S13213.84BIBLIOGRAPHY[6] A. Rae-Grant, G. S. Day, R. A. Marrie, et al., “Comprehensive sys-tematic review summary: Disease-modifying therapies for adults withmultiple sclerosis,” Neurology, vol. 90, no. 17, pp. 789–800, 2018. doi:10.1212/WNL.0000000000005345. eprint: https://www.neurology.org/doi/pdf/10.1212/WNL.0000000000005345. [Online]. Avail-able: https : / / www . neurology . org / doi / abs / 10 . 1212 / WNL .0000000000005345.[7] H. S. Ng, F. Zhu, E. Kingwell, et al., “Characteristics of a population-based multiple sclerosis cohort treated with disease-modifying drugsin a universal healthcare setting,” Expert Review of Neurotherapeutics,vol. 21, no. 1, pp. 131–140, 2021.[8] R. A. Marrie, A. Miller, M. P. Sormani, et al., “The challenge of co-morbidity in clinical trials for multiple sclerosis,” Neurology, vol. 86,no. 15, pp. 1437–1445, 2016. doi: 10.1212/WNL.0000000000002471.eprint: https://www.neurology.org/doi/pdf/10.1212/WNL.0000000000002471. [Online]. Available: https://www.neurology.org/doi/abs/10.1212/WNL.0000000000002471.[9] R. Marrie, N. Yu, J. Blanchard, S. Leung, and L. Elliott, “The ris-ing prevalence and changing age distribution of multiple sclerosis inmanitoba,” Neurology, vol. 74, no. 6, pp. 465–471, 2010.[10] J. M. Valderas, B. Starfield, B. Sibbald, C. Salisbury, and M. Roland,“Defining comorbidity: Implications for understanding health and healthservices,” The Annals of Family Medicine, vol. 7, no. 4, pp. 357–363,2009, issn: 1544-1709. doi: 10.1370/afm.983. eprint: https://www.annfammed.org/content/7/4/357.full.pdf. [Online]. Available:https://www.annfammed.org/content/7/4/357.[11] H. S. Ng, F. Zhu, E. Kingwell, et al., “Disease-modifying drugs formultiple sclerosis and association with survival,” Neurology Neuroim-munology & Neuroinflammation, vol. 9, no. 5, e200005, 2022. doi:10.1212/NXI.0000000000200005. eprint: https://www.neurology.org/doi/pdf/10.1212/NXI.0000000000200005. [Online]. Avail-able: https : / / www . neurology . org / doi / abs / 10 . 1212 / NXI .0000000000200005.[12] H.-P. Hartung, S. G. Meuth, D. M. Miller, and G. Comi, “Stoppingdisease-modifying therapy in relapsing and progressive multiple sclero-sis,” Current Opinion in Neurology, vol. 34, no. 4, pp. 598–603, 2021.85BIBLIOGRAPHY[13] L. Prosperini, S. Haggiag, S. Ruggieri, C. Tortorella, and C. Gasperini,“Stopping disease-modifying treatments in multiple sclerosis: A sys-tematic review and meta-analysis of real-world studies,” CNS drugs,vol. 37, no. 10, pp. 915–927, 2023.[14] A. Ostolaza, J. Corroza, and T. Ayuso, “Multiple sclerosis and aging:Comorbidity and treatment challenges,”Multiple Sclerosis and RelatedDisorders, vol. 50, p. 102 815, 2021.[15] G. Giovannoni, E. Southam, and E. Waubant, “Systematic review ofdisease-modifying therapies to assess unmet needs in multiple sclerosis:Tolerability and adherence,” Multiple Sclerosis Journal, vol. 18, no. 7,pp. 932–946, 2012.[16] B. Hicks, J. A. Kaye, L. Azoulay, K. B. Kristensen, L. A. Habel,and A. Pottegard, “The application of lag times in cancer pharma-coepidemiology: A narrative review,” Annals of Epidemiology, vol. 84,pp. 25–32, 2023, issn: 1047-2797. doi: https://doi.org/10.1016/j . annepidem . 2023 . 05 . 004. [Online]. Available: https : / / www .sciencedirect.com/science/article/pii/S104727972300090X.[17] R. J. Cook and J. F. Lawless, Multistate models for the analysis of lifehistory data. CRC Press, 2018.[18] T. M. Therneau and P. M. Grambsch,Modeling Survival Data: Extend-ing the Cox Model. New York, NY: Springer New York, 2000, isbn:978-1-4757-3294-8. doi: 10.1007/978- 1- 4757- 3294- 8. [Online].Available: https://doi.org/10.1007/978-1-4757-3294-8.[19] D. R. Cox, “Partial likelihood,” Biometrika, vol. 62, no. 2, pp. 269–276, Aug. 1975, issn: 0006-3444. doi: 10.1093/biomet/62.2.269.eprint: https://academic.oup.com/biomet/article-pdf/62/2/269/964098/62-2-269.pdf. [Online]. Available: https://doi.org/10.1093/biomet/62.2.269.[20] J. F. Lawless, Statistical models and methods for lifetime data. JohnWiley & Sons, 2011.[21] P. M. Grambsch and T. M. Therneau, “Proportional hazards tests anddiagnostics based on weighted residuals,” Biometrika, vol. 81, no. 3,pp. 515–526, Sep. 1994, issn: 0006-3444. doi: 10.1093/biomet/81.3 . 515. eprint: https : / / academic . oup . com / biomet / article -pdf/81/3/515/714158/81-3-515.pdf. [Online]. Available: https://doi.org/10.1093/biomet/81.3.515.86BIBLIOGRAPHY[22] E. Kingwell, E. Leray, F. Zhu, et al., “Multiple sclerosis: Effect of betainterferon treatment on survival,” Brain, vol. 142, no. 5, pp. 1324–1333, 2019.[23] C. Danieli, S. Cohen, A. Liu, et al., “Flexible modeling of the associa-tion between cumulative exposure to low-dose ionizing radiation fromcardiac procedures and risk of cancer in adults with congenital heartdisease,” American Journal of Epidemiology, vol. 188, no. 8, pp. 1552–1562, 2019.[24] M.-P. Sylvestre and M. Abrahamowicz, “Flexible modeling of the cu-mulative effects of time-dependent exposures on the hazard,” Statisticsin Medicine, vol. 28, no. 27, pp. 3437–3453, 2009.[25] M.-P. Sylvestre, M.-E. Beauchamp, R. P. Kyle, and M. Abrahamowicz,Wce: Weighted cumulative exposure models, R package version 1.0.2,2018. [Online]. Available: https://CRAN.R-project.org/package=WCE.[26] T. Hastie, R. Tibshirani, and J. H. Friedman, The Elements of Sta-tistical Learning: Data Mining, Inference, and Prediction. Springer,2009, vol. 2.[27] M. Abrahamowicz, M.-E. Beauchamp, and M.-P. Sylvestre, “Com-parison of alternative models for linking drug exposure with adverseeffects,” Statistics in Medicine, vol. 31, no. 11-12, pp. 1014–1030, 2012.[28] R. C. Lucchetta, L. P. Leonart, M. V. Goncalves, et al., “Reliability inlong-term clinical studies of disease-modifying therapies for relapsing-remitting multiple sclerosis: A systematic review,” PLoS One, vol. 15,no. 6, e0231722, 2020.[29] H. S. Ng, F. Zhu, Y. Zhao, et al., “Adverse events associated withdisease-modifying drugs for multiple sclerosis: A multiregional population-based study,” Neurology, vol. 102, no. 3, e208006, 2024.[30] J. G. Le-Rademacher, R. A. Peterson, T. M. Therneau, B. L. San-ford, R. M. Stone, and S. J. Mandrekar, “Application of multi-statemodels in cancer clinical trials,” Clinical Trials, vol. 15, no. 5, pp. 489–498, 2018, PMID: 30035644. doi: 10.1177/1740774518789098. eprint:https://doi.org/10.1177/1740774518789098. [Online]. Available:https://doi.org/10.1177/1740774518789098.[31] J. G. Le-Rademacher, T. M. Therneau, and F.-S. Ou, “The utilityof multistate models: A flexible framework for time-to-event data,”Current Epidemiology Reports, vol. 9, no. 3, pp. 183–189, 2022.87BIBLIOGRAPHY[32] L. E. Lovblom, L. Briollais, B. A. Perkins, and G. Tomlinson, “Mod-eling multiple correlated end-organ disease trajectories: A tutorial formultistate and joint models with applications in diabetes complica-tions,” Statistics in Medicine, vol. 43, no. 5, pp. 1048–1082, 2024.doi: https : / / doi . org / 10 . 1002 / sim . 9984. eprint: https : / /onlinelibrary.wiley.com/doi/pdf/10.1002/sim.9984. [On-line]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.9984.[33] C. R. Jack Jr, T. M. Therneau, E. S. Lundt, et al., “Long-term associa-tions between amyloid positron emission tomography, sex, apolipopro-tein E and incident dementia and mortality among individuals with-out dementia: hazard ratios and absolute risk,” Brain Communica-tions, vol. 4, no. 2, fcac017, Feb. 2022, issn: 2632-1297. doi: 10 .1093/braincomms/fcac017. eprint: https://academic.oup.com/braincomms / article - pdf / 4 / 2 / fcac017 / 42902090 / fcac017 .pdf. [Online]. Available: https://doi.org/10.1093/braincomms/fcac017.[34] R. J. Cook and J. F. Lawless, “Life history analysis with multistatemodels: A review and some current issues,” Canadian Journal of Statis-tics, vol. 50, no. 4, pp. 1270–1298, 2022. doi: https://doi.org/10 . 1002 / cjs . 11711. eprint: https : / / onlinelibrary . wiley .com / doi / pdf / 10 . 1002 / cjs . 11711. [Online]. Available: https ://onlinelibrary.wiley.com/doi/abs/10.1002/cjs.11711.[35] T. M. Therneau, A package for survival analysis in r, R package ver-sion 3.5-5, 2023. [Online]. Available: https://CRAN.R-project.org/package=survival.[36] P. Kragh Andersen, M. Pohar Perme, H. C. van Houwelingen, et al.,“Analysis of time-to-event for observational studies: Guidance to theuse of intensity models,” Statistics in Medicine, vol. 40, no. 1, pp. 185–211, 2021. doi: https://doi.org/10.1002/sim.8757. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.8757. [On-line]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.8757.[37] M. A. Herna´n, “The hazards of hazard ratios,” Epidemiology, vol. 21,no. 1, pp. 13–15, 2010.[38] M. A. Herna´n and J. M. Robins, Causal inference. Boca Raton: Chap-man Hall/CRC., 2020.88BIBLIOGRAPHY[39] E. M. M. Strijbis, A. Kerbrat, and J. R. Corboy, “Discontinuationof Disease-Modifying Therapy in Multiple Sclerosis: Should We Stayor Should We Go?” JAMA Neurology, vol. 78, no. 7, pp. 787–788,Jul. 2021, issn: 2168-6149. doi: 10.1001/jamaneurol.2021.0764.eprint: https : / / jamanetwork . com / journals / jamaneurology /articlepdf/2778356/jamaneurology\_strijbis\_2021\_vp\_210005\_1624467407.35917.pdf. [Online]. Available: https://doi.org/10.1001/jamaneurol.2021.0764.[40] H. Putter and H. C. van Houwelingen, “Understanding landmarkingand its relation with time-dependent cox regression,” Statistics in Bio-sciences, vol. 9, pp. 489–503, 2017.[41] M. K. Grand and H. Putter, “Regression models for expected lengthof stay,” Statistics in Medicine, vol. 35, no. 7, pp. 1178–1192, 2016.doi: https : / / doi . org / 10 . 1002 / sim . 6771. eprint: https : / /onlinelibrary.wiley.com/doi/pdf/10.1002/sim.6771. [Online].Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6771.[42] P. K. Andersen and M. P. Perme, “Pseudo-observations in survivalanalysis,” Statistical Methods in Medical Research, vol. 19, no. 1, pp. 71–99, 2010, PMID: 19654170. doi: 10.1177/0962280209105020. eprint:https://doi.org/10.1177/0962280209105020. [Online]. Available:https://doi.org/10.1177/0962280209105020.[43] P. K. Andersen and H. Ravn, Models for multi-state survival data:rates, risks, and pseudo-values. CRC Press, 2023.[44] M. Overgaard, P. K. Andersen, and E. T. Parner, “Pseudo-observationsin a multistate setting,” The Stata Journal, vol. 23, no. 2, pp. 491–517,2023. doi: 10.1177/1536867X231175332. eprint: https://doi.org/10.1177/1536867X231175332. [Online]. Available: https://doi.org/10.1177/1536867X231175332.[45] R. Davidson and E. Flachaire, “The wild bootstrap, tamed at last,”Journal of Econometrics, vol. 146, no. 1, pp. 162–169, 2008, issn: 0304-4076. doi: https://doi.org/10.1016/j.jeconom.2008.08.003.[Online]. Available: https://www.sciencedirect.com/science/article/pii/S0304407608000833.[46] A. Chesher and I. Jewitt, “The bias of a heteroskedasticity consistentcovariance matrix estimator,” Econometrica, vol. 55, no. 5, pp. 1217–1222, 1987, issn: 00129682, 14680262. [Online]. Available: http://www.jstor.org/stable/1911269 (visited on 04/27/2024).89BIBLIOGRAPHY[47] E. Flachaire, “Bootstrapping heteroskedastic regression models: Wildbootstrap vs. pairs bootstrap,” Computational Statistics Data Analy-sis, vol. 49, no. 2, pp. 361–376, 2005, 2nd CSDA Special Issue on Com-putational Econometrics, issn: 0167-9473. doi: https://doi.org/10.1016/j.csda.2004.05.018. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0167947304001537.[48] E. Mammen, “Bootstrap and wild bootstrap for high dimensional lin-ear models,” The Annals of Statistics, vol. 21, no. 1, pp. 255–285, 1993,issn: 00905364. [Online]. Available: http://www.jstor.org/stable/3035590 (visited on 04/27/2024).[49] K. Pavlicˇ, T. Martinussen, and P. K. Andersen, “Goodness of fit testsfor estimating equations based on pseudo-observations,” Lifetime DataAnalysis, vol. 25, pp. 189–205, 2019.[50] T. Zhang, H. Tremlett, S. Leung, et al., “Examining the effects ofcomorbidities on disease-modifying therapy use in multiple sclerosis,”Neurology, vol. 86, no. 14, pp. 1287–1295, 2016. doi: 10.1212/WNL.0000000000002543. eprint: https://www.neurology.org/doi/pdf/10.1212/WNL.0000000000002543. [Online]. Available: https://www.neurology.org/doi/abs/10.1212/WNL.0000000000002543.[51] R. A. Marrie, J. D. Fisk, K. Fitzgerald, et al., “Etiology, effects andmanagement of comorbidities in multiple sclerosis: Recent advances,”Frontiers in immunology, vol. 14, p. 1 197 195, 2023.[52] A. J. Palmer, I. van der Mei, B. V. Taylor, P. M. Clarke, J. SteveSimpson, and H. Ahmad, “Modelling the impact of multiple sclerosison life expectancy, quality-adjusted life years and total lifetime costs:Evidence from australia,” Multiple Sclerosis Journal, vol. 26, no. 4,pp. 411–420, 2020, PMID: 30806569. doi: 10.1177/1352458519831213.eprint: https://doi.org/10.1177/1352458519831213. [Online].Available: https://doi.org/10.1177/1352458519831213.[53] R. Gani, R. M. Nixon, S. Hughes, and C. H. Jackson, “Estimatingthe rates of disability progression in people with active relapsing-remitting multiple sclerosis,” Journal of Medical Economics, vol. 10,no. 2, pp. 79–89, 2007. doi: 10.3111/200710079089. eprint: https://doi.org/10.3111/200710079089. [Online]. Available: https://doi.org/10.3111/200710079089.90[54] N. Bergvall, H. Rathi, R. Nixon, H. Thom, J. Alsop, and L. Dun-sire, “Modeling the impact of disease modifying treatment on timeto disability health states in multiple sclerosis: An evaluation of oraltherapies through indirect comparisons of 6-month confirmed disabil-ity progression,” Value in Health, vol. 16, no. 7, A619, 2013.[55] P. K. Andersen, R. B. Geskus, T. de Witte, and H. Putter, “Competingrisks in epidemiology: possibilities and pitfalls,” International Journalof Epidemiology, vol. 41, no. 3, pp. 861–870, Jan. 2012, issn: 0300-5771. doi: 10.1093/ije/dyr213. eprint: https://academic.oup.com/ije/article-pdf/41/3/861/18481551/dyr213.pdf. [Online].Available: https://doi.org/10.1093/ije/dyr213.[56] E. Pagani, M. A. Rocca, A. Gallo, et al., “Regional brain atrophyevolves differently in patients with multiple sclerosis according to clin-ical phenotype,” American Journal of Neuroradiology, vol. 26, no. 2,pp. 341–346, 2005, issn: 0195-6108. eprint: https://www.ajnr.org/content/26/2/341.full.pdf. [Online]. Available: https://www.ajnr.org/content/26/2/341.[57] M. J. Hohol, E. J. Orav, and H. L. Weiner, “Disease steps in multiplesclerosis: A longitudinal study comparing disease steps and edss toevaluate disease progression,” Multiple Sclerosis Journal, vol. 5, no. 5,pp. 349–354, 1999, PMID: 10516779. doi: 10.1177/135245859900500508.eprint: https://doi.org/10.1177/135245859900500508. [Online].Available: https://doi.org/10.1177/135245859900500508.[58] N. Tzenios, “The determinants of access to healthcare: A review ofindividual, structural, and systemic factors,” Journal of Humanitiesand Applied Science Research, vol. 2, no. 1, pp. 1–14, Dec. 2019. [On-line]. Available: https://journals.sagescience.org/index.php/JHASR/article/view/23.[59] K. P. Burnham and D. R. Anderson, “Multimodel inference: Under-standing aic and bic in model selection,” Sociological Methods & Re-search, vol. 33, no. 2, pp. 261–304, 2004. doi: 10.1177/0049124104268644.eprint: https://doi.org/10.1177/0049124104268644. [Online].Available: https://doi.org/10.1177/0049124104268644.91Appendix AAppendix for Chapter 1Table A.1: Time of first prescriptionsDisease-modifying drug Month and year of first prescriptionBeta-interferon July, 1996Glatiramer acetate December, 1997Natalizumab March, 2007Fingolimod July, 2011Dimethyl fumarate June, 2013Teriflunomide March, 2014Alemtuzumab July, 2014Daclizumab June, 2017Ocrelizumab December, 201792Appendix BAppendix for Chapter 2B.1 Survival modelsB.1.1 Candidate modelsThe unstratified model has the highest AIC (14,581.66) followed by themodel stratified by index period (12,418.75) and lastly, the model stratifiedby calendar year of the follow-up segment (12,291.31). All three producedsimilar estimated HRs for all covariates (differences of less than 0.02). Thetentative final model was selected from these three candidate models basedupon the AIC.Table B.1: Unstratified survival modelHR LB UB pDMD1 1.0186 0.8031 1.2918 0.8793SEXM 1.9172 1.6851 2.1811 0.0000Index_Age 1.0924 1.0870 1.0979 0.0000SES2 0.7625 0.6291 0.9241 0.0057SES3 0.7677 0.6344 0.9289 0.0065SES4 0.6436 0.5292 0.7827 0.0000SES5 0.5682 0.4627 0.6977 0.0000Baseline_Comorbidity1 1.3618 1.1622 1.5958 0.0001Baseline_Comorbidity2 1.7556 1.4720 2.0938 0.0000Num_Prev_Hosp 1.0590 1.0421 1.0762 0.000093B.1. Survival modelsTable B.2: Survival model with fixed stratification by index periodHR LB UB pDMD1 1.0149 0.8001 1.2874 0.9030SEXM 1.9192 1.6867 2.1837 0.0000Index_Age 1.0926 1.0870 1.0981 0.0000SES2 0.7639 0.6301 0.9261 0.0061SES3 0.7699 0.6362 0.9317 0.0072SES4 0.6450 0.5302 0.7845 0.0000SES5 0.5675 0.4621 0.6969 0.0000Baseline_Comorbidity1 1.3623 1.1624 1.5967 0.0001Baseline_Comorbidity2 1.7554 1.4699 2.0964 0.0000Num_Prev_Hosp 1.0590 1.0420 1.0763 0.0000Table B.3: Model with time-dependent stratification by calendar periodHR LB UB pDMD1 1.0169 0.8016 1.2900 0.8904SEXM 1.9056 1.6748 2.1683 0.0000Index_Age 1.0928 1.0873 1.0984 0.0000SES2 0.7643 0.6305 0.9267 0.0062SES3 0.7720 0.6379 0.9342 0.0078SES4 0.6448 0.5301 0.7843 0.0000SES5 0.5720 0.4657 0.7025 0.0000Baseline_Comorbidity1 1.3732 1.1718 1.6092 0.0001Baseline_Comorbidity2 1.7594 1.4732 2.1012 0.0000Num_Prev_Hosp 1.0591 1.0422 1.0762 0.0000B.1.2 Final model for primary exposure measureFunctional formFigure B.1 and Figure B.2 show age at the index date and pre-index hospi-talizations respectively modelled using a natural spline with five degrees offreedom. There was no evidence that age had a non-linear relationship withthe hazard but the number of pre-index hospitalizations seemed to be mod-elled better using a log(x+ 1) function. This lowered the AIC to 12,276.41,so it was incorporated into an updated tentative final model.94B.1. Survival modelsFigure B.1: Age at index date modelled with a natural spline versus a linearfit in red.Figure B.2: Previous hospitalizations modelled with a natural spline versusa log(x+ 1) transformation in red.95B.1. Survival modelsProportional hazards assumptionFigure B.3: Scaled Schoenfeld residuals for the coefficient of the primaryexposure. Model estimate in red.The estimate of the coefficient of DMD exposure from the model is plottedin red as shown in Figure B.3. This is done to illustrate how a time-varyingeffect is diagnosed but is left out of subsequent scaled Schoenfeld residualplots. The coefficient of exposure was more or less constant. After t = 6000,the coefficient increased slightly and was outside of the confidence interval.Of the 965 individuals that died, 78 died after 6000 days of follow-up, andof that subset, only 18 died while exposed to DMDs. Estimates based on18 events should be considered skeptically so this time-dependent effect wasnot implemented in the model. As a final check to be sure that no time-dependent coefficient is necessary, a step function was used to create twocoefficients of the effect of exposure: one before t = 6000 and another after.This resulted in a model with an AIC of 12,274.99, a 1.42 reduction. Thisis a minor difference so the simpler model without the time-dependent co-efficient is used [59].The other coefficients in the model did not exhibit time-varying effects.Assessment of final model for primary exposure measureNo sum of scaled dfbeta residuals was larger than 1/3 in this updated fit.The largest sum for any coefficient was 0.26.96B.1. Survival modelsThe deviance residuals are shown in Figure B.4. The model seems to over-estimate survival times more than under-estimate them. This is reasonablesince all-cause mortality likely includes some deaths that are not predictablefrom the covariates such as accidental deaths. There are 34 deviance residu-als greater than three. These tended to be individuals that were young at theindex date (e.g., 20-39 years old) or died very soon into the study (t < 300days). There are 14 deviance residuals less than -1.75 which were olderindividuals that had comorbidities but lived a long time (mostly t > 6000).Figure B.4: Deviance residuals for the final model using the primary expo-sure measureResults of final model for primary exposure measureMost of the estimates of the HRs were as one would expect. The older anindividual was, the greater the hazard of death and generally, the higher SESan individual had, the lower their hazard of death. More hospitalizationsbefore the index date was likely an indication of worse overall health prior toMS diagnosis and was associated with a greater hazard. However, the effectof sex on the hazard of mortality was striking. Male individuals had a 93%greater hazard than female individuals, more than double the increase inhazard associated with a baseline comorbidity of one and greater still thanthe increase in mortality associated with a comorbidity of two.97B.1. Survival modelsTable B.4: Final model with primary exposure measureHR LB UB pDMD 1.0183 0.8025 1.2920 0.8814SEXM 1.9317 1.6973 2.1984 0.0000Index_Age 1.0911 1.0855 1.0966 0.0000SES2 0.7677 0.6332 0.9307 0.0071SES3 0.7729 0.6387 0.9353 0.0081SES4 0.6497 0.5341 0.7904 0.0000SES5 0.5789 0.4715 0.7109 0.0000Baseline_Comorbidity1 1.3289 1.1329 1.5588 0.0005Baseline_Comorbidity2 1.6562 1.3843 1.9815 0.0000log(Num_Prev_Hosp + 1) 1.4150 1.2888 1.5536 0.0000B.1.3 Final model for secondary exposure measureFigure B.5: Scaled Schoenfeld residuals for the coefficient of low exposure.98B.1. Survival modelsFigure B.6: Scaled Schoenfeld residuals for the coefficient of high exposure.The initial working model was specified similar to the final model for theprimary measure of exposure but with separate effects for low and high ex-posure. The functional form of age at index date and number of pre-indexhospitalizations was not re-examined.No individuals had a sum of scaled dfbeta residuals greater than 1/3; thelargest was 0.27.A natural spline with four degrees of freedom was fit to the scaled Schoenfeldresiduals of the Low exposure and High exposure coefficients (Figures B.5and B.6 respectively). Similar indications of possible time-dependent effectsthat were seen earlier for the primary definition of exposure now appear morepronounced for low exposure. In the final model, an interaction between thecoefficient of low exposure and t when t > 4000 days was included.Results of final model for secondary exposure measureThe estimates were close to identical between the primary and secondarymodels: men still had a hazard ratio that was nearly double that of women;higher socioeconomic status was associated with a reduction in the hazard;higher comorbidity levels were associated with a higher hazard; more pre-index hospitalizations was associated with a greater chance of death.99B.2. Cox model results using the cohort from Ng et al.Table B.5: Final model with secondary exposure measureHR LB UB pDMDLow 0.8167 0.5743 1.1615 0.2598DMDHigh 0.8742 0.5855 1.3053 0.5110SEXM 1.9058 1.6749 2.1685 0.0000Index_Age 1.0925 1.0870 1.0980 0.0000SES2 0.7657 0.6316 0.9283 0.0066SES3 0.7703 0.6365 0.9322 0.0073SES4 0.6451 0.5303 0.7847 0.0000SES5 0.5729 0.4664 0.7036 0.0000Baseline_Comorbidity1 1.3696 1.1687 1.6050 0.0001Baseline_Comorbidity2 1.7545 1.4690 2.0954 0.0000Num_Prev_Hosp 1.0593 1.0424 1.0765 0.0000DMDLow:I(time > 4000):(time - 4000) 1.1769 1.0756 1.2877 0.0004B.2 Cox model results using the cohort from Nget al.Ng et al. [1] included individuals that had missing socioeconomic statusinformation, denoted SES9 in Table B.6, and this fitted model does notinclude pre-index hospitalizations. Its exclusion from our primary exposuremodel caused only a negligible change in the estimate of the HR of exposure,so this should not hinder comparison; see Table B.6. With the exception ofthe exposure HR (see discussion in the main text), the HRs were broadlysimilar to those obtained in the models of this study: the impact of sex isless prominent albeit still large; more income yielded a lower hazard; anycomorbidity increased the hazard and a comorbidity of two raised it evenfurther; older individuals were more likely to die.100B.3. Time-to-exposure modelTable B.6: Cox model using the cohort from Ng et al. [1]HR LB UB pDMD 0.8408 0.7469 0.9464 0.0041SEXM 1.6348 1.5224 1.7555 0.0000Index_Age 1.0835 1.0804 1.0866 0.0000SES2 0.8833 0.7948 0.9818 0.0214SES3 0.7786 0.6988 0.8674 0.0000SES4 0.7394 0.6646 0.8226 0.0000SES5 0.6629 0.5933 0.7407 0.0000SES9 1.3680 0.9692 1.9308 0.0747Baseline_Comorbidity1 1.2148 1.1065 1.3336 0.0000Baseline_Comorbidity2 1.8230 1.6336 2.0344 0.0000B.3 Time-to-exposure modelThe AIC for the initial time-to-exposure model (model summary shown inTable B.7) with calendar period stratification was 33,669.45 and the log-likelihood was -16,825.73.Table B.7: Initial time-to-exposure model with time-dependent stratificationby calendar periodHR LB UB pSEXM 1.0784 0.9786 1.1883 0.1275Index_Age 0.9502 0.9466 0.9539 0.0000SES2 0.9293 0.8082 1.0686 0.3034SES3 0.9769 0.8528 1.1191 0.7363SES4 1.0258 0.8984 1.1713 0.7066SES5 1.0275 0.8959 1.1784 0.6984Comorbidity1 0.8230 0.7237 0.9360 0.0030Comorbidity2 0.6988 0.5723 0.8533 0.0004Num_Prev_Hosp 0.9646 0.9418 0.9880 0.0032B.3.1 Functional formThe number of pre-index hospitalizations was modelled using a log(x + 1)transform since a more appropriate functional form was unclear from thespline fit as shown in Figure B.7. All but one percent of individuals hadhospitalizations lower than 11, and while the log(x + 1) transform is im-perfect, it is consistent with the spline fit in the region where the data is101B.3. Time-to-exposure modelconcentrated as well as with shape of pre-index hospitalizations in othermodels.Figure B.7: Previous hospitalizations modelled with a natural spline versusa log(x+ 1) transformation in red.Figure B.8: Age at index date modelled with a natural spline versus a linearfit in red.Age was left untransformed in spite of evidence that a spline would bemore suitable as shown in Figure B.8. When age was less than 60, theeffect appeared quite linear. It was only when age was greater than 60 that102B.3. Time-to-exposure modelthe effect appeared to deviate from a linear relationship with the hazard.However, only 21 individuals older than 60 at the index date were everexposed. Moreover, as shown below, age has a time-dependent effect. Foreasy interpretation, a linear functional form was used to incorporate a time-dependent effect.B.3.2 Model diagnosticsThe proportional hazards assumption was investigated using scaled Schoen-feld residuals. The proportional hazards assumption appeared violated forthe age variable but Figure B.9 suggested incorporating an interaction withlog(t). log(t) tends to be suitable when most of the events occur early onthe time scale, and the relatively linear form of the natural spline in FigureB.9 indicates that log(t) captures the time-varying effect well.Figure B.9: Scaled Schoenfeld residuals with time on a log scale for thecoefficient of age.B.3.3 Time-to-exposure final model resultsThe AIC of the final time-to-exposure model was 33,638.43 and the log-partial likelihood was -16,809.21. The model summary may be found inTable B.8.103B.3. Time-to-exposure modelTable B.8: Time-to-exposure model with stratification by calendar periodand time-dependent coefficient for ageHR LB UB pSEXM 1.0764 0.9767 1.1864 0.1378Index_Age 0.9583 0.9537 0.9629 0.0000SES2 0.9330 0.8114 1.0729 0.3305SES3 0.9765 0.8524 1.1186 0.7313SES4 1.0319 0.9037 1.1782 0.6426SES5 1.0326 0.9003 1.1843 0.6466Comorbidity1 0.8250 0.7253 0.9384 0.0034Comorbidity2 0.7044 0.5768 0.8602 0.0006log(Num_Prev_Hosp + 1) 0.8792 0.8194 0.9434 0.0003Index_Age:log(time) 0.9870 0.9825 0.9916 0.0000104Appendix CAppendix for Chapter 3C.1 Transition OHC.1.1 Model with primary measure of exposureCandidate modelsThe unstratified model has the highest AIC (435,519.3) followed by themodel stratified by index period (375,045.5) and lastly, the model stratifiedby calendar year of the follow-up segment (372,355.4), respectively summa-rized in Tables C.1, C.2, and C.3. All three produced similar estimatedHRs. The log partial likelihood was -217,747.70 for the unstratified model,-187,510.7 for the model stratified by index period, and -186,165.70 for themodel stratified by calendar period.Table C.1: Unstratified modelHR LB UB pDMD 0.9308 0.8607 1.0067 0.0728SEXM 1.0310 0.9677 1.0984 0.3453Index_Age 1.0209 1.0185 1.0233 0.0000SES2 0.9828 0.8942 1.0802 0.7189SES3 0.9275 0.8517 1.0101 0.0836SES4 0.8955 0.8203 0.9776 0.0136SES5 0.8351 0.7634 0.9136 0.0001Comorbidity1 1.4487 1.3686 1.5335 0.0000Comorbidity2 2.1963 2.0431 2.3611 0.0000Num_Prev_Hosp 1.0787 1.0652 1.0924 0.0000DMD:Comorbidity1 0.9857 0.8563 1.1345 0.8405DMD:Comorbidity2 0.9710 0.8088 1.1656 0.7521105C.1. Transition OHTable C.2: Model with fixed stratification by index periodHR LB UB pDMD 0.9291 0.8587 1.0052 0.0670SEXM 1.0277 0.9647 1.0948 0.3968Index_Age 1.0209 1.0184 1.0233 0.0000SES2 0.9847 0.8963 1.0817 0.7477SES3 0.9268 0.8512 1.0091 0.0799SES4 0.8966 0.8215 0.9786 0.0145SES5 0.8367 0.7655 0.9146 0.0001Comorbidity1 1.4419 1.3623 1.5261 0.0000Comorbidity2 2.1771 2.0250 2.3406 0.0000Num_Prev_Hosp 1.0791 1.0656 1.0926 0.0000DMD:Comorbidity1 0.9847 0.8554 1.1335 0.8298DMD:Comorbidity2 0.9665 0.8051 1.1602 0.7145Table C.3: Model with time-dependent stratification by calendar periodHR LB UB pDMD 0.9291 0.8588 1.0050 0.0665SEXM 1.0286 0.9659 1.0954 0.3790Index_Age 1.0208 1.0184 1.0233 0.0000SES2 0.9844 0.8959 1.0817 0.7435SES3 0.9276 0.8519 1.0100 0.0835SES4 0.8974 0.8223 0.9793 0.0151SES5 0.8372 0.7657 0.9155 0.0001Comorbidity1 1.4417 1.3621 1.5260 0.0000Comorbidity2 2.1722 2.0208 2.3351 0.0000Num_Prev_Hosp 1.0793 1.0660 1.0929 0.0000DMD:Comorbidity1 0.9868 0.8579 1.1350 0.8524DMD:Comorbidity2 0.9670 0.8061 1.1600 0.7178Functional formPre-index hospitalizations did not have a linear relationship with the hazardfunction as shown in Figure C.1 and it had a time-dependent effect as shownin Figure C.4. Two solutions were considered to handle the non-linearity—the spline fit or a log(x+ 1) transformation. Pre-index hospitalizations wasmodelled using the latter instead of the former for three reasons: simulta-neously modelling a time-dependent continuous coefficient and a spline for106C.1. Transition OHa variable is convoluted, and it is unclear how one should best model thatand then interpret it; it captures much of the non-linearity while being moreinterpretable than the spline; it is consistent with the functional forms sug-gested for the other transitions. As demonstrated in the other sections ofthis appendix, pre-index hospitalizations does not have a linear relationshipto the hazard and may be modelled as a log(x + 1) transformation in theO → D, H → O, and H → D Cox models.Figure C.1: Previous hospitalizations modelled with a natural spline (5 d.f.)versus log(x+ 1) in red.Figure C.2: Age at index date modelled with a natural spline (3 d.f.) versusa linear fit in red.107C.1. Transition OHAge was modelled using a natural spline with 3 d.f. and knots placed at the25th, 50th, and 75th percentiles as shown in Figure C.2. Splines with moreknots appeared to overfit the data.Proportional hazards assumptionThe proportional hazards assumption was investigated using Schoenfeldresiduals calculated using the model with time-dependent stratification bycalendar period, age modelled as a spline with 3 d.f., and pre-index hospital-izations modelled using the log(x+ 1) transform. As shown in Figures C.3,C.4, and C.5, the coefficients for exposure and pre-index hospitalizationsappeared to change linearly in time whereas the coefficients for comorbiditychanged during the first 2000 days and then ceased to change. As a result,the model was stratified on comorbidity and an interaction between timeand exposure as well as time and pre-index hospitalizations were specifiedin the model.Figure C.3: Scaled Schoenfeld residuals for the coefficient of exposure.108C.1. Transition OHFigure C.4: Scaled Schoenfeld residuals for the coefficient of pre-index hos-pitalizations.(a) Comorbidity1 (b) Comorbidity2Figure C.5: Scaled Schoenfeld residuals for the coefficients of comorbidity.Model resultsThe final model was specified in Section 3.5.1 and the R output of the coxphobject may be found in Table C.4. The AIC of the final O → H model is320,517.1 and the log-likelihood is -160,244.50. Interpretation of the modelresults may be found after the diagnostics.109C.1. Transition OHTable C.4: Final model for the transition from state O to state HHR LB UB pDMD 0.7480 0.6606 0.8468 0.0000SEXM 1.0600 0.9996 1.1240 0.0517ns(Index_Age, df = 3)1 3.0334 2.5988 3.5408 0.0000ns(Index_Age, df = 3)2 2.7420 1.7240 4.3609 0.0000ns(Index_Age, df = 3)3 1.6345 1.1684 2.2864 0.0041SES2 1.0112 0.9272 1.1027 0.8018SES3 0.9368 0.8626 1.0174 0.1209SES4 0.9206 0.8472 1.0004 0.0511SES5 0.8651 0.7966 0.9396 0.0006log(Num_Prev_Hosp + 1) 2.0006 1.8855 2.1229 0.0000DMD:time 1.0283 1.0165 1.0402 0.0000log(Num_Prev_Hosp + 1):time 0.9760 0.9695 0.9825 0.0000DMD:strata(Comorbidity)1 0.9713 0.8502 1.1096 0.6676DMD:strata(Comorbidity)2 0.8629 0.7160 1.0399 0.1214Model diagnosticsScaled dfbeta residuals indicated that there were 20 potentially influentialindividuals. The model was re-estimated with a dataset that excluded theseindividuals. Most estimates changed very little. However, the interactionterm between exposure and a comorbidity of two changed [HR = 0.80, CI =(0.68, 0.94), p-value = 0.008]. It decreased by about 0.06 but the CI becamenarrower, and the Wald test now showed that there was an effect.Of the 20 individuals that were excluded, 18 were hospitalized with a co-morbidity of two at least once during follow-up. In total 1,730 individualswere hospitalized at least once while they had a comorbidity of two. Theaverage age among these individuals was about 52 and the average time athospitalization was 2,832 days post-index. On average, the 1,730 individualswere hospitalized with a comorbidity of two 3.2 times.The 20 influential individuals were about 15 years younger (average age atthe index date was 37) and had many more hospitalizations (on average,these individuals were hospitalized 14.8 times while having a comorbidity oftwo). These hospitalizations tended to occur later as well (average time-to-hospitalization was 4,018 days post-index).110C.1. Transition OHExcluding these individuals and claiming that there is a favourable associ-ation between individuals with a comorbidity of two and reduced mortalitywould be too data-dependent to be a reasonable conclusion, especially sincethere is no effect among individuals with a comorbidity of one. One wouldexpect some benefit there if the benefit among individuals with more/worsecomorbidities is real. Nevertheless, these individuals are noted here.Figure C.6: Deviance residuals for final O → H modelThe plot of deviance residuals is in Figure C.6. The x-axis is labelled ‘Ob-servation’ which refers to one sojourn. Twelve sojourns had deviances lessthan -4. Most of these sojourns ended with an individual being censored atthe end of their follow-up. Individuals with deviances less than -4 were older(the mean age of these twelve individuals at the index date was 62.8), hadhigher pre-index hospitalizations than normal (average of 9.3), and werenot exposed to DMDs. Sixteen sojourns had deviances greater than 3.5.These sixteen individuals were young (they had an average age at indexdate of 38.8) and fewer pre-index hospitalizations (average of 0.4). Extremedeviances seem to be determined by age at index date and pre-index hospi-talizations. None of the deviance residuals indicated serious problems withthe model.Model interpretationThe number of pre-index hospitalizations was modelled using a log(x + 1)transformation. It had a time-dependent effect as well which was imple-111C.1. Transition OHmented through an interaction between pre-index hospitalizations and timet. At t = 0, the HR of hospitalization was 2.00. Every year afterwards wasassociated with a decrease in the hazard of 2.40%. After ten years, the HRdecreased to 1.57. Pre-index hospitalizations was included in the model asa proxy for the general health of an individual. At the beginning of follow-up, this may work quite well. Over the course of follow-up, an individual’scondition likely changed, and the number of pre-index hospitalizations maynot describe their condition as accurately years or even decades post-index.This loss of relevancy may be reflected in the decreasing effect on the hazard.Individuals in the top quintile of socioeconomic status experienced fewerhospitalizations corresponding to a 13% decrease in the hazard of hospital-ization compared to individuals in the lowest quintile. Men did not appearto have a meaningful change in hazard of hospitalization but older individu-als had a higher hazard than younger individuals. Age was modelled using aspline. The effect of age appears roughly linear until age 70, at which point,the effect plateaus. Individuals older than 70 had a HR of approximately1.75 compared to those at the median index age of 44.C.1.2 Model with secondary measure of exposureProportional hazards assumptionFigure C.7 shows the effect of low DMD exposure, which seems to increaseuntil about 4000 days post-index and remain constant afterwards. Con-versely, Figure C.8 shows that the effect of high exposure before 4000 dayspost-index was roughly constant and then increased afterwards. An interac-tion with time before 4000 days and low exposure as well as an interactionwith time after 4000 days and high exposure was added to the model.The time-dependent effect of sex as shown in Figure C.9 appears similar tothat of high exposure. Thus, an interaction with time after 4000 days andsex was added to the model. Figure C.10 shows the time-dependent effectof the coefficient for one of the basis functions of age at the index date (thecoefficients of the other two basis functions showed a similar shape) whichseems to rise until t = 2000 and then plateau afterwards. An indicator wasused to create two splines in the final model: one for the effect of age up toand including t = 2000 and one afterwards.Lastly, Figure C.11 shows a dip in the estimated coefficient of SES5. Given112C.1. Transition OHthe complexity of the model, that this time-dependent effect can not bemodelled with a simple interaction with time, and that the difference isrelatively small, this suggested time-dependent effect was ignored in thefinal model.Figure C.7: Scaled Schoenfeld residuals for the coefficient of low exposure.Figure C.8: Scaled Schoenfeld residuals for the coefficient of high exposure.113C.1. Transition OHFigure C.9: Scaled Schoenfeld residuals for the coefficient of sex.Figure C.10: Scaled Schoenfeld residuals for the coefficient of the first basisfunction of the spline for age at the index date.114C.1. Transition OHFigure C.11: Scaled Schoenfeld residuals for the coefficient of SES5.Model results and interpretationThe final model has an AIC of 319,785.9 and a log partial likelihood of-159,871.0, and the model summary is shown in Table C.5.The effect of sex did not appear to vary significantly over time despite thesuggestion by the Schoenfeld residuals and overall, sex did not seem tochange the hazard of hospitalization. Pre-index hospitalizations were asso-ciated with a much greater hazard of hospitalization, although this effectattenuated over time. The top two quintiles of socioeconomic status wereassociated with fewer hospitalizations.115C.2. Transition ODTable C.5: Final model using the secondary exposure measure for the tran-sition from state O to state HHR LB UB pDMDLow 0.9127 0.8085 1.0302 0.1392DMDHigh 0.7903 0.6923 0.9023 0.0005SEXM 1.0404 0.9798 1.1047 0.1962log(Num_Prev_Hosp + 1) 2.0048 1.8888 2.1280 0.0000SES2 1.0109 0.9269 1.1025 0.8070SES3 0.9363 0.8621 1.0168 0.1179SES4 0.9199 0.8465 0.9997 0.0491SES5 0.8643 0.7958 0.9388 0.0005ns(Index_Age, df = 3)1:I(time < 2000) 3.2049 2.6961 3.8097 0.0000ns(Index_Age, df = 3)2:I(time < 2000) 2.0882 1.2524 3.4818 0.0048ns(Index_Age, df = 3)3:I(time < 2000) 1.4248 0.9723 2.0880 0.0694ns(Index_Age, df = 3)1:I(time >= 2000) 2.8368 2.2846 3.5225 0.0000ns(Index_Age, df = 3)2:I(time >= 2000) 3.5692 1.8290 6.9650 0.0002ns(Index_Age, df = 3)3:I(time >= 2000) 2.0768 1.2110 3.5618 0.0079DMDLow:strata(Comorbidity)1 0.9975 0.8515 1.1686 0.9758DMDLow:strata(Comorbidity)2 0.9655 0.7736 1.2049 0.7559DMDHigh:strata(Comorbidity)1 0.9492 0.7689 1.1718 0.6278DMDHigh:strata(Comorbidity)2 0.7429 0.5656 0.9757 0.0326DMDLow:I(time <= 4000):time 1.0141 0.9954 1.0331 0.1402DMDHigh:I(time > 4000):(time - 4000) 1.0589 1.0280 1.0906 0.0001SEXM:I(time > 4000):(time - 4000) 1.0202 0.9952 1.0458 0.1136log(Num_Prev_Hosp + 1):time 0.9758 0.9693 0.9824 0.0000C.2 Transition ODC.2.1 Model with primary measure of exposureCandidate modelsThe unstratified model has the highest AIC (8,699.8) followed by the modelstratified by index period (7,363.1) and lastly, the model stratified by cal-endar year of the follow-up segment (7,307.9). All three produced similarestimated HRs. The log partial likelihood was -4,339.90 for the unstratifiedmodel, -3,671.55 for the model stratified by index period, and -3,643.96 forthe model stratified by calendar period.116C.2. Transition ODTable C.6: Unstratified modelHR LB UB pDMD 1.0877 0.7936 1.4907 0.6013SEXM 1.8121 1.5359 2.1379 0.0000Index_Age 1.0944 1.0873 1.1015 0.0000SES2 0.7781 0.6116 0.9899 0.0411SES3 0.7062 0.5535 0.9011 0.0052SES4 0.6439 0.5035 0.8234 0.0004SES5 0.5173 0.3956 0.6764 0.0000Comorbidity1 1.9419 1.5717 2.3993 0.0000Comorbidity2 3.2177 2.6241 3.9456 0.0000Num_Prev_Hosp 1.0411 1.0186 1.0640 0.0003Table C.7: Model with fixed stratification by index periodHR LB UB pDMD 1.0799 0.7876 1.4805 0.6332SEXM 1.8040 1.5287 2.1287 0.0000Index_Age 1.0944 1.0873 1.1016 0.0000SES2 0.7777 0.6109 0.9900 0.0412SES3 0.7090 0.5556 0.9048 0.0057SES4 0.6451 0.5044 0.8252 0.0005SES5 0.5167 0.3950 0.6758 0.0000Comorbidity1 1.9287 1.5600 2.3846 0.0000Comorbidity2 3.2118 2.6155 3.9441 0.0000Num_Prev_Hosp 1.0412 1.0186 1.0643 0.0003117C.2. Transition ODTable C.8: Model with time-dependent stratification by calendar periodHR LB UB pDMD 1.0791 0.7871 1.4794 0.6364SEXM 1.7923 1.5188 2.1151 0.0000Index_Age 1.0947 1.0875 1.1018 0.0000SES2 0.7793 0.6123 0.9920 0.0428SES3 0.7111 0.5572 0.9075 0.0061SES4 0.6445 0.5038 0.8244 0.0005SES5 0.5195 0.3971 0.6797 0.0000Comorbidity1 1.9160 1.5498 2.3687 0.0000Comorbidity2 3.2168 2.6199 3.9496 0.0000Num_Prev_Hosp 1.0420 1.0197 1.0648 0.0002Functional formPre-index hospitalizations was modelled using a log(x + 1) transformationand age at the index date was kept with a linear functional form. Thesefunctional forms are compared to splines in Figure C.12 and Figure C.13.Figure C.12: Previous hospitalizations modelled with a natural spline (5d.f.) versus log(x+ 1) in red.118C.2. Transition ODFigure C.13: Age at index date modelled with a natural spline (5 d.f.) versusa linear fit in red.Proportional hazards assumptionThere appeared to be a violation in the proportional hazards assumptionfor exposure and pre-index hospitalizations as shown in Figures C.14 andC.15. Two interactions with time were added to the model: an interactionbetween time and exposure beginning at t = 4000 and an interaction withtime and pre-index hospitalizations.Figure C.14: Scaled Schoenfeld residuals for the coefficient of exposure.119C.2. Transition ODFigure C.15: Scaled Schoenfeld residuals for the coefficient of pre-indexhospitalizations.Model resultsThe final model has an AIC of 7287.78 and a log partial likelihood of-3,631.89, and the model summary is shown in Table C.9.Table C.9: Final model for the transition from state O to state DHR LB UB pDMD 0.7085 0.4306 1.1659 0.1751SEXM 1.8040 1.5285 2.1293 0.0000Index_Age 1.0917 1.0845 1.0990 0.0000SES2 0.7877 0.6188 1.0027 0.0526SES3 0.7233 0.5667 0.9232 0.0093SES4 0.6581 0.5144 0.8420 0.0009SES5 0.5284 0.4040 0.6911 0.0000Comorbidity1 1.8805 1.5210 2.3251 0.0000Comorbidity2 3.0998 2.5243 3.8064 0.0000log(Num_Prev_Hosp + 1) 1.6594 1.3654 2.0167 0.0000DMD:I(time > 4000):(time - 4000) 1.0488 1.0077 1.0915 0.0193log(Num_Prev_Hosp + 1):time 0.9689 0.9485 0.9897 0.0036120C.2. Transition ODModel diagnosticsThere did not appear to be any influential individuals in the dataset of thismodel. The largest sum of scaled dfbeta residuals was 0.29.There were some deviance residuals (shown in Figure C.16) that required fur-ther investigation: 30 individuals had deviance residuals greater than 3 andfewer than six individuals (exact count not given due to privacy agreementwith Population Data BC) had deviance residuals less than -2. The latterindividuals did not experience the event. They were old at the index daterelative to other individuals, had comorbidities, and higher-than-averagepre-index hospitalizations but did not die during their follow-up (which wasover 15 years).Figure C.16: Deviance residuals for the final O → D modelThe individuals with high deviance residuals are slightly more complex. Theaverage age of an individual that died outside the hospital was over 61 but16 (out of the 30 individuals with a deviance residual greater than 3) were45 years old or younger (the mean age at the index date). The remaining 14individuals (who were older than 45 years old) died quite quickly—in the firstyear and a half of follow-up—but did not have ‘high risk’ characteristics (e.g.,if they were older then they had no comorbidities and high SES or otherwisewere young). Overall, the model appears to explain time-to-death relativelywell. Deviations from the model are inevitable since all-cause mortality willinclude random and unexpected deaths especially in a dataset this large.121C.2. Transition ODModel interpretationMen had a HR of 1.80 compared to women. This striking increase in thehazard of death is consistent with the results presented in Chapter 2. Highersocioeconomic status was associated with a lower hazard. A comorbidity of1 was associated with an 88% increase in the hazard of death while a comor-bidity of 2 was associated with a HR of 3.10 compared to no comorbidities.More pre-index hospitalizations was also associated with a greater chanceof dying although this effect was tempered as follow-up went on.C.2.2 Model with secondary measure of exposureLow exposure to any DMD showed a roughly linear increase in the hazardof death after 4000 days post-index (Figure C.17). As in the model with aprimary measure of exposure, an interaction between time after t = 4000and DMDLow was added to the model. High exposure to any DMD didnot show a clear time-dependent effect (Figure C.18). The effect of pre-index hospitalizations appeared to decrease gradually as the study went on(Figure C.19) so an interaction between pre-index hospitalizations and timewas added for the entire duration of follow-up.Figure C.17: Scaled Schoenfeld residuals for the coefficient of low exposure.122C.2. Transition ODFigure C.18: Scaled Schoenfeld residuals for the coefficient of high exposure.Figure C.19: Scaled Schoenfeld residuals for the coefficient of pre-indexhospitalizations.Model results and interpretationThe final model has an AIC of 7,282.9 and a log partial likelihood of -3,628.5,and the model summary is shown in Table C.10.Being male, older age, having more pre-index hospitalizations, higher co-123C.3. Transitions HOmorbidity, and lower SES were associated with a greater hazard of death.The effect of pre-index hospitalizations decreased over time, as in the modelsfor the O → H transition.Table C.10: Final model using the secondary exposure measure for the tran-sition from state O to state DHR LB UB pDMDLow 0.6954 0.4156 1.1637 0.1667DMDHigh 0.9858 0.5956 1.6319 0.9558SEXM 1.8019 1.5322 2.1190 0.0000Index_Age 1.0917 1.0841 1.0993 0.0000log(Num_Prev_Hosp + 1) 1.6614 1.3680 2.0178 0.0000SES2 0.7888 0.6237 0.9977 0.0478SES3 0.7202 0.5686 0.9122 0.0065SES4 0.6542 0.5145 0.8317 0.0005SES5 0.5311 0.4059 0.6948 0.0000Comorbidity1 1.8833 1.5208 2.3322 0.0000Comorbidity2 3.1014 2.5066 3.8374 0.0000DMDLow:I(time > 4000):(time - 4000) 1.2465 1.1063 1.4044 0.0003log(Num_Prev_Hosp + 1):time 0.9686 0.9481 0.9895 0.0034C.3 Transitions HOC.3.1 Model with primary measure of exposureCandidate modelsThe unstratified model has the highest AIC (484,242.5) followed by themodel stratified by index period (411,795.9) and lastly, the model stratifiedby calendar period of the follow-up segment (411,028.1). The log-likelihoodof the unstratified model, model stratified by index period, and model strat-ified by calendar period was -242,109.2, -205,8885.9, -205,502.0 respectively.124C.3. Transitions HOTable C.11: Unstratified modelHR LB UB pDMD 0.9673 0.9232 1.0136 0.1631SEXM 0.9132 0.8802 0.9474 0.0000Admission_Age 0.9935 0.9923 0.9947 0.0000SES2 1.0801 1.0274 1.1356 0.0025SES3 1.1146 1.0605 1.1715 0.0000SES4 1.1131 1.0575 1.1717 0.0000SES5 1.1503 1.0930 1.2105 0.0000Comorbidity1 0.9479 0.9143 0.9826 0.0036Comorbidity2 0.8966 0.8630 0.9315 0.0000Num_Prev_Hosp 0.9876 0.9834 0.9918 0.0000DMD:Comorbidity1 1.1070 1.0195 1.2020 0.0155DMD:Comorbidity2 1.0285 0.9126 1.1590 0.6451Table C.12: Model with fixed stratification by index periodHR LB UB pDMD 0.9684 0.9240 1.0149 0.1798SEXM 0.9124 0.8797 0.9463 0.0000Admission_Age 0.9936 0.9924 0.9948 0.0000SES2 1.0799 1.0276 1.1348 0.0024SES3 1.1132 1.0595 1.1696 0.0000SES4 1.1132 1.0582 1.1711 0.0000SES5 1.1515 1.0945 1.2115 0.0000Comorbidity1 0.9468 0.9135 0.9812 0.0027Comorbidity2 0.8944 0.8608 0.9294 0.0000Num_Prev_Hosp 0.9875 0.9832 0.9917 0.0000DMD:Comorbidity1 1.1097 1.0220 1.2050 0.0132DMD:Comorbidity2 1.0323 0.9166 1.1628 0.5999125C.3. Transitions HOTable C.13: Model with time-dependent stratification by calendar periodHR LB UB pDMD 0.9489 0.9053 0.9945 0.0287SEXM 0.9113 0.8782 0.9457 0.0000Admission_Age 0.9930 0.9918 0.9943 0.0000SES2 1.0782 1.0255 1.1336 0.0032SES3 1.1134 1.0594 1.1702 0.0000SES4 1.1131 1.0574 1.1718 0.0000SES5 1.1485 1.0908 1.2093 0.0000Comorbidity1 0.9431 0.9099 0.9775 0.0014Comorbidity2 0.8877 0.8546 0.9222 0.0000Num_Prev_Hosp 0.9881 0.9838 0.9924 0.0000DMD:Comorbidity1 1.1113 1.0242 1.2057 0.0113DMD:Comorbidity2 1.0358 0.9212 1.1647 0.5566Functional formThe functional forms of pre-index hospitalizations and age at admission tothe hospital are explored in Figures C.20 and C.21 for the model with onlyone set of coefficients. For the final model, pre-index hospitalizations wasmodelled using a log(x + 1) transformation and age at admission to thehospital was modelled with a 3 d.f. spline. With fewer internal knots, thespline did not adequately capture the relationship, while with more knots,the spline appeared to overfit the data.126C.3. Transitions HOFigure C.20: Previous hospitalizations modelled with a natural spline (5d.f.) versus log(x+ 1) in red.Figure C.21: Age at admission date modelled with a natural spline (3 d.f.)versus a linear fit in red.Proportional hazards assumptionTo accomodate potentially different effects for very short (< 1 day) andvery long (> 90 days) sojourns, the final model was fit with three separatesets of coefficients for the three types of sojourns. Time-dependent effects127C.3. Transitions HOwere only investigated for the second set, when 1 ≤ t∗ ≤ 90. One day wastoo short an interval for a time-dependent effect to exist. Time-dependentcovariate effects after 90 days were ignored for two reasons. First, therewere comparably fewer events (223) over quite a long time after 90 days(longest discharge time was 1,998 days after admission). It is unlikely thattime-dependent effects could be adequately modelled given the relative in-frequency of discharge at many times after 90 days. Second, many of thesehospitalizations were likely long-term care and are not the focus of the study.Because there are three sets of coefficients, any possible time-dependenceafter 90 days that is left unmodelled does not interfere with estimates ofexposure for shorter hospitalizations. The time spent in hospital after 90days was not excluded from the model so that individuals who experiencedthese hospitalizations remained under observation and in a state at all timesduring their follow-up.A Cox model was fit using only observations from 1 ≤ t∗ ≤ 90 and the scaledSchoenfeld residuals were calculated from that model. The plot of scaledSchoenfeld residuals for exposure are shown in Figure C.22. The spline doesnot indicate a time-dependent effect since one can draw a horizontal lineand stay within the bounds of the pointwise confidence intervals.Figure C.22: Scaled Schoenfeld residuals for the coefficient of exposure.The following covariates appeared to have a time-dependent effect: co-morbidity (Figure C.23), pre-index hospitalizations (Figure C.24), socioe-conomic status (but only for the highest quintile as shown in Figure C.25),128C.3. Transitions HOsex (Figure C.26), and age (Figure C.27).(a) Comorbidity1 (b) Comorbidity2Figure C.23: Scaled Schoenfeld residuals for the coefficients of comorbidity.Figure C.24: Scaled Schoenfeld residuals for the coefficient of pre-indexhospitalizations.129C.3. Transitions HOFigure C.25: Scaled Schoenfeld residuals for the coefficient of SES5.Figure C.26: Scaled Schoenfeld residuals for the coefficient of sex.130C.3. Transitions HOFigure C.27: Scaled Schoenfeld residuals for the coefficient of the third basisfunction for the spline modelling age at admission.The model was stratified on comorbidity. An interaction between sojourntime and sex was added but was not significant so it was excluded in the finalmodel. An interaction with log sojourn time was added between pre-indexhospitalizations and the highest quintile of socioeconomic status. It was lessobvious how to handle the time-dependent effect for age at admission, since itwas modelled with a spline and the scaled Schoenfeld residuals indicated thatonly the third basis function had a time-dependent effect. The effect appearsto change fairly linearly but to avoid adding even more complexity to themodel, an indicator function was used to fit separate splines for 1 ≤ t∗ < 10and 10 ≤ t∗ ≤ 90 instead of a continuous time-dependent effect. Thus, therewere four different splines to estimate the effect of age at admission on thehazard.Model resultsThe final model has an AIC of 361,198.80 and a log-likelihood of -180,558.40.Summaries of all three sets of coefficients for the HO model may be foundin Tables C.14, C.15, and C.16.131C.3. Transitions HOTable C.14: First set of coefficients where t∗ < 1HR LB UB pDMD 0.9394 0.8920 0.9894 0.0182SEXM 0.9016 0.8582 0.9472 0.0000ns(Admission_Age, df = 3)1 1.1859 1.0690 1.3156 0.0013ns(Admission_Age, df = 3)2 0.7684 0.5314 1.1111 0.1615ns(Admission_Age, df = 3)3 0.3174 0.2411 0.4177 0.0000SES2 1.0824 1.0059 1.1648 0.0342SES3 1.1336 1.0585 1.2141 0.0003SES4 1.1390 1.0637 1.2197 0.0002SES5 1.1431 1.0684 1.2230 0.0001log(Num_Prev_Hosp + 1) 0.9049 0.8763 0.9345 0.0000DMD:strata(Comorbidity)1 1.0645 0.9624 1.1774 0.2242DMD:strata(Comorbidity)2 1.0243 0.8863 1.1837 0.7453Table C.15: Second set of coefficients where 1 ≤ t∗ ≤ 90HR LB UB pDMD 0.9435 0.8764 1.0159 0.1232SEXM 0.9193 0.8729 0.9681 0.0014ns(Admission_Age, df = 3)1:I(1 <= t* < 10) 0.7129 0.6102 0.8330 0.0000ns(Admission_Age, df = 3)2:I(1 <= t* < 10) 0.8819 0.4781 1.6266 0.6874ns(Admission_Age, df = 3)3:I(1 <= t* < 10) 0.4869 0.3667 0.6465 0.0000ns(Admission_Age, df = 3)1:I(10 <= t* <= 90) 0.6650 0.5144 0.8598 0.0019ns(Admission_Age, df = 3)2:I(10 <= t* <= 90) 0.4611 0.1703 1.2481 0.1276ns(Admission_Age, df = 3)3:I(10 <= t* <= 90) 0.7185 0.5202 0.9925 0.0449SES2 1.0781 1.0023 1.1597 0.0432SES3 1.0872 1.0098 1.1705 0.0265SES4 1.0783 1.0026 1.1598 0.0425SES5 1.2265 1.1092 1.3563 0.0001SES5:log(t*) 0.9395 0.8987 0.9821 0.0058log(Num_Prev_Hosp + 1) 0.9267 0.8819 0.9737 0.0026log(Num_Prev_Hosp + 1):log(t*) 1.0419 1.0164 1.0680 0.0012DMD:strata(Comorbidity)1 1.1588 1.0070 1.3335 0.0396DMD:strata(Comorbidity)2 1.0004 0.8200 1.2203 0.9972132C.3. Transitions HOTable C.16: Third set of coefficients where t∗ > 90HR LB UB pDMD 1.2574 0.7931 1.9934 0.3300SEXM 0.6918 0.5155 0.9285 0.0141ns(Admission_Age, df = 3)1 0.9215 0.4976 1.7063 0.7947ns(Admission_Age, df = 3)2 1.1737 0.1853 7.4336 0.8649ns(Admission_Age, df = 3)3 2.2636 0.7081 7.2360 0.1682SES2 0.8942 0.5680 1.4076 0.6290SES3 1.1801 0.7591 1.8346 0.4619SES4 0.9591 0.6217 1.4796 0.8502SES5 0.9362 0.5591 1.5676 0.8021log(Num_Prev_Hosp + 1) 0.8989 0.7535 1.0724 0.2365DMD:strata(Comorbidity)1 1.0183 0.2362 4.3895 0.9806DMD:strata(Comorbidity)2 0.6860 0.2509 1.8758 0.4628Model diagnosticsEight individuals had scaled dfbeta residuals greater than one-third. Ex-cluding them from the dataset and refitting the model changed most of thecoefficients by less than 0.02. This suggests the individuals are not influen-tial.As can be seen in Figure C.28, most deviance residuals were clustered be-tween 0.4 and 0.8; few were above one. The largest positive deviance residualwas 1.4 but sixteen were below -3 and the largest negative deviance resid-ual was -3.7. These negative residuals were hospitalizations longer than 90days; the mean length was 232 days, and they were likely a different typeof hospital visit than this study aims to explore. They were handled in themodel through their own set of coefficients.133C.3. Transitions HOFigure C.28: Deviance residuals for the final H → O modelModel interpretationOn the first day, male patients had a 10% decrease in the hazard of beingdischarged. Higher socioeconomic status (3rd, 4th, and 5th quintiles) wasassociated with a HR of 1.13 - 1.14 compared to the lowest income quintile.The effect was less prominent for the second lowest quintile. Higher pre-index hospitalizations was also associated with a lower chance of leaving thehospital.The effect of sex was roughly the same in the second interval of sojourn time(between one day and 90 days). The effect of SES2, SES3, and SES4 weresimilar: about an 8% greater hazard of being discharged. Compared to thelowest income quintile, SES5 had a greater effect. Being in the top incomequintile was associated with a 23% greater hazard of being discharged att∗ = 1. As the number of days spent in the hospital increased, the effect ofbeing in the highest SES decreased.For hospitalizations lasting longer than 90 days, only sex was associatedwith discharge time. Men were still less likely to be discharged and thiseffect was stronger than in the two prior sets of coefficients (HR of 0.69).134C.3. Transitions HOC.3.2 Model with secondary measure of exposureAs in the model with the primary measure of exposure, time-dependentcovariate effects were checked only for overnight hospitalizations that wereno longer than 90 days. The results were similar to the model with theprimary measure of exposure for sex, the highest socioeconomic status, pre-index hospitalizations, and comorbidity. Those plots are not shown to avoidredundancy. The Schoenfeld residual plots for the exposure coefficients areshown in Figure C.29 and Figure C.30. The latter shows a very slightdecrease but it still seems that a horizontal line could be drawn withinthe bounds of the CIs. No time-dependent effect was modelled for thosecoefficients. The coefficient for the third basis function of the spline usedto model age at hospital admission again showed a time-dependent effectas shown in Figure C.31. An interaction term was used to fit two separatesplines for the second set of coefficients, one when 1 ≤ t∗ < 10 and anotherwhen 10 ≤ t∗ ≤ 90.Figure C.29: Scaled Schoenfeld residuals for the coefficient of low exposure.135C.3. Transitions HOFigure C.30: Scaled Schoenfeld residuals for the coefficient of high exposureFigure C.31: Scaled Schoenfeld residuals for the coefficient of the third basisfunction for the spline modelling age at admission.Model results and interpretationThe final model has an AIC of 361,183.3 and a log partial likelihood of-180,550.6, and the model summary is shown in Tables C.17, C.18, and C.19.On the first day, male patients and patients with more pre-index hospital-136C.3. Transitions HOizations were discharged less. Relative to the lowest socioeconomic quintile,higher socioeconomic status was associated with a greater chance of dis-charge.These effects appeared to be relatively consistent with those for betweenday one and day 90 in the hospital although the hazard of discharge forindividuals in the highest socioeconomic status did increase relative to theother quintiles.Table C.19 has dashes instead of an estimated coefficient for the interactionterm between high exposure to DMDs and a comorbidity of one becausefewer than six individuals (the exact number is hidden for privacy reasons)had high exposure and a comorbidity of one during a hospitalization lastingat least 90 days. The model failed to return an estimate. After 90 days,compared to women, men were discharged even less than in the first 90days of hospitalization (HR of 0.69). There was no evidence that SES orpre-index hospitalizations were associated with any change in the hazard ofdischarge after 90 days.Table C.17: First set of coefficients where t∗ < 1HR LB UB pDMDLow 0.9144 0.8612 0.9709 0.0034DMDHigh 0.9942 0.9175 1.0773 0.8867SEXM 0.9016 0.8582 0.9471 0.0000ns(Admission_Age, df = 3)1 1.1814 1.0645 1.3111 0.0017ns(Admission_Age, df = 3)2 0.7611 0.5261 1.1012 0.1475ns(Admission_Age, df = 3)3 0.3174 0.2411 0.4177 0.0000SES2 1.0831 1.0065 1.1654 0.0329SES3 1.1341 1.0589 1.2146 0.0003SES4 1.1385 1.0633 1.2191 0.0002SES5 1.1432 1.0685 1.2232 0.0001log(Num_Prev_Hosp + 1) 0.9050 0.8764 0.9345 0.0000DMDLow:strata(Comorbidity)1 1.1167 0.9911 1.2583 0.0699DMDLow:strata(Comorbidity)2 1.0115 0.8393 1.2190 0.9046DMDHigh:strata(Comorbidity)1 0.9571 0.8142 1.1250 0.5947DMDHigh:strata(Comorbidity)2 1.0494 0.8734 1.2609 0.6064137C.3. Transitions HOTable C.18: Second set of coefficients where 1 ≤ t∗ ≤ 90HR LB UB pDMDLow 0.9597 0.8855 1.0401 0.3165DMDHigh 0.9023 0.7840 1.0385 0.1518SEXM 0.9194 0.8730 0.9682 0.0015ns(Admission_Age, df = 3)1:I(1 <= t* < 10) 0.7141 0.6108 0.8348 0.0000ns(Admission_Age, df = 3)2:I(1 <= t* < 10) 0.8853 0.4798 1.6336 0.6968ns(Admission_Age, df = 3)3:I(1 <= t* < 10) 0.4869 0.3666 0.6465 0.0000ns(Admission_Age, df = 3)1:I(10 <= t* <= 90) 0.6653 0.5147 0.8601 0.0019ns(Admission_Age, df = 3)2:I(10 <= t* <= 90) 0.4634 0.1711 1.2551 0.1303ns(Admission_Age, df = 3)3:I(10 <= t* <= 90) 0.7187 0.5203 0.9928 0.0451SES2 1.0779 1.0020 1.1596 0.0440SES3 1.0868 1.0094 1.1701 0.0272SES4 1.0790 1.0033 1.1604 0.0405SES5 1.2262 1.1089 1.3559 0.0001log(Num_Prev_Hosp + 1) 0.9267 0.8819 0.9737 0.0026DMDLow:strata(Comorbidity)1 1.1264 0.9487 1.3373 0.1744DMDLow:strata(Comorbidity)2 1.0062 0.7821 1.2944 0.9617DMDHigh:strata(Comorbidity)1 1.2425 1.0154 1.5205 0.0350DMDHigh:strata(Comorbidity)2 0.9841 0.7487 1.2935 0.9087SES5_1to90:log(t*) 0.9395 0.8988 0.9821 0.0058log(Num_Prev_Hosp + 1):log(t*) 1.0419 1.0164 1.0680 0.0012138C.4. Transition HDTable C.19: Third set of coefficients where t∗ > 90HR LB UB pDMDLow 1.1743 0.6596 2.0908 0.5851DMDHigh 1.4692 0.7209 2.9946 0.2896SEXM 0.6855 0.5093 0.9227 0.0128ns(Admission_Age, df = 3)1 0.9160 0.4805 1.7463 0.7899ns(Admission_Age, df = 3)2 1.1225 0.1732 7.2764 0.9035ns(Admission_Age, df = 3)3 2.2674 0.7030 7.3125 0.1706SES2 0.8689 0.5460 1.3827 0.5532SES3 1.2049 0.7684 1.8893 0.4166SES4 1.0176 0.6617 1.5648 0.9368SES5 0.9495 0.5653 1.5948 0.8447log(Num_Prev_Hosp + 1) 0.8906 0.7454 1.0641 0.2020DMDLow:strata(Comorbidity)1 1.0551 0.2393 4.6527 0.9436DMDLow:strata(Comorbidity)2 1.5455 0.6946 3.4388 0.2861DMDHigh:strata(Comorbidity)1 ------ ------ ------ ------DMDHigh:strata(Comorbidity)2 0.1424 0.0137 1.4843 0.1032C.4 Transition HDC.4.1 Model with primary measure of exposureCandidate modelsThe unstratified model has the highest AIC (5,395.6) and lowest log-likelihood(-2,687.79) followed by the model stratified by index period (AIC = 4,375.7,log-likelihood = -2,177.87) and lastly, the model stratified by calendar yearof the follow-up segment (AIC = 4,328.8, log-likelihood = -2,154.42).139C.4. Transition HDTable C.20: Unstratified modelHR LB UB pDMD 1.1112 0.7746 1.5940 0.5668SEXM 1.6257 1.3225 1.9984 0.0000Admission_Age 1.0285 1.0205 1.0367 0.0000SES2 0.9261 0.6734 1.2734 0.6364SES3 1.0622 0.7813 1.4441 0.7002SES4 0.8552 0.6177 1.1841 0.3462SES5 1.1680 0.8461 1.6124 0.3451Comorbidity1 1.4739 1.1109 1.9555 0.0072Comorbidity2 2.4181 1.8904 3.0930 0.0000Num_Prev_Hosp 0.9571 0.9286 0.9865 0.0045Table C.21: Model with fixed stratification by index periodHR LB UB pDMD 1.0630 0.7386 1.5300 0.7422SEXM 1.6573 1.3426 2.0459 0.0000Admission_Age 1.0273 1.0191 1.0356 0.0000SES2 0.9259 0.6720 1.2757 0.6377SES3 1.0630 0.7804 1.4480 0.6985SES4 0.8420 0.6057 1.1706 0.3064SES5 1.1653 0.8430 1.6109 0.3544Comorbidity1 1.4571 1.0954 1.9382 0.0097Comorbidity2 2.4398 1.9050 3.1246 0.0000Num_Prev_Hosp 0.9590 0.9302 0.9887 0.0071140C.4. Transition HDTable C.22: Model with time-dependent stratification by calendar periodHR LB UB pDMD 1.1384 0.7915 1.6373 0.4847SEXM 1.6207 1.3159 1.9961 0.0000Admission_Age 1.0287 1.0206 1.0369 0.0000SES2 0.9361 0.6795 1.2896 0.6862SES3 1.0748 0.7885 1.4650 0.6480SES4 0.8748 0.6299 1.2151 0.4250SES5 1.1793 0.8526 1.6313 0.3191Comorbidity1 1.5180 1.1417 2.0184 0.0041Comorbidity2 2.5050 1.9503 3.2174 0.0000Num_Prev_Hosp 0.9573 0.9286 0.9870 0.0051Functional formPre-index hospitalizations was modelled using a log(x + 1) transformationand age at admission was kept with a linear functional form. These func-tional forms are compared to splines in Figures C.32 and C.33.Figure C.32: Previous hospitalizations modelled with a natural spline (5d.f.) versus log(x+ 1) in red.141C.4. Transition HDFigure C.33: Age at admission modelled with a natural spline (5 d.f.) versusa linear fit in red.Proportional hazards assumptionThe proportional hazards assumption did not seem to be violated for anyof the covariates. The scaled Schoenfeld residual plot for exposure is never-theless shown in Figure C.34.Figure C.34: Scaled Schoenfeld residuals for the coefficient of exposure.142C.4. Transition HDModel resultsThe final model has an AIC of 4,327.12 and a log-likelihood of -2,153.56.Table C.23 shows a summary of the model.Table C.23: Final model for transition from state H to state DHR LB UB pDMD 1.1284 0.7844 1.6233 0.5150SEXM 1.6239 1.3183 2.0002 0.0000Admission_Age 1.0293 1.0212 1.0376 0.0000SES2 0.9277 0.6734 1.2781 0.6463SES3 1.0697 0.7848 1.4581 0.6696SES4 0.8684 0.6252 1.2061 0.3999SES5 1.1586 0.8380 1.6020 0.3731Comorbidity1 1.5354 1.1543 2.0422 0.0032Comorbidity2 2.5641 1.9912 3.3019 0.0000log(Num_Prev_Hosp + 1) 0.7934 0.6924 0.9090 0.0009Model diagnosticsThere did not appear to be any influential individuals in the dataset of thismodel. The largest scaled dfbeta residuals was 0.18.Deviance residuals are pictured in Figure C.35. Fewer than six individualshad a deviance residual below -2. Similarly to the O → D model, they didnot die, were male, and tended to be older at follow-up. These individualstended to have longer hospitalizations (greater than 90 days). 11 individualshad deviance residuals greater than 3 and the largest residual was 3.72. Ninedied in the first or second day after hospital admission, and these individualstended not to have comorbidities. This model appears to fit the individualsfairly well. Given the semi-Markov time scale, many events occurred on thevery first day which can skew the residuals.143C.4. Transition HDFigure C.35: Deviance residuals for the final H → D modelModel interpretationMale sex increased the hazard by a lower amount than in the O → D tran-sition and age at admission had less of an effect than age at the index datein the O → D transition as well. Some effects differed more dramaticallyas well. Higher socioeconomic status was not associated with reduced mor-tality. Within the hospital, more pre-index hospitalizations was associatedwith a decrease in the hazard of death. This last point is unusual and it isnot clear why this is the case. It does not appear to be a time-dependenteffect since an interaction with t∗, the time since entry to the hospital, wasnot significant.C.4.2 Model with secondary measure of exposureThe proportional hazards assumption seems to hold for every covariate. Thescaled Schoenfeld residuals for the coefficients of high and low exposure areshown in Figures C.36 and C.37.144C.4. Transition HDFigure C.36: Scaled Schoenfeld residuals for the coefficient of low exposureFigure C.37: Scaled Schoenfeld residuals for the coefficient of high exposureModel results and interpretationThe final model has an AIC of 4,329.02 and a log partial likelihood of-2,153.51, and the model summary is shown in Table C.24.The coefficients for sex, age at admission, socioeconomic status, comorbidity,and pre-index hospitalizations are nearly identical as in the model with the145C.4. Transition HDprimary measure of exposure.Table C.24: Final model for transition from state H to state DHR LB UB pDMDLow 1.1691 0.7899 1.7303 0.4347DMDHigh 1.0341 0.5383 1.9865 0.9198SEXM 1.6239 1.3169 2.0024 0.0000Admission_Age 1.0293 1.0212 1.0376 0.0000SES2 0.9275 0.6782 1.2683 0.6373SES3 1.0688 0.7866 1.4522 0.6705SES4 0.8675 0.6311 1.1926 0.3815SES5 1.1575 0.8478 1.5804 0.3572Comorbidity1 1.5350 1.1458 2.0563 0.0041Comorbidity2 2.5624 1.9596 3.3507 0.0000log(Num_Prev_Hosp + 1) 0.7933 0.6976 0.9021 0.0004146Appendix DAppendix for Chapter 4D.1 Pseudo-value Model O for Time Period 1Table D.1: Summary of Pseudo-value Model O for Time Period 1Beta LB UB p(Intercept) 735.7340 730.5594 740.9087 0.0000exposure 0.5394 -1.5925 2.6712 0.6200ns(Index_Age, df = 3)1 9.1423 -1.0055 19.2901 0.0774ns(Index_Age, df = 3)2 -84.2834 -112.4157 -56.1511 0.0000ns(Index_Age, df = 3)3 -139.6290 -182.9331 -96.3249 0.0000SEXM -5.6694 -9.1449 -2.1939 0.0014SES2 -0.1681 -4.4597 4.1235 0.9388SES3 -2.2502 -6.6211 2.1207 0.3130SES4 -1.0580 -5.2181 3.1021 0.6182SES5 2.0767 -1.6165 5.7700 0.2704Comorbidity1 -5.4757 -10.4004 -0.5511 0.0293Comorbidity2 -15.0450 -23.7201 -6.3698 0.0007log(Num_Prev_Hosp + 1) -4.4008 -6.7274 -2.0741 0.0002exposure:Comorbidity1 4.3855 -2.0666 10.8377 0.1828exposure:Comorbidity2 16.8474 7.8413 25.8534 0.0002HC SE: 2.64, 1.09, 5.18, 14.35, 22.09, 1.77, 2.19, 2.23,2.12, 1.88, 2.51, 4.43, 1.19, 3.29, 4.59.Wild boot: 2.64, 1.10, 5.25, 14.37, 22.21, 1.77, 2.21, 2.22,2.08, 1.86, 2.49, 4.51, 1.18, 3.27, 4.68.Table D.1 shows the summary of pseudo-value Model O for Time Period 1.Figure D.1 indicates that once an individual is older than 50, they begin toexperience fewer days alive and out of the hospital between 2 and 4 yearspost-index. The oldest individuals in the dataset (over 90 years old at theindex date) appeared to have about 5 months less time in state O. Being147D.1. Pseudo-value Model O for Time Period 1male and having more pre-index hospitalizations was associated with lesstime spent in state O. Higher comorbidity at 2 years post-index was associ-ated with fewer days spent in state O in the next two years, and the effectwas greater for individuals with a comorbidity of two versus a comorbidityof one. Socioeconomic status did not seem to have an effect on time spentin state O.The plot of residuals versus fitted values for Model O for Time Period 1is shown in Figure D.2. The red line is a natural spline with five d.f. Itgets dragged down at the lower end of the x-axis similar to the pseudo-valueresidual plots in [42]. The histogram of the fitted values, shown in FigureD.3, makes clear that the residuals are highly left-skewed. The smallest 1%of the fitted values were removed and the spline fit again. This is shown inFigure D.4. The line is now nearly perfectly horizontal, which suggests thatthe pattern in Figure D.2 was due to the few fitted values below 600. Thesplines in the residual plots in Figure D.5 are centred at zero so the trans-formations of age and pre-index hospitalizations appear suitable. Overall,the residuals plots do not indicate any serious issue with the model.Figure D.1: Age at index date modelled with a natural spline148D.1. Pseudo-value Model O for Time Period 1Figure D.2: Residual plot for Model O for time period 1Figure D.3: Residual plot for Model O for time period 1, smallest 1% offitted values removed149D.2. Pseudo-value Model O for Time Period 2Figure D.4: Histogram of fitted pseudo-values for Model O for time period1(a) Age at index date (b) Pre-index hospitalizationsFigure D.5: Residual plots for covariates in Model O for Time Period 1D.2 Pseudo-value Model O for Time Period 2Table D.2 shows the summary of pseudo-value Model O for Time Period 2.Age in Model O for Time Period 2 exhibits a similar pattern to Model Ofor Time Period 1—older individuals are in state O for a shorter amountof time. The effect of sex, pre-index hospitalizations, and comorbidity ap-peared qualitatively consistent as well, e.g., the effect of Comorbidity2 is stilllarger than the effect of Comorbidity1 but the magnitude of the difference isnow much greater. SES now appears to have an effect; higher socioeconomic150D.2. Pseudo-value Model O for Time Period 2status (top two quintiles) was generally associated with more time spent outof hospital and alive. The HC standard errors seem appropriate as well.Table D.2: Summary of Pseudo-value Model O for Time Period 2Beta LB UB p(Intercept) 1107.7878 1098.0149 1117.5608 0.0000exposure -0.6775 -4.0586 2.7037 0.6945ns(Index_Age, df = 3)1 29.1023 8.1938 50.0108 0.0064ns(Index_Age, df = 3)2 -206.8690 -265.5956 -148.1423 0.0000ns(Index_Age, df = 3)3 -331.3283 -419.4102 -243.2465 0.0000SEXM -15.2068 -21.8000 -8.6137 0.0000SES2 3.5883 -5.6044 12.7810 0.4442SES3 6.5211 -2.0255 15.0677 0.1348SES4 8.6451 0.1066 17.1836 0.0472SES5 11.9130 3.9057 19.9203 0.0035Comorbidity1 -8.3911 -17.9916 1.2095 0.0867Comorbidity2 -48.6235 -68.1471 -29.0998 0.0000log(Num_Prev_Hosp + 1) -9.8291 -14.7787 -4.8796 0.0001exposure:Comorbidity1 -3.2331 -21.3097 14.8434 0.7259exposure:Comorbidity2 54.2558 33.9771 74.5346 0.0000HC SE: 4.99, 1.73, 10.67, 29.96, 44.94, 3.36, 4.69, 4.36,4.36, 4.09, 4.90, 9.96, 2.53, 9.22, 10.35Wild boot: 4.97, 1.71, 10.64, 29.99, 44.94, 3.37, 4.61, 4.40,4.41, 4.05, 4.88, 9.87, 2.55, 9.24, 10.21151D.2. Pseudo-value Model O for Time Period 2Figure D.6: Age at index date modelled with a natural spline for Model Ofor Time Period 2The residual plots do not indicate an issue with the model. Figure D.7 showsthe spline falling quite far below zero at the lower end of fitted values. This islikely because there are few data points in the lower half of the range. Whenthe smallest 1% of fitted values were removed, as shown in Figure D.8, thespline indicates that the residuals are centred at zero, demonstrating thatit is sensitive to the relatively few fitted values in this range. The residual-covariate plots in Figure D.9 show splines that are centred at zero as well.Figure D.7: Residual plot for pseudo-value Model O for Time Period 2152D.3. Pseudo-value Model O for Time Period 3Figure D.8: Residual plot for pseudo-value Model O for Time Period 2 withsmallest 1% of fitted values removed(a) Age at index date (b) Pre-index hospitalizationsFigure D.9: Residual plots for covariates in Model O for Time Period 2D.3 Pseudo-value Model O for Time Period 3Table D.3 shows the summary of pseudo-value Model O for Time Period3. The spline for age has the same shape as it did in earlier time periodsas shown in Figure D.10. Once again, higher socioeconomic status wasgenerally associated with more time spent in state O relative to the lowestquintile, and any comorbidity was associated with fewer days out of thehospital. Being male and more pre-index hospitalizations were associatedwith fewer days in state O.153D.3. Pseudo-value Model O for Time Period 3Table D.3: Summary of Pseudo-value Model O for Time Period 3Beta LB UB p(Intercept) 1105.8909 1093.6983 1118.0835 0.0000exposure -3.6360 -8.4954 1.2234 0.1425ns(Index_Age, df = 3)1 26.0183 6.3293 45.7073 0.0096ns(Index_Age, df = 3)2 -214.2332 -273.3812 -155.0853 0.0000ns(Index_Age, df = 3)3 -332.9181 -420.5948 -245.2414 0.0000SEXM -9.7286 -17.3060 -2.1511 0.0119SES2 12.7041 2.2746 23.1335 0.0170SES3 9.6038 -0.3550 19.5626 0.0587SES4 6.8895 -3.6491 17.4282 0.2001SES5 10.9405 0.9953 20.8858 0.0311Comorbidity1 -23.0750 -35.2453 -10.9047 0.0002Comorbidity2 -40.2183 -58.5729 -21.8637 0.0000log(Num_Prev_Hosp + 1) -4.7650 -9.7302 0.2003 0.0600exposure:Comorbidity1 22.1514 4.4096 39.8933 0.0144exposure:Comorbidity2 28.2239 -4.2701 60.7179 0.0887HC SE: 6.22 2.48 10.05 30.18 44.73 3.87 5.32 5.08,5.38 5.07 6.21 9.36 2.53 9.05 16.58Wild boot: 6.27 2.49 9.99 30.13 44.97 3.93 5.29 5.00,5.38 5.02 6.18 9.38 2.53 9.07 16.34Figure D.10: Age at index date for Model O for Time Period 3154D.3. Pseudo-value Model O for Time Period 3The residual plots do not indicate a serious issue with the model. There issome deviation from zero in the splines of Figures D.11, D.12, and D.13 butnone that are problematic since they occur near very few fitted values.Figure D.11: Residual plot for pseudo-value Model O for Time Period 3Figure D.12: Residual plot for pseudo-value Model O for Time Period 3with smallest 1% of fitted values removed155D.4. Pseudo-value Model O for Time Period 4(a) Age at index date (b) Pre-index hospitalizationsFigure D.13: Residual plots for covariates in Model O for Time Period 3D.4 Pseudo-value Model O for Time Period 4Table D.4: Summary of Pseudo-value Model O for Time Period 4Beta LB UB p(Intercept) 4357.2399 4214.2143 4500.2656 0.0000exposure -71.0846 -144.1609 1.9917 0.0566ns(Index_Age, df = 3)1 -2.9982 -168.2919 162.2956 0.9716ns(Index_Age, df = 3)2 -2024.2279 -2441.8738 -1606.5820 0.0000ns(Index_Age, df = 3)3 -3092.6717 -3630.1160 -2555.2273 0.0000SEXM -182.7556 -264.6739 -100.8374 0.0000SES2 128.6718 13.6792 243.6644 0.0283SES3 162.5939 53.6001 271.5877 0.0035SES4 178.8661 70.6770 287.0553 0.0012SES5 180.7003 70.4673 290.9333 0.0013Comorbidity1 -49.6819 -163.1646 63.8008 0.3909Comorbidity2 -210.5318 -377.0431 -44.0206 0.0132log(Num_Prev_Hosp + 1) -70.3386 -129.2738 -11.4033 0.0193exposure:Comorbidity1 82.9098 -99.3891 265.2087 0.3727exposure:Comorbidity2 261.3797 -0.9835 523.7429 0.0509HC SE: 72.97, 37.28, 84.33, 213.08, 274.21, 41.80, 58.67, 55.61,55.20, 56.24, 57.90, 84.95, 30.07, 93.01, 133.86Wild boot: 73.98, 38.33, 84.70, 211.33, 274.10, 41.52, 58.22, 56.19,54.88, 56.85, 57.91, 84.95, 29.59, 93.11, 133.50156D.4. Pseudo-value Model O for Time Period 4Table D.4 shows the summary of pseudo-value Model O for Time Period 4.The spline for age shown in Figure D.14 has (roughly) the same shape as theearlier time periods. The effect was greater since the time period coveredby this model was 12 years (10 to 22 years post-index) versus two or threeyears. A higher SES was associated with more time spent out of the hospitaland alive. Being in the highest socioeconomic quintile was associated with181 extra days in state O compared to the lowest quintile. A comorbidity of2 while unexposed was associated with 211 fewer days in state O. More pre-index hospitalizations was also associated with less time in state O. Menexperienced about 183 fewer days in state O holding all other covariatesconstant.Figure D.14: Age at index date for Model O for Time Period 4The HC standard errors seem reasonable, and the residual plots (FiguresD.15 and D.16) did not indicate any serious problems with model specifica-tion: the residuals are centred at zero in each figure.157D.5. Pseudo-value Model H for Time Period 1Figure D.15: Residual plot for pseudo-value Model O for Time Period 4(a) Age at index date (b) Pre-index hospitalizationsFigure D.16: Residual plots for covariates in Model O for Time Period 4D.5 Pseudo-value Model H for Time Period 1The summary of Model H for time period 1 may be found in Table D.5.Figure D.17 shows the spline for age at the index date in Model H for TimePeriod 1. Being younger than 40 years old was not associated with anyeffect on time spent in state H between two and four years post-index. Theeffect was linear after age 40 and older age was associated with more daysspent in the hospital. Being male and having more pre-index hospitaliza-tions were associated with more days in the hospital. Being in the highestsocioeconomic quintile was associated with spending 1.85 fewer days in the158D.5. Pseudo-value Model H for Time Period 1hospital relative to the lowest quintile.As shown in Figures D.18 and D.19, the residual plots did not indicate anyserious problems with the model.Table D.5: Summary of Pseudo-value Model H for Time Period 1Beta LB UB p(Intercept) 0.4675 -0.9408 1.8759 0.5153exposure 0.1588 -0.6584 0.9760 0.7033ns(Index_Age, df = 3)1 4.5036 2.1234 6.8837 0.0002ns(Index_Age, df = 3)2 7.2688 1.9408 12.5969 0.0075ns(Index_Age, df = 3)3 11.7046 3.1499 20.2592 0.0073SEXM 1.0533 0.0842 2.0224 0.0331SES2 -0.7897 -2.4298 0.8504 0.3453SES3 -0.9989 -2.5736 0.5759 0.2138SES4 -1.2652 -2.8223 0.2920 0.1113SES5 -1.8589 -3.2796 -0.4382 0.0103Comorbidity1 0.9276 -0.4157 2.2709 0.1759Comorbidity2 0.4270 -1.6735 2.5275 0.6903log(Num_Prev_Hosp + 1) 2.8480 1.7343 3.9616 0.0000exposure:Comorbidity1 -0.5979 -3.6415 2.4458 0.7002exposure:Comorbidity2 -1.8570 -4.1692 0.4551 0.1154HC SE: 0.72, 0.42, 1.21, 2.72, 4.36, 0.49, 0.84, 0.80,0.79, 0.72, 0.69, 1.07, 0.57, 1.55, 1.18Wild boot: 0.72, 0.42, 1.21, 2.68, 4.31, 0.49, 0.84, 0.80,0.79, 0.72, 0.69, 1.09, 0.57, 1.54, 1.20159D.5. Pseudo-value Model H for Time Period 1Figure D.17: Age at index date for Model H for Time Period 1Figure D.18: Residual plot for pseudo-value Model H for Time Period 1160D.6. Pseudo-value Model H for Time Period 2(a) Age at index date (b) Pre-index hospitalizationsFigure D.19: Residual plots for covariates in Model H for Time Period 1D.6 Pseudo-value Model H for Time Period 2Table D.6: Summary of Pseudo-value Model H for Time Period 2Beta LB UB p(Intercept) 1.4545 -0.9098 3.8187 0.2279exposure 0.3853 -0.6173 1.3879 0.4513ns(Index_Age, df = 3)1 3.2527 -0.8390 7.3444 0.1192ns(Index_Age, df = 3)2 13.4897 2.6021 24.3773 0.0152ns(Index_Age, df = 3)3 24.0514 7.8219 40.2810 0.0037SEXM 2.1052 0.6710 3.5394 0.0040SES2 -1.4001 -3.8993 1.0991 0.2722SES3 -3.1549 -5.5398 -0.7700 0.0095SES4 -2.4040 -4.9341 0.1262 0.0626SES5 -3.3163 -5.6826 -0.9500 0.0060Comorbidity1 3.0953 0.6539 5.5367 0.0130Comorbidity2 4.5519 1.0010 8.1028 0.0120log(Num_Prev_Hosp + 1) 4.2806 2.4664 6.0947 0.0000exposure:Comorbidity1 -1.8270 -5.7830 2.1289 0.3653exposure:Comorbidity2 -6.5071 -10.6747 -2.3395 0.0022HC SE: 1.21 0.51 2.09 5.55 8.28 0.73 1.28 1.22,1.29 1.21 1.25 1.81 0.93 2.02 2.13Wild boot: 1.21, 0.52, 2.07, 5.52, 8.26, 0.73, 1.28, 1.22,1.30, 1.20, 1.24, 1.83, 0.92, 2.04, 2.14161D.6. Pseudo-value Model H for Time Period 2The spline of age at the index date for Model H for Time Period 2 (shownin D.20) has the same shape as for time period 1, and pre-index hospitaliza-tions was again associated with more time spent in the hospital. In ModelH for Time Period 2, higher socioeconomic status was associated with fewerdays spent in the hospital: individuals in the third, fourth, and fifth quin-tile spent around three fewer days in state H than individuals in the lowestquintile. The most noteworthy change from time period 1 is that comorbid-ity is associated with more time spent in state H.The residual plots (Figures D.21 and D.22) do not indicate any problemswith the model.Figure D.20: Age at index date for Model H for Time Period 2162D.7. Pseudo-value Model H for Time Period 3Figure D.21: Residual plot for pseudo-value Model H for Time Period 2(a) Age at index date (b) Pre-index hospitalizationsFigure D.22: Residual plots for covariates in Model H for Time Period 2D.7 Pseudo-value Model H for Time Period 3Comorbidity was associated with more time spent in the hospital between 7and 10 years post-index as shown in Table D.7. More pre-index hospitaliza-tions was again associated with more days in the hospital. Individuals over50 spent more time in the hospital than younger individuals. Figures D.23to D.26 show the residual plots for the model. The pre-index hospitaliza-tion residual plot in Figure D.26(a) may appear to indicate a problem butas shown in Figure D.26(b), the pattern indicated by the spline seems tobe the result of very few individuals with pre-index hospitalizations greater163D.7. Pseudo-value Model H for Time Period 3than 25. The other residual plots (Figures D.24 and D.25) did not indicateany serious issues.Table D.7: Summary of Pseudo-value Model H for Time Period 3Beta LB UB p(Intercept) 1.9736 -0.7574 4.7046 0.1567exposure 0.9977 -0.1711 2.1664 0.0943ns(Index_Age, df = 3)1 1.4854 -1.5297 4.5005 0.3343ns(Index_Age, df = 3)2 12.3192 4.1304 20.5080 0.0032ns(Index_Age, df = 3)3 24.8440 14.4076 35.2804 0.0000SEXM 1.1967 -0.0506 2.4440 0.0601SES2 -1.6963 -3.5006 0.1080 0.0654SES3 -1.3687 -3.1765 0.4390 0.1378SES4 -1.7459 -3.3674 -0.1244 0.0348SES5 -0.9940 -2.8987 0.9106 0.3063Comorbidity1 3.0454 1.1042 4.9866 0.0021Comorbidity2 5.9554 2.7029 9.2080 0.0003log(Num_Prev_Hosp + 1) 2.8489 1.2858 4.4120 0.0004exposure:Comorbidity1 -2.2826 -6.2370 1.6719 0.2579exposure:Comorbidity2 -0.9712 -7.4542 5.5119 0.7691HC SE: 1.39, 0.60, 1.54, 4.18, 5.32, 0.64, 0.92, 0.92,0.83, 0.97, 0.99, 1.66, 0.80, 2.02, 3.31Wild boot: 1.40, 0.59, 1.55, 4.10, 5.25, 0.64, 0.94, 0.93,0.83, 0.98, 1.00, 1.67, 0.78, 2.01, 3.28164D.7. Pseudo-value Model H for Time Period 3Figure D.23: Age at index date for Model H for Time Period 3Figure D.24: Residual plot for pseudo-value Model H for Time Period 3165D.8. Pseudo-value Model H for Time Period 4Figure D.25: Residuals against age for pseudo-value Model H for TimePeriod 3(a) Pre-index hospitalizations(b) Pre-index hospitalizations (fivelargest points greater than 25 removed)Figure D.26: Residual plots for pre-index hospitalization in Model H forTime Period 3D.8 Pseudo-value Model H for Time Period 4The pseudo-value regression model H for Time Period 4 is summarized inTable D.8 and the natural spline for age may be found in Figure D.27. Theeffect of age on days spent in the hospital after ten years post-index was lessclear than in the other time periods. It appears that between 40 and 60,the time spent in the hospital increased and plateaued after 60 years old. It166D.8. Pseudo-value Model H for Time Period 4also seems that most baseline variables are no longer as relevant in this timeperiod. Moreover, it appears that the time-dependent variables measuredat their landmark time of 10 years post-index were also not relevant for thelast few years of the study. The residual plots (Figures D.28 and D.29) donot indicate any serious issues with the model.Table D.8: Summary of Pseudo-value Model H for Time Period 4Beta LB UB p(Intercept) 14.1669 1.7316 26.6022 0.0256exposure 6.5374 0.8613 12.2136 0.0240ns(Index_Age, df = 3)1 32.1939 18.5239 45.8638 0.0000ns(Index_Age, df = 3)2 11.5912 -13.5000 36.6824 0.3652ns(Index_Age, df = 3)3 30.5475 5.0946 56.0005 0.0187SEXM 2.0995 -3.8719 8.0709 0.4907SES2 -0.0929 -7.5921 7.4064 0.9806SES3 1.4872 -8.7733 11.7477 0.7763SES4 -1.0817 -8.7767 6.6133 0.7829SES5 -2.7738 -10.5640 5.0165 0.4853Comorbidity1 -1.4948 -9.1358 6.1462 0.7014Comorbidity2 4.0781 -5.7139 13.8700 0.4143log(Num_Prev_Hosp + 1) 3.1419 -1.2628 7.5466 0.1621exposure:Comorbidity1 6.9361 -8.8157 22.6879 0.3881exposure:Comorbidity2 -9.2236 -26.7599 8.3128 0.3026HC SE: 6.34, 2.90, 6.97, 12.80, 12.99, 3.05, 3.83, 5.23,3.93, 3.97, 3.90, 5.00, 2.25, 8.04, 8.95Wild boot: 6.39, 2.92, 6.96, 12.92, 13.14, 2.99, 3.86, 5.28,3.93, 4.00, 3.94, 5.04, 2.24, 8.15, 8.76167D.8. Pseudo-value Model H for Time Period 4Figure D.27: Age at index date for Model H for Time Period 4Figure D.28: Residual plot for pseudo-value Model H for Time Period 4168D.9. Summary of covariate effects(a) Age at index date (b) Pre-index hospitalizationsFigure D.29: Residual plots for covariates in Model H for Time Period 4D.9 Summary of covariate effectsAge was the covariate with the strongest association to days spent in eitherstate. Older age was associated with fewer days in state O and more daysspent in state H. Living with comorbidities was often but not always as-sociated with more time spent in state H whereas it quite consistently ledto a shorter mean sojourn time in state O. Pre-index hospitalizations alsoled to more time spent in the hospital and fewer days in state O. Highersocioeconomic status was associated with more time in state O.169
Thesis/Dissertation
eng
Statistics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Petkau, A. John, 1950-
Zhao, Yinshan
The impact of disease-modifying drugs for multiple sclerosis on hospitalizations and mortality in British Columbia : a retrospective study using an illness-death multi-state model
Text