A study on the factors influencing the transfer of COVID-19 severe illness patients out of the ICU based on generalized linear mixed effect model

1.   Introduction

SARS-CoV-2, a novel coronavirus pneumonia, was discovered in Wuhan, China, in December 2019 ^[1,2], and it spread fast to other parts of the world. The World Committee on Virus Classification gave it the formal designation SARS-CoV-2 ^[3], and the World Health Organization designated the virus-caused illness COVID-19 ^[4]. The People's Republic of China's Law on the Prevention and Control of Infectious Diseases has classified COVID-19 as a category B infectious disease since it is a highly contagious acute respiratory infectious disease, and it is treated as such ^[5]. Serious COVID-19 instances can cause severe pneumonia, acute respiratory distress syndrome (ARDS), and multiple organ failure that can result in death, but mild cases only cause typical respiratory system infection symptoms ^[6]. The majority of COVID-19-infected patients some of whom may be asymptomatic are the source of the virus. According to available data, COVID-19 can spread from person to person by respiratory droplets, direct touch, and even fecal-oral transmission ^[7]. Additionally, it can be spread through aerosolization if exposed to a high aerosol concentration for an extended period of time. Generally speaking, the SARS-CoV-2 virus can infect people of any age ^[8].

There might not yet be a treatment for COVID-19 that works. Numerous techniques are required to stop and limit the epidemic's progress. On the one hand, prioritizing strong source control, personal protection measures, early diagnosis, and isolation. On the other hand, a focus on symptomatic treatment is also crucial to raise the disease's cure rate and lower the mortality rate. As a result, treating COVID-19 patients who are severely ill and in critical condition is very vital and challenging. Clinical traits of COVID-19 patients have been published so far ^[9,10].

Patients with COVID-19 can be categorized into severe and non-severe groups based on their clinical symptoms. Older age, male sex, and comorbidities are risk factors for COVID-19, according to a number of retrospective cohort studies from Wuhan ^[11,12,13]. There have also been reports of cohorts from other nations ^[14,15]. Although all survivors were pooled together in the majority of studies that evaluated risk factors, there were notable disparities between survivors of non-severe and severe groups.

Severe and critical subgroups could be created from the severe group. Patients who are considered severe or critical must be transferred to the ICU for treatment because they have a quickly progressing course and may experience multiple organ dysfunction, including respiratory failure ^[16].

However, some COVID-19 patients cannot be transferred out of the ICU in a short time after effective treatment. By using the Laplace approximation parameter estimation method based on maximum likelihood estimation, the generalized linear mixed effect model (GLMM) was established to analyze the characteristics of clinical indicators in critical patients, and to screen the main influencing factors of COVID-19 critical patients's inability to be transferred out of the ICU in a short time.

2.   Material and methods

2.1.   Study design and period

An institution-based, two-month retrospective follow-up study was conducted among COVID-19 critical patients admitted to the ICU to start treatment in Urumqi, Xinjiang Uygur Autonomous Region between July 10, 2020 and August 31, 2020.

2.2.   Inclusion criteria

According to the Diagnostic Criteria of COVID-19 Diagnosis and Treatment Protocol (Trial Edition 8) issued by the National Health Committee of the People's Republic of China on August 19, 2020, the confirmed cases must conform to the clinical manifestations in the diagnosis and treatment protocol, and the pathogenic or serological test is positive ^[17]. Common type: having fever, respiratory symptoms, etc., and imaging findings of pneumonia. Adult-severe: Adults meeting any of the following conditions: shortness of breath, Respiration Rate (RR) $\geq$ 30 beats/min; In the resting state, finger oxygen saturation $\leq93\%$ when inhaling air; The arterial oxygen partial pressure (PaO$_{2}$)/oxygen concentration (FiO$_{2}$) $\leq$ 300 mmHg (1 mmHg = 0.133 kPa); In areas with high altitude (over 1000 m above sea level), PaO$_{2}$/FiO$_{2}$ should be corrected according to the following formula: PaO$_{2}$/FiO$_{2}$$\times$[760/ mmHg)]; Progressive worsening of clinical symptoms and significant progression of lesions greater than 50% within 24 to 48 hours on lung imaging are managed as severe cases. Critical type: patients meeting one of the following conditions: respiratory failure requiring mechanical ventilation; shock; combined with other organ failure requiring ICU monitoring and treatment.

