A mathematical model for predicting glucose levels in critically-ill patients: the PIGnOLI model

Design

The retrospective study encompassed analysis of a Multiparameter Intelligent Monitoring in Intensive Care II (MIMIC-II), a large clinical database of critically ill patients. Because the study utilized an open access clinical database, formal IRB approval was not required.

Data source

MIMC-II is a large registry of intensive care unit patients treated at Beth Israel Deaconess Medical Center, Boston, Massachusetts. Patient information on demographics, laboratory findings, imaging study, vital signs and progress notes was available (Saeed et al., 2011). MIMIC contains data on over 30,000 patients admitted during the period of 2001–2008. The database comprised varieties of ICUs, including the medical, surgical, coronary, and cardiac surgery recovery care units. ICU stays separated by less than 24 h were considered as one episode of ICU stay. Data were collected from electronic health data and all information produced during hospital stay were stored in the database.

All data were extracted by using structural query language (SQL) programming language from the database (Zhang, 2015). The Institutional Review Boards of the Massachusetts Institute of Technology (Cambridge, MA) and Beth Israel Deaconess Medical Center (Boston, MA) approved the establishment of the database. De-identification was performed to ensure patients’ confidentiality. Access to the database was approved after completion of the NIH web-based training course named “Protecting Human Research Participants” by the author Z.Z. (certification number: 1132877).

Selection of subjects

All adult patients were considered potentially eligible for this study. Patients actually received continuous insulin infusion were included. Children and neonates were excluded.

Outcomes

The primary outcome was retested BG level (mg/dl). The value, date and time of each BG were recorded in medical record. BG values can be those from blood chemistry and fingerstick, and the differentiation between venous glucose measurements and fingerstick levels was not performed in current analysis.

Clinical variables

Specific SQL programming languages for data extraction are shown in Supplemental Information 1. Comorbidities including diabetes, liver failure, congestive heart failure and renal failure were extracted because I felt that they may influence the sensitivity to insulin therapy. Laboratory parameters including bilirubin, C-reactive protein, serum creatinine and lactate were extracted. A total number of 1,117,076 BG measured with finger stick and in chemistry were extracted. A total number of 480,560 episodes of insulin rate were extracted.

Simultaneous use of intravenous (IV) total parental nutrition (TPN) and dextrose were extracted from the database. Different concentrations of dextrose were transformed to 5% dextrose (e.g., a volume of 10 ml 10% dextrose equals to 20 ml 5% dextrose). A total of 558,634 episodes of oral feeding (PO) containing glucose were extracted for its volume and time. All events were based upon charted time.

Data analysis

The objective of the analysis was to establish a linear regression equation between retested BG and insulin rate, controlling for other potential confounders. A data-driven approach means that the form of the equation was determined by data, depending on statistical significance. All variables thought to be associated with insulin sensitivity were extracted from the database and were considered for their inclusion in the model at outset.

I employed a multivariable fractional polynomial (MFP) method to construct the main effect model. The method combines backward elimination of statistically non-significant covariates with an iterative examination of the scale of continuous variables. MFP specifies two levels of significance levels: α1 = 0.15 for the test for exclusion and addition of variables to the equation and α2 = 0.05 to assess significance of fractional polynomial transforms of continuous variables. One degree of freedom was assigned to dichotomous variables and two-term fractional polynomials with 4 degrees of freedom were assigned to continuous variables. Continuous variables were modeled using closed test procedure, determining whether the covariate should be dropped from model at α1. Then, α2 = 0.05 was employed to test the need for transformation of the variable. The best two-term transformation was compared to the linear term by employing a closed test procedure. If the two term model is significantly better than the linear one at α2 = 0.05, the two term model is then compared to the one-term model. Otherwise, the linear term was retained in the model. Interactions were explored and terms with p < 0.05 were retained in the model.

The overall fit of the model was assessed by using R-squared, which is a reflection of the variance that can be explained by the model. Influential observations were evaluated by examining the leverage, Cook’s D and DFITS. Influential observations were excluded and the model was refitted by using MFPIGEN module. If the new model was significantly different from the original one, the original model would be reconsidered for the more parsimonious one. For example, some fractional transformation would be dropped and interaction terms could be dropped if the likelihood ratio test showed p > 0.05. R-squared values of the new parsimonious model would be compared to the original one to see whether the fitness was good enough. The whole dataset was split into two subsets: the training subset and the validation subset. Observed values of covariates were substituted into the fitted model to derive linear prediction. I then performed a regression model with linear prediction of the training subset as dependent variable and linear prediction of the validation subset as independent variable. The regression coefficient should be close to 1 and statistically significant at p < 0.05 if the model fits well to the validation subset.

All statistical analyses were performed by using Stata 13.1 (StataCorp, College Station, Texas, USA) and R software (R 3.1.1). Statistical significance was considered at p < 0.05.

Model exploration and development

The results of initial model fitting are shown in Table 1. The continuous variables including glucose, interval, dextrose rate and insulin rate were FP transformed and there were significant interactions between insulin rate and two terms of glucose. Glucose was transformed by two-term FP with the power of −0.5 and 1. Interval was transformed by two-term FP with the power of −2 and 1. Dextrose rate was transformed by one term FP with the power of 0.5. There were two interaction terms between insulin rate and glucose because glucose was modeled with two terms. The overall fit of the model was thought to be good with an adjusted R-squared value of 0.8449.

Influential observations were examined by using regression diagnostics (Supplemental Information 2). By excluding these influential observations, I refitted the model and found that FP terms and coefficients were remarkably changed (Table 2). Glucose was transformed by two-term FP with the power of −2 and 1. Interval was transformed by two-term power of 3 and 3. The results showed that the model was not stable, most likely due to complexity of the FP assignment and multiple testing during model fitting. The FP terms were influenced by several influential observations.

The parsimonious model was fitted to address the problem of instability. Graphical presentation showed that although the interaction term was statistically significant, the magnitude was of marginal clinical significance (Fig. 3 in Supplemental Information 2). Therefore, I opted not to incorporate interaction terms in the parsimonious model. Figure 1 shows the scatter points predicted by FP model and simple linear model, and the two lines were close to each other. Visual inspection of the graph indicates the use of parsimonious model would not compromise the prediction accuracy of the model.

Final model and model validation

The final Predicting Glucose Levels in ICU (PIGnOLI) model is shown in Table 3. Insulin rate was significantly associated with blood glucose reduction (coefficient: −0.52, 95% CI [−1.03, −0.01]). Initial blood glucose was the most important determinant of retested blood glucose (coefficient: 0.89, 95% CI [0.88, 0.90]). Oral intake, TPN and dextrose infusion were all associated with blood glucose control. Furthermore, serum lactate and BMI were positively associated with retested blood glucose. The time interval was negatively associated with retested blood glucose level (coefficient: −0.18; 95% CI [−0.22, −0.14]). The PIGnOLI model showed an adjusted R-squared value of 0.84, which was not significantly different from the FP model with interaction terms (R-squared = 0.84). The PIGnOLI model was tested in the validation subset and the result showed that the coefficient between estimated retest glucose and observed retest glucose was 0.99 (95% CI [0.97–1.01]; p < 0.001). The adjusted R-squared value was 0.8259, suggesting that the model was well calibrated with the validation subset.