{'fraction_non_alphanumeric': 0.028719665271966528, 'fraction_numerical': 0.01502928870292887, 'mean_word_length': 4.9752, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multistate Models for the Analysis of Time to Type II Chronic Diabetic Complications in Debre Markos Referral Hospital, Northwest Ethiopia

Copyright © 2021 Tehran University of Medical Sciences. Published by Tehran University of Medical Sciences. This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International license (https://creativecommons.org/licenses/by-nc/4.0/). Noncommercial uses of the work are permitted, provided the original work is properly cited. *.Corresponding Author: Journal of Biostatistics and Epidemiology


Introduction
Diabetes mellitus (DM) is a chronic, noncommunicable disease characterized by elevated levels of blood glucose that occurs due to failure in secretion, action, or both of insulin. 1 It is commonly classified as type I (T1DM), type II (T2DM), or gestational diabetes, of which T2DM is the most common form. 2 T1DM, where the possible risk factors are autoimmunity, genetic, and environmental factors, accounts for 5% to 10% of all cases of diabetes, whereas type II diabetes occurs when the body cannot produce enough insulin or cannot use insulin, which accounts for 90% to 95% of all diagnosed diabetes cases. 3 T2DM may be unrecognized until complications become evident as the disease develops gradually. Both the number of new cases and the prevalence of diabetes have been steadily increasing over the past few decades. 4 According to the International Diabetes Federation report, nearly half a billion people live with diabetes worldwide and among these, almost (80%) of diabetes cases are found in low-and middleincome countries. This happens because of rapid urbanization, unhealthy diets, increasingly sedentary lifestyles, and inadequate resources to provide preventive or medical care for their populations. 5 This is a major cause of death in most countries. As a result, diabetic complications have a significant economic impact on countries, healthcare systems, and above all, individuals with diabetes and their families. [6][7][8] Ethiopia is challenged by the growing magnitude of chronic non-communicable diseases such as diabetes mellitus. 9 According to the 2017 International Diabetes Federation report, the number of people living with diabetes in Ethiopia was nearly 2.6 million, with an estimated prevalence of 3.8%. 6 The same report estimated that Ethiopia would be the 9th highest country globally, with 14.1 million diabetic cases by 2045. Diabetes is not only a serious condition but also a risk factor for both short-and long-term complications. Microvascular complications (MiVasC's) and macrovascular complications (MaVasC's) are associated with long-term diabetes-related complications that can result in morbidity and mortality. 2,8,10 Diabetic nephropathy, neuropathy, and retinopathy are the main types of MiVasC's. 11 Diabetes mellitus and its complications have increased over the years and are among the common reasons for inpatient admissions in Ethiopia. 12 It is the ninth leading killer in the country. 13 The prevalence of diabetes-related vascular complications has not been well identified and documented, although its impact continues to rise. In studying chronic diseases such as diabetic complications, clinical interest lies in both the final outcome (death) and the dynamics of the complications process. Multistate models are useful in describing the progression of diseases with several possible states over time. 14-18 These models provide an appealing framework to study the movement of disease progression across a certain number of states defined by specific disease conditions and their prognostic factors. 19,20 Multistate models were introduced in 1978 by Aalen and Johansen. 21 However, these models are not commonly used despite their potential applications. Although multistate models are widely applied in clinical studies and other fields, they are often overlooked by investigators in time-to-event analysis. 22 Compared with the Cox regression model with a single endpoint, multistate models provide a more detailed insight into the disease process. 23, 24 Survival models in the multistate framework provide methods to explore the effect of covariates on multiple responses, taking into account and evaluating the effects of intermediate events. Therefore, the purpose of this study was to explore how T2DM patients progress in different disease status, and to determine the prevalence and predictors of vascular complications as well as mortality status of type II diabetes mellitus (T2DM) patients based on a retrospective cohort study of patients in Debre Markos Referral Hospital (DMRH). Specifically, this research sought to answer the following questions: What is the prevalence of microvascular, and macrovascular complications among patients with T2DM? To what extent do the possible risk factors affect the transition intensity of T2DM to different clinical states? Is the effect of prognostic factors the same across transitions?

Study Area and Period
The study was conducted at the Debre Markos Referral Hospital (DMRH) in the town of Debre Markos. DMRH is located in the East Gojjam, Amhara National Regional State, which is 3000 KM in the northwest of Addis Ababa, capital of Ethiopia. The data for this study were collected from January 1, 2019, to March 30, 2019.

Study Design, Population and Period
A hospital-based retrospective study design was used on T2DM patients. All T2DM patients attending out-patient follow-up clinic at DMRH were taken as the source population. The study population was T2DM patients whose age was 18 years and above who have been on follow-up at DMRH from January 1, 2005, to December 31, 2017. All T2DM patients were retrospectively followed-up for any medical condition (such as developing vascular complications or death) from entry to the end of follow up. Once the individual diagnosed to have T2DM, the individual are required to visit T2DM clinics every month under normal conditions or absence of comorbidities. The T2DM patient was considered right censored when she/he did not experience the event of interest or lost to follow-up in the follow-up period.

Data Collection Tools and Procedures
A data abstraction sheet was developed based on information documented in the patients' medical records. To ensure the quality of the data, nurses working in the same hospital were trained on data collection procedures and data abstraction sheets. The members of the research team supervised the data collection process to verify the completeness of the data. The data were checked for completeness during data entry and cleaning processes. The information obtained from the patient medical cards was entered into Excel sheets and subsequently transferred for analysis into R software.

Study Variables
The primary outcome (state) was defined as the time taken to develop vascular complications (MiVasC or MaVasC) or death from the patients' first T2DM illness but free of any complications. Covariates such as patient sex, diabetes duration with T2DM, age at the onset of diabetes, residence (rural or urban), blood pressure (<140/90mmHg or ≥140/90 mmHg ), family history of diabetes, including fasting serum lipids (HDL, LDL, and triglycerides), glycated hemoglobin percentage (HbA1c%), serum creatinine level and serum electrolytes (sodium, potassium, and chloride levels) were 288 Vol 7 No 3 (2021) extracted. A patient was considered in the microvascular state if he/she developed at least one of the MiVasC's such as diabetic neuropathy, retinopathy, or nephropathy. Similarly, a patient was considered in the MaVasC state when a diabetic patient experienced an event of at least one MaVasC, such as peripheral artery disease, coronary artery disease, or stroke. The individual status of vascular complications was recorded only during clinical visits. Therefore, the exact date of an event is unknown (interval censoring). The right censoring time for each of the states (events) was measured for individuals who did not develop any of the complications during the follow-up period or were lost to follow-up.

Statistical Analyses
The baseline characteristics of the study participants were reported using descriptive statistics. A semi-parametric multistate model was applied to assess the relationship between the clinical conditions of patients and to determine the effects of covariates on individual states. The data were first entered into Microsoft Excel sheets and exported to R statistical software. 25 The mstate package in R software, which was developed by de Wreede et al., 14 and the survival package written by Terry Therneau and Lumley 26 were used to carry out all analyses and generate graphs.

State Structure for Diabetic Complications
States were distinguished by vascular complication type and death. To formally describe the state structure of this study, we let "MiVasC", "MaVasC", and "died" denote T2DM patients who experienced microvascular compilation, macrovascular complications, and death, respectively. The state structure, which specifies the states and possible transitions, is shown in Figure 1. A patient who was complication-free was at risk of developing micro-or macro-vascular complications, or died with or without developing vascular complications at some time during the retrospective follow-up period.

Missing Data Analysis
The potential risk factors associated with T2DM were included in the transition intensity of the model as covariates. However, some values of the covariate were missing due to different reasons. Ignoring the presence of missing data might produce biased estimates that lead to a false conclusion when missing data are not missing completely at random (MCAR), 27 and even if the data are MCAR, dropping a large proportion of the data results in a substantial loss of information. 28 To account for missing values in the covariates, we employed multiple imputations 28 to create completed data sets with missing information filled in. This is an established and broadly accepted method to enhance data sets with missing information. 29 There are several methods that can be used to handle missing data that are based on the assumptions of the type and degree of 'missingness' in the dataset. R packages including Amelia II (uses Bootstrap EM ), 30 Multivariate Imputation using Chained Equations (mice), 31 VIM, 32 and missForest 33 can be used to perform imputing missing data. For this analysis, the missForest R package was chosen as it was shown by Waljee et al.(2013) that it outperforms other well-known methods, and it is suitable for mixed-type data. 34 In missForest, missing values are imputed following the procedures and guidelines outlined in. 35

Multistate Model for Diabetic Complication data
Multistate models are routinely used in research where a change in status over time is of interest. They have a wide range of applications, notably in chronic disease progression studies such as chronic vascular complications in diabetes patients where the study units reversibly or irreversibly move through a succession of events or states. 15,24 It provides the opportunity to give a detailed description of the several states that individuals usually go through during a life course with type 2 diabetes. 36 Let Mi(t) be the assumed state of individual i, where i = 1, 2, . . . , n at the instant t = 1, 2, . . ., γ with respect to patient's T2DM complication status. Under the Markov assumption, the probability of an individual i moving from an h state to a j state in an interval δt given the patient's history (Hh−) is defined as: If state h is not persistent, then P hj (s, t) ≥ 0 for and P hj (s, t) = 0 for h=j. Otherwise, if state h is persistent, P hj = 0. The transition intensity between states h and j is defined by equation (2) (2)

Semi-parametric Multistate Model
The covariates in multistate models are often incorporated through transition intensity functions to explain differences among individuals during disease progression. One popular choice is the semi-parametric Cox regression model. 37 Under the proportional hazard assumption, the effect of a unit increase in a covariate is multiplicable with respect to the hazard rate. 38 According to de Wreede, Fiocco, 14 Andersen PK, 39 a time-homogeneous multistate model with proportional intensity having the transition-specific hazard, has the following form: (3) Here, λ ij0 (t) is the baseline i → j transition intensity, which is common for all individuals and specifies how the intensity changes with time, while β ij and Z are vectors of unknown regression coefficients and covariates, respectively, from the ith state to the j th state transition. The model specified in Equation (3) is a parametric model when the baseline hazard specified by λ ij0 (t) takes a functional 290 Vol 7 No 3 (2021) form. However, the covariates may be time-dependent, and the model becomes a nonhomogeneous Markov process. 40 The model in Equation (3) assumes that the process is Markovian and is known as a Cox Markov model. Several limitations are encountered when ignoring the disease history. Therefore, an alternative approach is to use Cox semi-Markov models in which the future of the process depends on the duration of the patient in the current state rather than the current time. 29 Such models are also referred to as "clock reset" in which the time t in λ hj (t) is considered the time duration from the entrance of the patient to the state h. 17 The semi-Markov process considers that the sojourn time of transition from one state to another depends on the time already spent in the current state, which is the most important assumption in medical data. 41 The semi-Markov model for the progression of diabetic complications is specified in Equation. 4 (4) where T ij is the sojourn time spent in the i th state and 1 ≤ i ˂ j ≤ 4. 40

Estimation Techniques
Under the Markov assumption, the nonparametric estimation of the transition probability matrix can be expressed as a function of the transition intensity in the form of a product integral 42 (5) where î j   Λ =   Ë is called the Aalen-Johansen estimator. The Aalen-Johansen estimator, which uses the mathematical framework of multivariate counting processes, reduces to simple empirical proportions for the complete data. Let N hj ,i be a counting process that counts the number of i's patient direct transitions (without visiting another state in between) from state h to state j up to time t, where h, j ϵ 1, . . . , 4 and h ≠ j. Here, time t represents the time measured from the initial state, which is the complication-free state. Let Yj (t) refer to the number of individuals to be observed at risk in a state just prior to time t. Then the conditional transition probability is specified using Equation (6). (6) where ∆N hj gives the number of h → j transitions observed exactly at time t. Kalbfleisch and Lawless 43 and later Kay R 44 described a general method for evaluating the likelihood of a general multistate model in continuous time, applicable to any form of transition matrix. The maximum likelihood that can be used to estimate the model parameters specified in Equation (3) is presented in equation (7) below. The partial likelihood function takes the product overall individuals and time points, as shown in equation (7). (7) where t hj,k is the failure or censoring time of individual k for transition h→j, d hj,k = 1, the individual k has an event for transition i → j, 0 otherwise, and Rh(t) is the risk set of state h at time t, that is, the set of individuals who are in state h at time t (t in this case is the time since the entry in state h). The estimate of the cumulative baseline hazard of transition i→j is the Nelson-Aalen estimate 45 (8) which is the proportion of the number of observed h → j transitions at tk to the number of individuals at risk in state h just prior to t k with h, j = r, h ƒ= j. In (4), the summation is over all event times tk, which are less than or equal to time t. Left truncation and right censoring are naturally accounted for in the denominator. The risk set in Equation (8) includes all individuals who have entered state h before time t and who have not yet moved out of state h again, provided that they have already entered the study (left-truncation) and are still under observation (right-censoring).

Baseline Demographic and Clinical Characteristics of the Study Participants
The characteristics of 524 T2DM patients who were followed retrospectively from January 1, 2005, to December 31, 2017, are presented in Table 1 to assess their time for vascular complications or the date of death, whichever occurred first. All patients were free of T2DM-related complications when they first entered into follow-up. The median follow-up time was 7.4 with an interquartile range (IQR) of 6.01-10.02, implying that 50% of the patients had been under followup for at least 7.4 years. According to the results in Table1, the majority (66.27%) of patients with T2DM were men. Microvascular complications were present in 19% (n= 33) and 18% (n=66) of female and male diabetic patients, respectively; while 23% (n=40) of female and 21% (n=70) of male patients developed MaVasC's; 13% (n=23) and 18% (n= 62) of female and male diabetic patients, respectively died due to diabetes-related complications. In addition, 10% (n=18) of female and 11% (n=37) of male patients developed both MiVasC and MaVasC. The percentage of females with T2DM who experienced MaVasC was higher than that of males. The percentage of female patients with diabetes who experienced MaVasC was higher than that of their male counterparts. This revealed that female patients with diabetes were more likely than their male counterparts to develop chronic diabetic complications than their male counterparts. However, men were more likely to die from diabetes-related complications than their female counterparts. The majority (66.6%) of the patients were from urban areas. As for the distribution of diabetic patients regarding their vascular complication status and their death status, the percentage of rural dwellers who developed MiVasC, MaVasC, and died was consistently higher than that of urban dwellers. However, the effects of sex and residence on the differences in proportions were not statistically significant. More than half of the patients, 278(58.3%), had baseline blood pressure values of ≥140/90mmHg, whereas the remaining 199 (41.7%) patients had <140/90mmHg baseline blood pressure values. Moreover,

Exploring Patterns of Missingness
The heatmap produced using vis_miss() from Visdat package 46 provides a pattern and percentage summary of missingness for each variable or column is shown in Figure 3. This plot helps the visualization of the pattern of percentages of missing values distribution. 47 It can be seen from Figure 3 that 8 covariates have a missing value with a very small proportion, whereas 7 variables including fasting serum lipids (HDL, LDL, total cholesterol), glycated hemoglobin percentage (HbA1%), and serum electrolyte (potassium, sodium, chlorine levels) have an extraordinary proportion of missing data. The proportion of missing data should not be used to guide decisions on multiple imputation. 48 Therefore, we decided to take out the seven covariates and rather work with the other covariates that have less proportion missing data and can be successfully imputed. 29 Imputation can be more precise if other variables are non-missing for those cases that are to be imputed. 49  It was also revealed in the figure that the patient had an increased probability to move into a micro-and macrovascular complication state from the beginning of the follow-up. However, the curve becomes stable for the fact that patients will understand the danger of the disease as time passes and modify their lifestyle, revealed that the probability of staying in the same state monotonically decreased with time, whereas the probability that DM patients move to the next state increased with time.

Multistate Model Results
Model Building and Assumption Assessment After imputation, a multistate Markov model was fitted to the imputed data. The four states of the multistate Markov model include complication-free, MiVasC, MaVasC, and death ( Figure 1). The death state is an absorbing state. To begin the model fitting process, we used the likelihood ratio test if the fitness of the model was improved after incorporating all sojourn times in different states. The likelihood ratio test showed that the fitness of the model to explain the data was significantly improved after incorporating the sojourn times, revealing that the Markov assumption was violated, and hence the transition hazards/intensities dependent on the time spent in the current state.
To improve the prediction of the transition probabilities, baseline predictors were considered. To select important predictors, covariates that were statistically significant at 10% level of significance in the univariable model with at least one transition were included in the higher (multivariable). Accordingly, age at the onset of diabetes, residence (rural or urban); blood pressure category (<140/90mmHg or ≥ 140/90 mmHg), family history of diabetes, triglyceride level, baseline hemoglobin level, and serum creatinine level were significant at the fixed level of significance in the univalent model. Hence, the remaining covariates were not included in the final model.
The assumption of proportional hazard or time homogeneity for the transitions needs to be checked to see if it is reasonably fulfilled for our data. To check the model assumption, both formal tests and graphical plots can be used. For this analysis, we used graphical plots (as shown in Figure 5) to check the assumption. The plots look fairly flat, revealing that the assumption of proportionality is not violated provided that the plots are decidedly different from flat. However, a formal test using the cox.zph  Table 3. Based on the output generated, all tests except for two covariates (residence and baseline creatine level) have shown non-violation of the proportional hazards assumption. The fitness, as well as the violation proportionality assumption, were improved when a model whose transition intensities to death (transitions specified by 4) are assumed to be proportional to the baseline hazards.

Checking the Markove Assumption
The duration of patients with T2DM (sojourn time in the diabetic complication-free state) was included as a prognostic covariate to test if it has a statistical significant relationship with the transition to death state from MiVasC or MaVas C. The fact that the sojourn time of patients had a significant association with the transition of the next transition of the model revealed an evidence that the Markov assumption is not valid. This calls for semi-Markov models as an alternative model for the analysis. The model's output is presented in Table 4 in which P-values and hazard ratio (HR) values of the predictors having statistically significant association (at 5% significance level) with the states are italicized. According to the model's results, keeping the effects of predictors in the model constant, an increased age at the time of diabetic diagnosis was associated with increased transition intensity into the death state. More specifically, a year increase in the age of patients at T2DM diagnosis corresponds to a 9% increase in the death of T2DM patients without experiencing any complications (or complications might not be detected). Moreover, a year increase in the age of T2DM patients at the time of diagnosis increases the risk of dying by a factor of 2.1 and 3.14 times for patients with microvascular and macrovascular, respectively, showing that an increased age at diagnosis corresponded with the risk of earlier vascular complications onset and death. The effect on the event of death was pronounced more after T2DM patients experienced vascular complications of either type. However, age at diagnosis showed no significant difference in the transition rate between type II diabetic patients in micro-or macrovascular complication states. Similarly, baseline blood glucose level was a statically significant predictor of the transition of patients into a microvascular complication state such that an increased baseline glucose level had an increased risk of developing microvascular complications. Furthermore, rural dwellers were more likely to move into the death state from the complication-free state (HR=2.1, P=0.01) and from the microvascular complication state (HR=1.31, P=0.05). Type II diabetic patients with a family history of diabetes had more likely to move from the macrovascular complication state into the death state. The presence of hypertension is a risk factor for the development of vascular complications or death. Hypertensive patients, for example, were 2.13 times more likely to move to the death state from the complication-free state compared to non-hypertensive DM patients after sole adjustment of predictors in the model. Besides, from the final model, we summarize that the time spent in the complication-free state (duration of diabetes) has a prognostic impact on the transition from macrovascular complication to death.

Discussion
Multistate models provide a flexible framework for understanding clinical events under consideration of the disease process as a whole, not only focusing on one single endpoint, like the classical Cox regression model. 23 In the applications of multistate, the dependence on the history of the process is negligible and gives sense to assume the Markov property, where the state sojourn times follow the exponential distribution with a constant hazard rate. 50 However, in the study of chronic disease, Markov assumption is unrealistic. For instance, in the study of T2DM disease progression, the disease stages may not follow a constant hazard rate. As a result, Semi-Markov models are applied for modeling a variety of phenomena in different areas, 51 such as the progression of breast cancer, 29 engineering, 52 and different chronic disease progression. 53 In this study, we used Semi-Markov model to model multistate chronic disease transitions by taking into consideration of the length of time that patients have stayed in a certain disease state, which is clinically more plausible. Semi-Markov models are applicable to degenerative diseases, while Markov models are applicable to epidemic diseases. 54 Therefore, a semi-Markov multistate model was adapted using the clock-reset approach to study the effects of different covariates on transit to the next state. During an average of 7.4 years follow-up period, nearly 28% of diabetic patients developed vascular complications of either type, which slightly higher than what has been previously reported by Wolde et al. 13 based on data from the University of Gondar Referral Hospital, Ethiopia, but lower than the study conducted in China. 55 The difference in the length of the follow-up period and the age groups of the study participants could be the reason for such discrepancy. Duration of diabetic illness was associated with an increased risk of death. The same findings were reported by. 56,57 An increased age at diagnosis corresponded with the increased likelihood of death. This feeding is in line with the findings by Kidanie, Alem, 12 Wolde, Atsedeweyen, 13 Zoungas, Woodward, 55 which was conducted in Ethiopia and Li, Chattopadhyay 56 in China. This could be due to the fact that aging results in relative insulin resistance, cell dysfunction, altered glucose metabolism, associated hormonal changes, and decreased physical activity, 58,59 which in turn leads to death. 56 In our study, women were more likely than men to develop chronic diabetic complications, although sex was not a significant predictor of diabetic complication status. This is because women are subjected to a partial loss of estrogen protection normally present in the premenopausal state. 60 In this study, being hypertensive was higher in the transition into microvascular and macrovascular complications as well as death state. This result is supported by Li, Chattopadhyay, 56 Agrawal, Ola, 61 and Tracey, McHugh. 62

Conclusion
Multistate models are good statistical methods for exploring predictors and multiple interrelated events by controlling their dependence. In this study, the semi-Markov model fitted the data well and could be taken as a convenient model for the analysis of time to diabetes-related complications or death in this particular data. Moreover, blood pressure, age at diagnosis of diabetes, and diabetes duration are were predictors of vascular complications and death. The percentage of diabetes-related complications and their consequences, clearly death, is not small, although the study needs to be confirmed by using an adequate sample. Thus, appropriate prevention and control strategies should be in place to limit further progress and impact. The current results could also be taken as a baseline for others who aim to explore the occurrence of vascular complications as a function of different risk factors in individuals with T2DM.

Conflict of interest
The authors declare that none of us has a financial or personal relationship with other people or organizations that could inappropriately influence or bias the content of the paper.

Ethics Approval and Consent to Participate
This study was approved by the Debre Markos University Research Ethics Committee and the medical director of the Hospital. As the study was retrospective, informed consent was not obtained from the study participants, but the anonymity and the confidentiality of the data were assured.

Data availability statement
The data set analyzed in this study is available