2.3.   Data

The demographic characteristics, epidemiological data, complications, clinical manifestations and laboratory examinations of 76 patients with COVID-19 admitted to the ICU ward of a hospital in Urumqi, Xinjiang Uygur Autonomous Region from July to August, 2020 were collected.

Longitudinal data refers to the data obtained by tracking the evolution of a series of test individuals over time. In this study, after COVID-19 critically ill patients were transferred to the ICU, they should undergo regular laboratory examinations to measure their clinical indicators for better treatment. The changes of clinical indicators of 76 COVID-19 patients during their transfer to ICU were recorded in this study. The number of test times per COVID-19 patient ranged from 10 to 25 times, and a total of 1097 data could be used. Figure 1 shows the changes of clinical indicators of one COVID-19 patient over time within 39 days after being transferred to the ICU.

2.4.   Descriptive statistics

For the description of baseline characteristics, count data in this data included Sex and Previous Underlying Diseases, and descriptive statistics were performed using frequency and composition ratios for count data. The measurement data included Age and Laboratory Examination Data, and the mean $\pm$ standard deviation was used to statistically describe the measurement data that conformed to a normal distribution. The results of the baseline characteristics were described as follows.

Among them, there were 38 males (50%) and 38 females (50%), whose age ranged from 24 to 83, with an average age of 61. There were 12 mild-to-moderate cases, 53 severe cases and 11 critical cases. See Table 1 for the composition of the patients with previous underlying diseases.

Clinical symptoms include fever, cough, expectoration, chest tightness, shortness of breath, muscle soreness, sore throat, dry mouth, headache, fatigue, anorexia, diarrhea, frequent micturition, weak consciousness, nausea and vomiting, etc. Laboratory examinations include routine blood, liver and kidney function, electrolyte, myocardial enzyme spectrum, inflammation indicators, and etiology detection. The laboratory examination data at the time of initial transfer to the ICU are shown in Table 2.

2.5.   Variable selection

The traditional methods of variable selection are the forward introduction method, the backward elimination method, and the stepwise regression method based on test statistics ^[18,19]. Usually, these three methods yield the desired model, but still have shortcomings. Tibshiranilea ^[20] pointed out that the variable selection process of these three methods is a discrete process, and there is a large variability in the selected models, i.e., small changes in the data may lead to the selection of a different model. He combined the advantages of the subset selection method and ridge regression and proposed the Lasso selection method, which accomplishes the selection of variables while continuously compressing the coefficients. In order to improve the accuracy of the model in this study, the Lasso selection method was also used in the process of variable selection in this study.

2.6.   Generalized linear mixed effect model (GLMM)

Generalized linear mixed effect model (GLMM) has a great application prospect in the medical field, which is obtained by adding random effect expansion to the fixed effect of generalized linear model, (GLM). Since the response variable of the generalized linear mixed effect model is the conditional mean function under the condition of the given environmental variable $X$ random effect $u$, it can be obtained as follows:

In Eq (1), $g(y)$ is the Link Function; $\beta$ is the Fixed Effects; $u$ is the Random Effects.

2.7.   Parameter estimation

The parameter estimation and performance of the generalized linear mixed effect model have attracted much attention, whose parameter estimation can be calculated by the maximum likelihood method, and whose likelihood function can be obtained according to the maximum likelihood method:

In Eq (2), $f(y|u)$ is the density function of response variable under the condition of random effect $u$, and $p(u)$ is the density function of random effect u. Due to the nonlinear relationship between response variables and independent variables and the existence of random effect u in the model, it is necessary to solve the problem of high-dimensional integration of random effect u. Many approximate inference methods by maximizing likelihood function have been proposed accordingly. At present, the commonly used methods are: Penalized Quasi-Likelihood (PQL) ^[21] method, Laplace approximation method and Gauss-Hermite Quadrature (GHQ) method.

In this study, Laplace approximation method based on maximum likelihood estimation is adopted, which is commonly used to estimate high calculus. Breslow and Lin (1995) used fourth-order Laplace approximation to establish a random effect model with only one random effect in each cluster. Raudenbush et al. (2000) extended this method to higher-order approximation and models with multiple related random effects in each cluster (repeated measurements for each patient) ^[22,23]. When using R software, the glmer function in package "lme4" ^[24] takes Laplace Approximation as the default method for regression analysis of generalized linear mixed effects model, but the glmer function will not converge when random effects are added into explanatory variables. At this time, the glmmTMB function in "glmmTMB" ^[25] package can be called to solve this problem, or the GLMM can be calculated by Laplace approximation.

2.8.   Intra-class correlation coefficient (ICC)

In GLMM multilayer statistical analysis model, another important concept is ICC (intra-class correlation coefficient) internal correlation coefficient, or called as intra-group correlation coefficient. According to the classical definition (Shrout $&$ Fleiss, 1979) ^[26], ICC is defined as the ratio of inter-group variance to total variance:

In Eq (3), $\sigma_t^2$ represents inter-group; $\sigma_h^2$ represents intra-group; $(\sigma_h^2+\sigma_t^2)$is the total variance.

In this study, since the response variable is whether COVID-19 critical patients are transferred out of ICU, and the response variable is the two-class discrete variable with (0, 1) distribution, and the fixed effect obeys logistic distribution, the intra-group variance $\sigma_h^2$can be expressed as $\pi^2/3$, and the inter-group variance $\sigma_t^2$ can be expressed by the random effect $\sigma_u^2$. Therefore, the intra-group correlation coefficient of response variable measurement in this study is as follows ^[27]:

Since ICC reflects the correlation degree among individuals in random effects, if ICC = 0, it means that there is no hierarchical structure among individuals in random effects, the random effects can be ignored and simplified to the traditional single-layer model; if ICC $>$ 0, it means that there is a hierarchical structure among individuals in random effects. If the random effect part is added, the mixed effect model can be used for analysis.

2.9.   Time-dependent covariates cox-proportional-hazard model

As a complement to the GLMM results and to test the stability of the GLMM model, , we also developed the Time-Dependent Covariates Cox-Proportional-Hazard model. A time-dependent Cox model is a Cox regression model with time-dependent covariates, also known as an extended Cox model. In a general Cox model, the covariates are fixed and do not change over time; a time-dependent Cox model takes into account the dynamics of the covariates over time and allows analysis of the case where the covariates are panel data.

The Cox model presented can be extended to handle exogenous time-dependent covariates using the counting process formulation investigated in detail by Andersen and Gill ^[28] and presented in greater detail in Fleming and Harrington ^[29] and Andersen et al. ^[30]. The time-dependent Cox regression model is an implementation that views the occurrence of an event as a very slow Poisson process. In counting process notation, the event process for subject i is written as $\{N_i (t), R_i (t)\}$, where $N_i (t)$ denotes the number of events for subject i before time t, and $R_i(t)$ is a left-sided continuous risk process with $R_i (t)$ = 1 if subject i is at risk at time t and $R_i (t)$ = 0 otherwise. The time-dependent Cox regression model is written as

where, in Eq (5), $\omega_i$ denotes a vector of baseline covariates, and $Y_i (t)$ denotes a vector of time-dependent covariates. Moreover, note that since $Y_i (t)$ is time-dependent, model no longer assumes that the hazard ratio is constant in time.

Estimation of $\gamma$ and $\alpha$ is based on the corresponding partial log-likelihood function that can be written as

3.   Model determination and result analysis

3.1.   Variable selection results

In this study, a total of 76 COVID-19 patients were recorded for changes in a total of 15 clinical indicators over time during their transfer to ICU. To improve the reliability and validity of our model, Lasso analysis was performed on 15 variables (Figure 2). Using the $\lambda$ with the Minimum Mean Standard Error (Log($\lambda$) = -4.025) method of Lasso analysis, 5 variables Age, CD8+, CRP, Cr, and LDH were finally selected for inclusion in the model for subsequent modeling, and their coefficients are shown in Table 3.

3.2.   Model determination

Using R (4.0.5) software, the glmmTMB function in "glmmTMB" ^[25] package was called to estimate the parameters of the generalized linear mixed effect model and analyze the fitting effect. In order to better screen the influencing factors of the transfer of COVID-19 critical patients out of ICU, it was necessary to determine an optimal model, and the following four steps could be carried out:

Step 1: Fit models with only random effects to calculate ICC, thereby determining whether the response variable could be analyzed using a generalized linear mixed-effects model.

In Eq (7), $\eta_{ij}$ was the response variable, and $u_{0j}$ was the random effect. Wherein, the response variable was the (0, 1) distribution variable of whether to transfer out of the ICU, and the random effects were COVID-19 critical patients in each state of transferring to and out of ICU. The parameter estimation results were obtained through the analysis of glmmTMB function of R software. Since the ICC value could not be calculated by glmmTMB function, the ICC function in "performance" ^[31] package was called to calculate the ICC value, and finally ICC = 0.24 was obtained.

From the obtained ICC value $>$ 0, it could be concluded that individuals in random effects had hierarchical structure, and GLMM model could be adopted.

Step 2: On the basis of Step 1, add explanatory variables to the fixed effect to obtain the random intercept Model 1:

In Eq (8), $\beta_0 +\beta_1x_{1ij}+\beta_2x_{2ij}+\cdots$ was the fixed effect part added with explanatory variables. After Lasso Cox analysis, the final fixed effect part included 5 variables, i.e., Age; CD8$+$; CRP; Cr and LDH. See Table 4 for the calculation results of the fixed effect by software glmmTMB function.

Step 3: On the basis of Step 2, add explanatory variables to the random effect to obtain the random slope and random intercept Model 2:

In Eq (9), $u_{1j}x_{1ij}$ was the explanatory variable added to the random effect, and the explanatory variable added in this study was obstime (July 10, 2020 as the observation start date and August 31, 2020 as the ending date, obstime is the time difference between the time of each measurement and the observation start date after the patient is admitted to the ICU until the end of the observation period). Since obstime was the time in each state of each COVID-19 critically patient's transferring to and out of ICU, and one patient corresponded to multiple times, it is included in the explanatory variable of the random effect. See Table 5 for the calculation results of random effects of Models 1 and 2 by the glmmTMB function of R software.

Step 4: Compare Model 1 with Model 2 to determine the optimal model and test.

The ANOVA function of R software was called to compare whether Model 2 with the addition of random effects is better than Model 1 without the addition of random effects. After performing Anova analysis on the two models, the AIC (Akaike Information Criterion) ^[32] of the Anova results were compared to select the optimal model. The results are shown in Table 6.

As can be seen from Table 6, when Models 1 and 2 differ by 2 degrees of freedom, the probability value of chi-square test is P $<$ 0.001, which is statistically significant. It can be considered that there is a significant difference between the two models. Since the AIC value of Model 2 is smaller than that of Model 1, it can be determined that Model 2 (Eq 9) is the best model, that is, the data can be better fitted after the random effect is added to the explanatory variables.

3.3.   Superiority test and stability test of the GLMM model

After the model was determined, the residual analysis of the model (Figure 3) and the normality test of random effects (Figure 4) were carried out by calling the binnedplot function in the "arm" ^[33] package and the qqplot function in the "car" ^[34] package in R software. As can be seen from Figure 3, the average residual value basically has a constant and uniform diffusion within the expected value range in a symmetrical mode, indicating that the overall fitting effect of the model is good. Figure 4 shows that both the random effects and the added explanatory variables obey normality.

In order to test the stability of the GLMM model, we also built the Time-Dependent Covariates Cox-Proportional-Hazard model using the "coxph" function in the "survival" ^[35] package. By comparing the AIC of these two models, we can see that the AIC of the GLMM model (AIC = 1018.3) is smaller than that of the Time-Dependent Covariates Cox-Proportional-Hazard model (AIC = 2147.3), so the GLMM model is better.

In addition, by comparing the results of these two models, it can be seen that the results of the GLMM model (Table 7) are similar to those of the Time-Dependent Covariates Cox-Proportional-Hazard model (Table 8), which indicates that the stability of the GLMM model is good.

3.4.   Result analysis

In our data collection process, each patient was clinically tested every day after he/she was admitted to the ICU, and the data of each clinical index was collected, and the value of "whether to transfer out of ICU" was 0 until the value of "whether to transfer out of ICU" became 1 after he/she was transferred out of the ICU. So the 0 means not to transfer out of ICU, and 1 means to transfer out of ICU.

Through the analysis of GLMM model 2, the results of influencing factors that affect COVID-19 severe illness patients can't be transferred out of intensive care unit in a short time are obtained, as shown in Table 7. From Table 7, it can be seen that there are four influencing factors with p-values less than 0.05: Age, CRP, Cr and LDH, among which Age, CRP and LDH have OR values less than 1, indicating that these three factors are protective factors, and the estimated values of these three factors are negative. On the contrary, Cr is the risk factors with OR value greater than 1, which indicates that with the increase of this factor, the probability of COVID-19 critically ill patients being transferred out of the intensive care unit in a short time is greater (Figure 5).

4.   Discussion

During the preliminary data collection, it was found that COVID-19 patients who were admitted to ICU were relatively older, with an average age of 62. This is because the immunity of middle-aged and elderly people was weak, and as shown in Table 1, that more than half of these elderly patients were accompanied by chronic basic diseases. Studies show that about 39.0–40.8% of elderly COVID-2019 patients suffer from hypertension, 16.0–17.7% suffering from diabetes, 4.4–15.7% suffering from cardiovascular diseases, 4.4–6.2% suffering from chronic basic diseases of liver and kidney, and 2.2–6.2% suffering from chronic obstructive pulmonary disease ^[36,37,38]. Therefore, after novel coronavirus infection, the human body's functional load is heavy and the immune function is weak, which leads to the aggravation of the disease, and further admission to the intensive care unit. The results also reveal that the younger the age, the greater the probability of being transferred out of ICU. The results still suggest that the inflammatory factors and immune function of COVID-19 patients need close attention during the treatment in ICU, such as the increase of C-reactive protein index indicates a relatively high inflammatory response; the decrease of lymphocyte count indicates possible immunosuppression ^[39]. Therefore, the changes of COVID-19 patient's condition during the treatment in ICU should be closely observed, especially the elderly patients. As the influencing factors, CRP, Cr, LDH and other indicators need to be paid close attention to, and active treatment should be taken to make the condition better, thereby transferring out of the ICU more quickly.

In addition, according to the results of the study, Cr is risk factor, indicating that as the levels of this factor increased, the probability of transferring patients with COVID-19 intensive care unit within a short period of time increased. Three factors, namely, "Age", "CRP", "LDH", were protective factors, and the estimated values of these three factors were negative values, indicating that as the levels of these three factors decrease, the greater the probability of transferring COVID-19 critically ill patients out of the intensive care unit within a short period of time. In the process of transferring COVID-19 patients to the ICU for treatment, clinicians can analyze the causes of decreasing Cr indexes and increasing Age, CRP, and LDH indexes in patients, and then can quickly diagnose the causes and treat them, so that patients can be transferred out of the ICU within a short period of time.

This study also has some limitations in that we did not collect Covid-19 death-cases admitted to ICU during data collection. Because there were no deaths in Urumqi during this outbreak according to the relevant reports of the Health Commission of Xinjiang Uygur Autonomous Region from July to August 2020, this resulted in missing information on death-cases admitted to ICU with Covid-19. The absence of death information may lead to some bias in the results, and in further studies later in this study, we will collect more data containing information on deceased patients to further improve the applicability of the model.

To sum up, this study mainly introduces the theoretical framework of generalized linear mixed effect model (GLMM), the method of parameter estimation, and the steps to realize and determine the optimal model by R (4.0.5) software. By analyzing and screening the factors influencing COVID-19 critical patients' inability to be transferred out of the ICU in a short time, the generalized linear mixed effect model is applied to the medical field for analysis. Since mutually independent data are often encountered in the research of medical field, the GLMM model is constructed to address the limitation of data independence required by the generalized linear model and the linear mixed model. The GLMM model can better analyze broader data types such as longitudinal data, time autocorrelation data and the like in the medical domain and thereby has a very good application prospect in the medical research field.

Author contributions

All authors conceived the study, carried out the analysis, discussed the results, drafted the first manuscript, critically read and revised the manuscript, and gave final approval for publication. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

The authors would like to thank the editor and the referees for their helpful comments. This research was supported by the Natural Science Foundation of Xinjiang (Grant No. 2021D01C268), the Program for Tianshan Innovative Research Team of Xinjiang Uygur Autonomous Region, China (2020D14020) and the National Natural Science Foundation of China (Grant No. 11961071).

Conflict of interest

The authors declare no conflict of interest.