NQF1626 Final Report
E.T. Moseley
Background
This is a report focusing on end of life care preferences and their implementation by the clinical care team. In particular, we will focus on the National Quality Forum Measure #1626:
Percentage of vulnerable adults admitted to ICU who survive at least 48 hours who have their care preferences documented within 48 hours OR documentation as to why this was not done.
For our purposeses:
Care preference documentation within 48hrs of admission, age
>74.
Care provider/Patient interactions from the MIMIC-III Intensive Care Unit Database are analyzed within the framework of a hospital admission utilizing a series of Mixed Effects Models. The patient is the random effect, so we can explain the difference in variance explained by patients compared to that explained by other effects (when the patient is introduced to other demographic and clinical variables) in the model.
The use of fixed and random effects is deliberate. A fixed effect is used when a categorical variable contains all levels of interest. A random effect is used when a categorical variable contains only a sample of all levels of interest, the samples which do not exist are unobserved variables and, since their distributions do not exist, they will be estimated. The levels of the random effect are assumed to be independent of each other; this assumption requires a priori knowledge of the phenomenon being modeled.
Sequential Mixed Model Design
Model 0 (No Covariates)
- Care Provider (Random)
Model 1a (Patient Demographics)
- Gender
- Age
- Ethnicity
- Language
- Marrital Status
- Care Provider (Random)
Model 1b (Clinical Variables)
- Sequential Organ Failure Assessment (SOFA)
- Care Provider (Random)
Model 2a (Demographics + ICU)
- Demographics (Model 1a)
- First Intensive Care Unit
- Care Provider (Random)
Model 2b (Clinical Variables + ICU)
- SOFA (Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
Model 3 (Demographics + Clinical Variables + ICU)
We will use the Sequential Organ Failure Assessment (SOFA) from Illness Severity Scores (github) to determine the severity of illness at admission.
- Demographics (as in Model 1)
- First Intensive Care Unit
- SOFA Score (As Model 2)
- Care Provider (Random)
Statistics of Interest
Intra Class Correlation calculated according to Wu, S. et al. (2012)
Standard Deviation of Random Effects calculated using
lme4, see: Bates, D. et al. (2015)Marginal (variance explained by fixed effects) R^2
\[R_{GLMM(m)}² = \frac{(σ_f²)}{(σ_f² + σ_α² + σ_ε²)}\]
- Conditional (variance explained by the entire model) R^2
\[R_{GLMM(c)}² = \frac{(σ_f² + σ_α²)}{(σ_f² + σ_α² + σ_ε²)}\]
Calculated according to Nakagawa, S. et al. (2017)
Libraries and Utility Functions
library("ggplot2") ## Graphing
library("scales") ## Graphing
library("GGally") ## Extension to ggplot
library("lme4") ## for Hierarchical Models
library("sjPlot") ## Lovely Presentation of Model Output
library("rcompanion") ## Pairwise nominal test
library("ROCR") ## ROC Curves
library("pROC") ## ROC Curves
library("MuMIn") ## MuMIn::r.squaredGLMM()Attending check
attending_check() will accept a data frame and return a data frame containing only those hospital admissions where an attending documented a note.
attending_check <- function(dat){
## Temporary data frame
tmp_frame <- data.frame()
## Results frame
res <- data.frame()
## For each hospital admission
for (name in unique(dat$HADM_ID)){
## Subset admission
tmp_frame <- dat[dat$HADM_ID == name, ]
## If any care providers are "Attending"
if (any("Attending" %in% tmp_frame$CG_DESCRIPTION)){
## add hospital admission to results
res <- rbind(res, tmp_frame)
}
}
## Return control to outer level
return(res)
}Provider Caremeasure Check
provider_caremeasure_check() will accept a data frame and, if a care provider implemented the care measure, that provider will be considered as having implemented the care measure for the entirety of the hospital admission.
provider_caremeasure_check <- function(dat){
## Temporary hospital admission data frame
tmp_frame <- data.frame()
## Temporary careprovider data frame
cg_frame <- data.frame()
## Results frame
res <- data.frame()
## For each hospital admission
for (name in unique(dat$HADM_ID)){
## Subset admission
tmp_frame <- dat[dat$HADM_ID == name, ]
## For each care provider
for (id in unique(tmp_frame$CGID)){
## Subset provider
cg_frame <- tmp_frame[tmp_frame$CGID == id, ]
## if care measure were implemented during the admission by the provider
if (any(cg_frame$NQF == 1)){
## Apply credit for entire admission
cg_frame$NQF <- rep(1, each = nrow(cg_frame))
}
## Bind
res <- rbind(res, cg_frame)
}
}
## Return control to outer level
return(res)
}Team Caremeasure Check
team_caremeasure_check() will accept a data frame and, if any care provider implemented the care measure, all providers in the care team will be considered as having implemented the care measure for the entirety of the hospital admission.
team_caremeasure_check <- function(dat){
## Temporary hospital admission data frame
tmp_frame <- data.frame()
## Results frame
res <- data.frame()
## For each hospital admission
for (name in unique(dat$HADM_ID)){
## Subset admission
tmp_frame <- dat[dat$HADM_ID == name, ]
## If ANY care provider documented care preferences
if (any(tmp_frame$NQF == 1)){
## Apply credit for entire admission
tmp_frame$NQF <- rep(1, each = nrow(tmp_frame))
}
## Bind
res <- rbind(res, tmp_frame)
}
## Return control to outer level
return(res)
}Expire Check
expire_check() will accept a data frame, and if the patient died in-hospital at any point that fact will be documented.
expire_check <- function(dat){
## Create Caremeasure Variable
dat$HOSP_DEATH <- rep(0, each = nrow(dat))
## Temporary data frame
tmp_frame <- data.frame()
## Results frame
res <- data.frame()
## For each hospital admission
for (name in unique(dat$HADM_ID)){
## Subset admission
tmp_frame <- dat[dat$HADM_ID == name, ]
## If any admission has a death
## Check if pt expired in hospital
if (any(tmp_frame$HOSPITAL_EXPIRE_FLAG == 1)){
tmp_frame$HOSP_DEATH <- rep(1, each = nrow(tmp_frame))
}
## add hospital admission to results
res <- rbind(res, tmp_frame)
}
## Return control to outer level
return(res)
}rocplot
rocplot() is a simple utility function for generating ROC curves.
rocplot <- function(pred, truth, ...) {
predob = prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf, ...)
area <- auc(truth, pred)
area <- format(round(area, 4), nsmall = 4)
text(x=0.7, y=0.1, labels = paste("AUC =", area))
# the reference x=y line
segments(x0=0, y0=0, x1=1, y1=1, col="gray", lty=2)
}model_info
model_info() is a simple utility function for summarizing model output and generating confidence intervals (95%).
model_info <- function(fit){
#Summary info
model_sum <- summary(fit)
#Odds ratio, confidence interval
odds_ratio <- cbind(OR = exp(fit$coef), exp(confint(fit)))
#Create list for return
my_list <- list(model_sum, odds_ratio)
#names
names(my_list) <- c("Model Summary","OR Summary")
return(my_list)
}Data Loading
dat <- read.csv("~/nqf_caregivers/data/nqf_study_data10Sep18.csv", header = T, stringsAsFactors = F)NQF Care Measure Check
Check for attending, then care measure implementation, then hospital expiration.
## Provider level
prov <- expire_check(provider_caremeasure_check(attending_check(dat)))
## Team level
team <- expire_check(team_caremeasure_check(attending_check(dat)))
Provider-Level Analysis
At this level, we will see if providers documented patient care preferences at least once in the first 48hrs of the hospital admission. Here, the clinician/patient interaction is the unit of analysis.
Initial Rates of Documentation Check
Here we will aggregate the care providers and ICUs to find care providers’ average rates of documentation.
tmp <- aggregate(NQF ~
CGID + FIRST_CAREUNIT,
data = prov,
FUN = mean)
par(mfrow=c(2,3))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "CCU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in CCU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "CSRU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in CSRU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "MICU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in MICU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "SICU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in SICU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "TSICU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in TSICU", breaks = 25, ylim = c(0,140))
## Clean tmp variable
rm(tmp)
Data Standardization and Factoring (for modeling)
Continuous data, Age and SOFA, will be standardized in the form:
\[z = \frac{x - \mu}{\sigma}\]
## Subset provider-level data to temp variable
temp <- prov
## Create vectors for holding data
lev <- vector()
mod <- vector()
icc <- vector()
std <- vector()
t_r2m <- vector() ## Marginal R_glmm^2
t_R2c <- vector() ## Conditional R^2
d_r2m <- vector() ## Marginal R_glmm^2
d_R2c <- vector() ## Conditional R^2
## Standardize
temp$AGE <- (temp$AGE - mean(temp$AGE))/sd(temp$AGE)
temp$SOFA <- (temp$SOFA - mean(temp$SOFA)/sd(temp$SOFA))Baseline Model (Provider Random Effects)
gm_initial <- glmer(NQF ~ (1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_initial)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 13715.3 13730.0 -6855.7 13711.3 11157
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2922 -0.8253 0.3867 0.7561 3.6611
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.138 1.462
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.16423 0.07617 2.156 0.0311 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sjstats::icc(gm_initial)##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: NQF ~ (1 | CGID)
##
## ICC (CGID): 0.3938## Results
lev <- c(lev, "Provider")
mod <- c(mod, "Model 0")
icc <- c(icc, as.numeric(sjstats::icc(gm_initial)))
std <- c(std, as.numeric(attr(summary(gm_initial)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_initial)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_initial)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_initial)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_initial)[2,2]) ## Conditional R^2
## Check assumption of normality with Q-Q Plot:
plot_model(gm_initial, type = "diag")## $CGID
Model 1a (Patient Demographics)
- Gender
- Age
- Ethnicity
- Language
- Marrital Status
- Care Provider (Random)
gm_a <- glmer(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_a)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 13472.8 13546.0 -6726.4 13452.8 11149
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2322 -0.7905 0.3480 0.7302 4.1883
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.113 1.454
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.22601 0.11618 1.945 0.05172 .
## GENDERM -0.33574 0.04934 -6.804 1.02e-11 ***
## AGE 0.25555 0.02352 10.867 < 2e-16 ***
## ETHNICITYOTHER 0.07824 0.11215 0.698 0.48541
## ETHNICITYWHITE -0.02203 0.08118 -0.271 0.78612
## LANGUAGEOTHER -0.08732 0.07153 -1.221 0.22214
## MARITAL_STATUSSINGLE 0.39320 0.06327 6.215 5.15e-10 ***
## MARITAL_STATUSUNKNOWN -0.14743 0.12465 -1.183 0.23690
## MARITAL_STATUSWIDOWED 0.17955 0.05818 3.086 0.00203 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GENDER AGE ETHNICITYO ETHNICITYW LANGUA
## GENDERM -0.266
## AGE 0.091 -0.054
## ETHNICITYOT -0.465 -0.039 -0.033
## ETHNICITYWH -0.645 -0.064 -0.062 0.654
## LANGUAGEOTH -0.036 0.012 -0.057 -0.282 -0.022
## MARITAL_STATUSS -0.300 0.247 0.023 0.102 0.096 -0.048
## MARITAL_STATUSU -0.130 0.066 -0.078 -0.072 0.047 0.025
## MARITAL_STATUSW -0.369 0.367 -0.176 0.097 0.096 0.000
## MARITAL_STATUSS MARITAL_STATUSU
## GENDERM
## AGE
## ETHNICITYOT
## ETHNICITYWH
## LANGUAGEOTH
## MARITAL_STATUSS
## MARITAL_STATUSU 0.159
## MARITAL_STATUSW 0.431 0.205tab_model(gm_a)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.25 | 1.00 – 1.57 | 0.052 |
| GENDERM | 0.71 | 0.65 – 0.79 | <0.001 |
| AGE | 1.29 | 1.23 – 1.35 | <0.001 |
| ETHNICITYOTHER | 1.08 | 0.87 – 1.35 | 0.485 |
| ETHNICITYWHITE | 0.98 | 0.83 – 1.15 | 0.786 |
| LANGUAGEOTHER | 0.92 | 0.80 – 1.05 | 0.222 |
| MARITAL_STATUSSINGLE | 1.48 | 1.31 – 1.68 | <0.001 |
| MARITAL_STATUSUNKNOWN | 0.86 | 0.68 – 1.10 | 0.237 |
| MARITAL_STATUSWIDOWED | 1.20 | 1.07 – 1.34 | 0.002 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 2.11 | ||
| ICC CGID | 0.39 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.024 / 0.406 | ||
MuMIn::r.squaredGLMM(gm_a)## R2m R2c
## theoretical 0.02444184 0.4060142
## delta 0.02161844 0.3591135## Results
lev <- c(lev, "Provider")
mod <- c(mod, "Model 1a")
icc <- c(icc, as.numeric(sjstats::icc(gm_a)))
std <- c(std, as.numeric(attr(summary(gm_a)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_a)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_a)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_a)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_a)[2,2]) ## Conditional R^2Model 1b (Clinical Variables)
- SOFA
- Care Provider (Random)
gm_b <- glmer(NQF ~ SOFA +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_b)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ SOFA + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 13717.3 13739.3 -6855.7 13711.3 11156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2906 -0.8251 0.3866 0.7564 3.6627
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.138 1.462
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1633086 0.0801375 2.038 0.0416 *
## SOFA 0.0002795 0.0075712 0.037 0.9706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SOFA -0.311tab_model(gm_b)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.18 | 1.01 – 1.38 | 0.042 |
| SOFA | 1.00 | 0.99 – 1.02 | 0.971 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 2.14 | ||
| ICC CGID | 0.39 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.394 | ||
## Results
lev <- c(lev, "Provider")
mod <- c(mod, "Model 1b")
icc <- c(icc, as.numeric(sjstats::icc(gm_b)))
std <- c(std, as.numeric(attr(summary(gm_b)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_b)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_b)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_b)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_b)[2,2]) ## Conditional R^2Model 2a (Demographics + ICU)
- Demographics (Model 1a)
- First Intensive Care Unit
- Care Provider (Random)
gm_2a <- glmer(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
FIRST_CAREUNIT +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_2a)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## FIRST_CAREUNIT + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 13425.3 13527.8 -6698.7 13397.3 11145
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1752 -0.7884 0.3457 0.7332 4.8049
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.181 1.477
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.085238 0.136196 -0.626 0.531412
## GENDERM -0.317860 0.049575 -6.412 1.44e-10 ***
## AGE 0.261045 0.023680 11.024 < 2e-16 ***
## ETHNICITYOTHER 0.126978 0.112951 1.124 0.260936
## ETHNICITYWHITE 0.009285 0.081506 0.114 0.909307
## LANGUAGEOTHER -0.094923 0.071804 -1.322 0.186175
## MARITAL_STATUSSINGLE 0.406708 0.063452 6.410 1.46e-10 ***
## MARITAL_STATUSUNKNOWN -0.193297 0.125557 -1.540 0.123679
## MARITAL_STATUSWIDOWED 0.199316 0.058500 3.407 0.000657 ***
## FIRST_CAREUNITCSRU -0.455543 0.171368 -2.658 0.007854 **
## FIRST_CAREUNITMICU 0.304684 0.078879 3.863 0.000112 ***
## FIRST_CAREUNITSICU 0.343549 0.116794 2.941 0.003266 **
## FIRST_CAREUNITTSICU 0.762857 0.142810 5.342 9.20e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need ittab_model(gm_2a)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.92 | 0.70 – 1.20 | 0.531 |
| GENDERM | 0.73 | 0.66 – 0.80 | <0.001 |
| AGE | 1.30 | 1.24 – 1.36 | <0.001 |
| ETHNICITYOTHER | 1.14 | 0.91 – 1.42 | 0.261 |
| ETHNICITYWHITE | 1.01 | 0.86 – 1.18 | 0.909 |
| LANGUAGEOTHER | 0.91 | 0.79 – 1.05 | 0.186 |
| MARITAL_STATUSSINGLE | 1.50 | 1.33 – 1.70 | <0.001 |
| MARITAL_STATUSUNKNOWN | 0.82 | 0.64 – 1.05 | 0.124 |
| MARITAL_STATUSWIDOWED | 1.22 | 1.09 – 1.37 | 0.001 |
| FIRST_CAREUNITCSRU | 0.63 | 0.45 – 0.89 | 0.008 |
| FIRST_CAREUNITMICU | 1.36 | 1.16 – 1.58 | <0.001 |
| FIRST_CAREUNITSICU | 1.41 | 1.12 – 1.77 | 0.003 |
| FIRST_CAREUNITTSICU | 2.14 | 1.62 – 2.84 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 2.18 | ||
| ICC CGID | 0.40 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.032 / 0.418 | ||
## Results
lev <- c(lev, "Provider")
mod <- c(mod, "Model 2a")
icc <- c(icc, as.numeric(sjstats::icc(gm_2a)))
std <- c(std, as.numeric(attr(summary(gm_2a)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_2a)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_2a)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_2a)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_2a)[2,2]) ## Conditional R^2Model 2b (Clinical Variables + ICU)
- Clinical Variables (Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_2b <- glmer(NQF ~ SOFA +
FIRST_CAREUNIT +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_2b)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ SOFA + FIRST_CAREUNIT + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 13673.1 13724.3 -6829.5 13659.1 11152
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3161 -0.8161 0.3780 0.7567 4.2014
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.243 1.498
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.077708 0.103960 -0.747 0.454771
## SOFA 0.003185 0.007625 0.418 0.676141
## FIRST_CAREUNITCSRU -0.466875 0.168642 -2.768 0.005633 **
## FIRST_CAREUNITMICU 0.241646 0.077318 3.125 0.001776 **
## FIRST_CAREUNITSICU 0.382926 0.115541 3.314 0.000919 ***
## FIRST_CAREUNITTSICU 0.734990 0.140761 5.222 1.77e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) SOFA FIRST_CAREUNITC FIRST_CAREUNITM
## SOFA -0.245
## FIRST_CAREUNITC -0.270 -0.020
## FIRST_CAREUNITM -0.570 -0.018 0.302
## FIRST_CAREUNITS -0.474 0.051 0.290 0.536
## FIRST_CAREUNITT -0.421 0.047 0.268 0.438
## FIRST_CAREUNITS
## SOFA
## FIRST_CAREUNITC
## FIRST_CAREUNITM
## FIRST_CAREUNITS
## FIRST_CAREUNITT 0.557tab_model(gm_2b)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.93 | 0.75 – 1.13 | 0.455 |
| SOFA | 1.00 | 0.99 – 1.02 | 0.676 |
| FIRST_CAREUNITCSRU | 0.63 | 0.45 – 0.87 | 0.006 |
| FIRST_CAREUNITMICU | 1.27 | 1.09 – 1.48 | 0.002 |
| FIRST_CAREUNITSICU | 1.47 | 1.17 – 1.84 | 0.001 |
| FIRST_CAREUNITTSICU | 2.09 | 1.58 – 2.75 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 2.24 | ||
| ICC CGID | 0.41 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.008 / 0.410 | ||
## Results
lev <- c(lev, "Provider")
mod <- c(mod, "Model 2b")
icc <- c(icc, as.numeric(sjstats::icc(gm_2b)))
std <- c(std, as.numeric(attr(summary(gm_2b)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_2b)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_2b)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_2b)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_2b)[2,2]) ## Conditional R^2Model 3 (Demographics + Clinical Variables + ICU)
- Demographics (as in Model 1a)
- SOFA Score (As Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_icu <- glmer(NQF ~
GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
SOFA +
FIRST_CAREUNIT +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_icu)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## SOFA + FIRST_CAREUNIT + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 13420.6 13530.4 -6695.3 13390.6 11144
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0518 -0.7934 0.3424 0.7329 5.0031
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.191 1.48
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.154490 0.138983 -1.112 0.266319
## GENDERM -0.332127 0.049908 -6.655 2.84e-11 ***
## AGE 0.265132 0.023754 11.162 < 2e-16 ***
## ETHNICITYOTHER 0.124913 0.113119 1.104 0.269480
## ETHNICITYWHITE 0.016891 0.081662 0.207 0.836133
## LANGUAGEOTHER -0.099402 0.071854 -1.383 0.166548
## MARITAL_STATUSSINGLE 0.412472 0.063514 6.494 8.35e-11 ***
## MARITAL_STATUSUNKNOWN -0.183439 0.125660 -1.460 0.144344
## MARITAL_STATUSWIDOWED 0.203645 0.058553 3.478 0.000505 ***
## SOFA 0.020376 0.007853 2.595 0.009463 **
## FIRST_CAREUNITCSRU -0.464589 0.171368 -2.711 0.006707 **
## FIRST_CAREUNITMICU 0.301162 0.078934 3.815 0.000136 ***
## FIRST_CAREUNITSICU 0.358239 0.117041 3.061 0.002207 **
## FIRST_CAREUNITTSICU 0.780900 0.143031 5.460 4.77e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need ittab_model(gm_icu)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.86 | 0.65 – 1.13 | 0.266 |
| GENDERM | 0.72 | 0.65 – 0.79 | <0.001 |
| AGE | 1.30 | 1.24 – 1.37 | <0.001 |
| ETHNICITYOTHER | 1.13 | 0.91 – 1.41 | 0.269 |
| ETHNICITYWHITE | 1.02 | 0.87 – 1.19 | 0.836 |
| LANGUAGEOTHER | 0.91 | 0.79 – 1.04 | 0.167 |
| MARITAL_STATUSSINGLE | 1.51 | 1.33 – 1.71 | <0.001 |
| MARITAL_STATUSUNKNOWN | 0.83 | 0.65 – 1.06 | 0.144 |
| MARITAL_STATUSWIDOWED | 1.23 | 1.09 – 1.37 | 0.001 |
| SOFA | 1.02 | 1.00 – 1.04 | 0.009 |
| FIRST_CAREUNITCSRU | 0.63 | 0.45 – 0.88 | 0.007 |
| FIRST_CAREUNITMICU | 1.35 | 1.16 – 1.58 | <0.001 |
| FIRST_CAREUNITSICU | 1.43 | 1.14 – 1.80 | 0.002 |
| FIRST_CAREUNITTSICU | 2.18 | 1.65 – 2.89 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 2.19 | ||
| ICC CGID | 0.40 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.033 / 0.419 | ||
## Results
lev <- c(lev, "Provider")
mod <- c(mod, "Model 3")
icc <- c(icc, as.numeric(sjstats::icc(gm_icu)))
std <- c(std, as.numeric(attr(summary(gm_icu)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_icu)[1,1]) ## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_icu)[1,2]) ## Conditional R^2## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_icu)[2,1]) ## Marginal R_glmm^2 ## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_icu)[2,2]) ## Conditional R^2## Warning: The null model is correct only if all variables used by the
## original model remain unchanged.Results
results <- as.data.frame(cbind(mod, icc, std, lev, t_r2m, t_R2c, d_r2m, d_R2c))
## Correct variables
results$mod <- as.character(results$mod)
results$icc <- as.numeric(as.character(results$icc))
results$std <- as.numeric(as.character(results$std))
results$lev <- as.character(results$lev)
print(results)## mod icc std lev t_r2m
## 1 Model 0 0.3938457 1.462045 Provider 0
## 2 Model 1a 0.3911323 1.453750 Provider 0.0244418379564992
## 3 Model 1b 0.3938446 1.462042 Provider 1.29522660237321e-07
## 4 Model 2a 0.3986534 1.476810 Provider 0.0322298072241591
## 5 Model 2b 0.4053823 1.497624 Provider 0.0076058384804099
## 6 Model 3 0.3997552 1.480206 Provider 0.0327984422997686
## t_R2c d_r2m d_R2c
## 1 0.393845690742604 0 0.34752823595739
## 2 0.406014175521916 0.0216184415541336 0.359113489717591
## 3 0.393844701204199 1.14290378300496e-07 0.347527296071577
## 4 0.418034637784055 0.0285735799337614 0.370611777313085
## 5 0.409904842781096 0.0067323447518925 0.362829255995969
## 6 0.419442255159451 0.0290856883952417 0.37196177250928Mixed Model ROC Curves (Provider Level)
par(mfrow=c(2,3))
rocplot(as.numeric(predict(gm_initial, type="response")), temp$NQF, col="blue", main = "Initial Model")
rocplot(as.numeric(predict(gm_a, type="response")), temp$NQF, col="blue", main = "Model 1a (Demographics)")
rocplot(as.numeric(predict(gm_b, type="response")), temp$NQF, col="blue", main = "Model 1b (Clinical Variables)")
rocplot(as.numeric(predict(gm_2a, type="response")), temp$NQF, col="blue", main = "Model 2a (Demographics + ICU)")
rocplot(as.numeric(predict(gm_2b, type="response")), temp$NQF, col="blue", main = "Model 2b (Clinical Variables + ICU)")
rocplot(as.numeric(predict(gm_icu, type="response")), temp$NQF, col="blue", main = "Model 3 (All Covariates)")
Generalized Logistic Regression Models (Provider Level)
Baseline Model (Provider Random Effects)
Not meaningful.
Model 1a (Patient Demographics)
- Gender
- Age
- Ethnicity
- Language
- Marrital Status
- Care Provider (Random)
gm_a <- glm(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_a)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS,
## family = binomial(link = "logit"), data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7122 -1.2237 0.8905 1.0774 1.5698
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.42029 0.07696 5.461 4.73e-08 ***
## GENDERM -0.31292 0.04191 -7.466 8.24e-14 ***
## AGE 0.25175 0.02012 12.515 < 2e-16 ***
## ETHNICITYOTHER 0.02276 0.09431 0.241 0.809266
## ETHNICITYWHITE -0.16037 0.06910 -2.321 0.020292 *
## LANGUAGEOTHER -0.06968 0.06123 -1.138 0.255150
## MARITAL_STATUSSINGLE 0.40719 0.05425 7.506 6.09e-14 ***
## MARITAL_STATUSUNKNOWN -0.39961 0.10615 -3.764 0.000167 ***
## MARITAL_STATUSWIDOWED 0.10637 0.04918 2.163 0.030538 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 15330 on 11158 degrees of freedom
## Residual deviance: 14979 on 11150 degrees of freedom
## AIC: 14997
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 1.5224080 1.3097351 1.7710370
## GENDERM 0.7313053 0.6736128 0.7938953
## AGE 1.2862732 1.2366173 1.3380790
## ETHNICITYOTHER 1.0230253 0.8503155 1.2307284
## ETHNICITYWHITE 0.8518327 0.7435769 0.9749569
## LANGUAGEOTHER 0.9326962 0.8272891 1.0517552
## MARITAL_STATUSSINGLE 1.5025887 1.3512257 1.6714201
## MARITAL_STATUSUNKNOWN 0.6705842 0.5442147 0.8252715
## MARITAL_STATUSWIDOWED 1.1122330 1.0100302 1.2247760Model 1b (Clinical Variables)
- SOFA
- Care Provider (Random)
gm_b <- glm(NQF ~ SOFA,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_b)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ SOFA, family = binomial(link = "logit"),
## data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.276 -1.274 1.083 1.084 1.087
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2272396 0.0285897 7.948 1.89e-15 ***
## SOFA -0.0007278 0.0063508 -0.115 0.909
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 15330 on 11158 degrees of freedom
## Residual deviance: 15330 on 11157 degrees of freedom
## AIC: 15334
##
## Number of Fisher Scoring iterations: 3
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 1.2551306 1.1867773 1.327527
## SOFA 0.9992725 0.9869175 1.011798Model 2a (Demographics + ICU)
- Demographics (Model 1a)
- First Intensive Care Unit
- Care Provider (Random)
gm_2a <- glm(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
FIRST_CAREUNIT,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_2a)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## FIRST_CAREUNIT, family = binomial(link = "logit"), data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7363 -1.2211 0.8301 1.0551 1.7691
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.36799 0.09278 3.966 7.30e-05 ***
## GENDERM -0.29393 0.04241 -6.931 4.18e-12 ***
## AGE 0.25718 0.02029 12.672 < 2e-16 ***
## ETHNICITYOTHER 0.12885 0.09586 1.344 0.17890
## ETHNICITYWHITE -0.13000 0.06959 -1.868 0.06173 .
## LANGUAGEOTHER -0.13815 0.06191 -2.231 0.02565 *
## MARITAL_STATUSSINGLE 0.39887 0.05455 7.312 2.64e-13 ***
## MARITAL_STATUSUNKNOWN -0.29833 0.10814 -2.759 0.00580 **
## MARITAL_STATUSWIDOWED 0.15096 0.04983 3.029 0.00245 **
## FIRST_CAREUNITCSRU -1.08710 0.12548 -8.664 < 2e-16 ***
## FIRST_CAREUNITMICU 0.13879 0.05708 2.432 0.01503 *
## FIRST_CAREUNITSICU -0.36409 0.07958 -4.575 4.76e-06 ***
## FIRST_CAREUNITTSICU -0.28484 0.09373 -3.039 0.00237 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 15330 on 11158 degrees of freedom
## Residual deviance: 14799 on 11146 degrees of freedom
## AIC: 14825
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 1.4448339 1.2050182 1.7336817
## GENDERM 0.7453311 0.6858651 0.8099145
## AGE 1.2932805 1.2429184 1.3458423
## ETHNICITYOTHER 1.1375187 0.9426591 1.3726823
## ETHNICITYWHITE 0.8780936 0.7657631 1.0059752
## LANGUAGEOTHER 0.8709679 0.7714980 0.9834440
## MARITAL_STATUSSINGLE 1.4901376 1.3392293 1.6585714
## MARITAL_STATUSUNKNOWN 0.7420538 0.5999307 0.9168668
## MARITAL_STATUSWIDOWED 1.1629461 1.0547389 1.2822969
## FIRST_CAREUNITCSRU 0.3371922 0.2628892 0.4300811
## FIRST_CAREUNITMICU 1.1488823 1.0271319 1.2847139
## FIRST_CAREUNITSICU 0.6948277 0.5943764 0.8120136
## FIRST_CAREUNITTSICU 0.7521322 0.6257859 0.9037148Model 2b (Clinical Variables + ICU)
- Clinical Variables (Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_2b <- glm(NQF ~ SOFA +
FIRST_CAREUNIT,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_2b)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ SOFA + FIRST_CAREUNIT, family = binomial(link = "logit"),
## data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.338 -1.323 1.027 1.036 1.579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.273813 0.054575 5.017 5.24e-07 ***
## SOFA -0.003461 0.006437 -0.538 0.591
## FIRST_CAREUNITCSRU -1.159929 0.123438 -9.397 < 2e-16 ***
## FIRST_CAREUNITMICU 0.090876 0.055739 1.630 0.103
## FIRST_CAREUNITSICU -0.346994 0.077920 -4.453 8.46e-06 ***
## FIRST_CAREUNITTSICU -0.357304 0.091384 -3.910 9.23e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 15330 on 11158 degrees of freedom
## Residual deviance: 15143 on 11153 degrees of freedom
## AIC: 15155
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 1.3149687 1.1817847 1.4637357
## SOFA 0.9965445 0.9840557 1.0092069
## FIRST_CAREUNITCSRU 0.3135084 0.2453797 0.3982432
## FIRST_CAREUNITMICU 1.0951331 0.9816228 1.2213724
## FIRST_CAREUNITSICU 0.7068092 0.6066043 0.8233301
## FIRST_CAREUNITTSICU 0.6995600 0.5847083 0.8366620Model 3 (Demographics + Clinical Variables + ICU)
- Demographics (as in Model 1a)
- SOFA Score (As Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_icu <- glm(NQF ~
GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
SOFA +
FIRST_CAREUNIT,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_icu)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## SOFA + FIRST_CAREUNIT, family = binomial(link = "logit"),
## data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7471 -1.2200 0.8355 1.0598 1.7589
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.328197 0.095340 3.442 0.000577 ***
## GENDERM -0.302979 0.042710 -7.094 1.30e-12 ***
## AGE 0.260069 0.020365 12.770 < 2e-16 ***
## ETHNICITYOTHER 0.128288 0.095918 1.337 0.181070
## ETHNICITYWHITE -0.124975 0.069683 -1.793 0.072896 .
## LANGUAGEOTHER -0.142354 0.061977 -2.297 0.021626 *
## MARITAL_STATUSSINGLE 0.401902 0.054588 7.362 1.81e-13 ***
## MARITAL_STATUSUNKNOWN -0.293350 0.108211 -2.711 0.006710 **
## MARITAL_STATUSWIDOWED 0.153202 0.049852 3.073 0.002118 **
## SOFA 0.012153 0.006663 1.824 0.068142 .
## FIRST_CAREUNITCSRU -1.095251 0.125564 -8.723 < 2e-16 ***
## FIRST_CAREUNITMICU 0.135679 0.057109 2.376 0.017511 *
## FIRST_CAREUNITSICU -0.358829 0.079653 -4.505 6.64e-06 ***
## FIRST_CAREUNITTSICU -0.277514 0.093789 -2.959 0.003087 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 15330 on 11158 degrees of freedom
## Residual deviance: 14796 on 11145 degrees of freedom
## AIC: 14824
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 1.3884625 1.1521673 1.6743439
## GENDERM 0.7386145 0.6792807 0.8030862
## AGE 1.2970196 1.2463410 1.3499227
## ETHNICITYOTHER 1.1368801 0.9420207 1.3720761
## ETHNICITYWHITE 0.8825186 0.7694794 1.0112434
## LANGUAGEOTHER 0.8673141 0.7681590 0.9794451
## MARITAL_STATUSSINGLE 1.4946647 1.3432087 1.6637255
## MARITAL_STATUSUNKNOWN 0.7457614 0.6028470 0.9215801
## MARITAL_STATUSWIDOWED 1.1655600 1.0570741 1.2852229
## SOFA 1.0122273 0.9991059 1.0255465
## FIRST_CAREUNITCSRU 0.3344556 0.2607124 0.4266638
## FIRST_CAREUNITMICU 1.1453145 1.0238761 1.2808048
## FIRST_CAREUNITSICU 0.6984937 0.5974332 0.8164097
## FIRST_CAREUNITTSICU 0.7576650 0.6303177 0.9104643Logistic Regression ROC Curves (Provider Level)
par(mfrow=c(2,3))
rocplot(as.numeric(predict(gm_a, type="response")), temp$NQF, col="blue", main = "Model 1a (Demographics)")
rocplot(as.numeric(predict(gm_b, type="response")), temp$NQF, col="blue", main = "Model 1b (Clinical Variables)")
rocplot(as.numeric(predict(gm_2a, type="response")), temp$NQF, col="blue", main = "Model 2a (Demographics + ICU)")
rocplot(as.numeric(predict(gm_2b, type="response")), temp$NQF, col="blue", main = "Model 2b (Clinical Variables + ICU)")
rocplot(as.numeric(predict(gm_icu, type="response")), temp$NQF, col="blue", main = "Model 3 (All Covariates)")
Admission-Level Care Team Analysis
For this analysis, we will look at the clinical care team. If any care provider documented the patient care preferences during the admission, then the entire care team will get credit for the documentation.
Initial Rates of Documentation Check
Here we will aggregate the care providers and ICUs to average rates of documentation as part of care teams.
tmp <- aggregate(NQF ~
CGID + FIRST_CAREUNIT,
data = team,
FUN = mean)
par(mfrow=c(2,3))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "CCU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in CCU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "CSRU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in CSRU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "MICU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in MICU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "SICU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in SICU", breaks = 25, ylim = c(0,140))
hist(tmp$NQF[tmp$FIRST_CAREUNIT == "TSICU"],
xlab = "Rates of Caremeasure Documentation",
main = "Rates of Caremeasure Documentation\n in TSICU", breaks = 25, ylim = c(0,140))
## Clean tmp variable
rm(tmp)
Continuous Data Standardization
Continuous data, Age and SOFA, will be standardized in the form:
\[z = \frac{x - \mu}{\sigma}\]
temp <- team
lev <- vector()
mod <- vector()
icc <- vector()
std <- vector()
t_r2m <- vector() ## Marginal R_glmm^2
t_R2c <- vector() ## Conditional R^2
d_r2m <- vector() ## Marginal R_glmm^2
d_R2c <- vector() ## Conditional R^2
temp$AGE <- (temp$AGE - mean(temp$AGE))/sd(temp$AGE)
temp$SOFA <- (temp$SOFA - mean(temp$SOFA)/sd(temp$SOFA))Baseline Model (Team Random Effects)
gm_initial <- glmer(NQF ~ (1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_initial)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 11907.8 11922.4 -5951.9 11903.8 11157
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6384 -0.7102 0.4009 0.5699 2.2754
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.011 1.418
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.05472 0.07638 13.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sjstats::icc(gm_initial)##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: NQF ~ (1 | CGID)
##
## ICC (CGID): 0.3794## Results
lev <- c(lev, "Team")
mod <- c(mod, "Model 0")
icc <- c(icc, as.numeric(sjstats::icc(gm_initial)))
std <- c(std, as.numeric(attr(summary(gm_initial)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_initial)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_initial)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_initial)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_initial)[2,2]) ## Conditional R^2
## Check assumption of normality:
plot_model(gm_initial, type = "diag")## $CGID
Model 1a (Patient Demographics)
- Gender
- Age
- Ethnicity
- Language
- Marrital Status
- Care Provider (Random)
gm_a <- glmer(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_a)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 11664.4 11737.6 -5822.2 11644.4 11149
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.3957 -0.6490 0.3757 0.5826 2.7628
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 1.938 1.392
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.07216 0.12272 8.737 < 2e-16 ***
## GENDERM -0.36439 0.05429 -6.712 1.92e-11 ***
## AGE 0.25466 0.02611 9.754 < 2e-16 ***
## ETHNICITYOTHER 0.52036 0.12715 4.092 4.27e-05 ***
## ETHNICITYWHITE 0.01858 0.09069 0.205 0.8377
## LANGUAGEOTHER -0.09035 0.07904 -1.143 0.2530
## MARITAL_STATUSSINGLE 0.50454 0.07127 7.079 1.45e-12 ***
## MARITAL_STATUSUNKNOWN -0.38865 0.12853 -3.024 0.0025 **
## MARITAL_STATUSWIDOWED 0.14116 0.06339 2.227 0.0260 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GENDER AGE ETHNICITYO ETHNICITYW LANGUA
## GENDERM -0.275
## AGE 0.111 -0.061
## ETHNICITYOT -0.479 -0.050 -0.016
## ETHNICITYWH -0.675 -0.078 -0.067 0.649
## LANGUAGEOTH -0.049 0.003 -0.055 -0.243 -0.011
## MARITAL_STATUSS -0.287 0.226 0.035 0.097 0.087 -0.032
## MARITAL_STATUSU -0.134 0.068 -0.079 -0.086 0.053 0.028
## MARITAL_STATUSW -0.372 0.361 -0.175 0.098 0.095 0.001
## MARITAL_STATUSS MARITAL_STATUSU
## GENDERM
## AGE
## ETHNICITYOT
## ETHNICITYWH
## LANGUAGEOTH
## MARITAL_STATUSS
## MARITAL_STATUSU 0.149
## MARITAL_STATUSW 0.390 0.202tab_model(gm_a)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.92 | 2.30 – 3.72 | <0.001 |
| GENDERM | 0.69 | 0.62 – 0.77 | <0.001 |
| AGE | 1.29 | 1.23 – 1.36 | <0.001 |
| ETHNICITYOTHER | 1.68 | 1.31 – 2.16 | <0.001 |
| ETHNICITYWHITE | 1.02 | 0.85 – 1.22 | 0.838 |
| LANGUAGEOTHER | 0.91 | 0.78 – 1.07 | 0.253 |
| MARITAL_STATUSSINGLE | 1.66 | 1.44 – 1.90 | <0.001 |
| MARITAL_STATUSUNKNOWN | 0.68 | 0.53 – 0.87 | 0.002 |
| MARITAL_STATUSWIDOWED | 1.15 | 1.02 – 1.30 | 0.026 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 1.94 | ||
| ICC CGID | 0.37 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.032 / 0.391 | ||
MuMIn::r.squaredGLMM(gm_a)## R2m R2c
## theoretical 0.03168954 0.3905952
## delta 0.02549255 0.3142131## Results
lev <- c(lev, "Team")
mod <- c(mod, "Model 1a")
icc <- c(icc, as.numeric(sjstats::icc(gm_a)))
std <- c(std, as.numeric(attr(summary(gm_a)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_a)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_a)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_a)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_a)[2,2]) ## Conditional R^2Model 1b (Clinical Variables)
- SOFA
- Care Provider (Random)
gm_b <- glmer(NQF ~ SOFA +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_b)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ SOFA + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 11898.9 11920.9 -5946.5 11892.9 11156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4794 -0.7124 0.3938 0.5783 2.4456
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 2.006 1.416
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.964591 0.080949 11.92 < 2e-16 ***
## SOFA 0.027701 0.008446 3.28 0.00104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## SOFA -0.334tab_model(gm_b)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.62 | 2.24 – 3.07 | <0.001 |
| SOFA | 1.03 | 1.01 – 1.05 | 0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 2.01 | ||
| ICC CGID | 0.38 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.001 / 0.380 | ||
## Results
lev <- c(lev, "Team")
mod <- c(mod, "Model 1b")
icc <- c(icc, as.numeric(sjstats::icc(gm_b)))
std <- c(std, as.numeric(attr(summary(gm_b)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_b)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_b)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_b)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_b)[2,2]) ## Conditional R^2Model 2a (Demographics + ICU)
- Demographics (Model 1a)
- First Intensive Care Unit
- Care Provider (Random)
gm_2a <- glmer(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
FIRST_CAREUNIT +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_2a)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## FIRST_CAREUNIT + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 11620.3 11722.7 -5796.1 11592.3 11145
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4094 -0.6266 0.3749 0.5797 3.1576
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 1.796 1.34
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.98383 0.14257 6.901 5.17e-12 ***
## GENDERM -0.34563 0.05450 -6.342 2.27e-10 ***
## AGE 0.25973 0.02619 9.915 < 2e-16 ***
## ETHNICITYOTHER 0.57812 0.12772 4.526 6.00e-06 ***
## ETHNICITYWHITE 0.03539 0.09058 0.391 0.69604
## LANGUAGEOTHER -0.11338 0.07922 -1.431 0.15239
## MARITAL_STATUSSINGLE 0.52265 0.07143 7.317 2.54e-13 ***
## MARITAL_STATUSUNKNOWN -0.40570 0.12905 -3.144 0.00167 **
## MARITAL_STATUSWIDOWED 0.17427 0.06369 2.736 0.00622 **
## FIRST_CAREUNITCSRU -0.76801 0.17089 -4.494 6.99e-06 ***
## FIRST_CAREUNITMICU 0.13875 0.08690 1.597 0.11033
## FIRST_CAREUNITSICU -0.21347 0.12372 -1.725 0.08444 .
## FIRST_CAREUNITTSICU 0.32692 0.14829 2.205 0.02748 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need ittab_model(gm_2a)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.67 | 2.02 – 3.54 | <0.001 |
| GENDERM | 0.71 | 0.64 – 0.79 | <0.001 |
| AGE | 1.30 | 1.23 – 1.36 | <0.001 |
| ETHNICITYOTHER | 1.78 | 1.39 – 2.29 | <0.001 |
| ETHNICITYWHITE | 1.04 | 0.87 – 1.24 | 0.696 |
| LANGUAGEOTHER | 0.89 | 0.76 – 1.04 | 0.152 |
| MARITAL_STATUSSINGLE | 1.69 | 1.47 – 1.94 | <0.001 |
| MARITAL_STATUSUNKNOWN | 0.67 | 0.52 – 0.86 | 0.002 |
| MARITAL_STATUSWIDOWED | 1.19 | 1.05 – 1.35 | 0.006 |
| FIRST_CAREUNITCSRU | 0.46 | 0.33 – 0.65 | <0.001 |
| FIRST_CAREUNITMICU | 1.15 | 0.97 – 1.36 | 0.110 |
| FIRST_CAREUNITSICU | 0.81 | 0.63 – 1.03 | 0.084 |
| FIRST_CAREUNITTSICU | 1.39 | 1.04 – 1.85 | 0.027 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 1.80 | ||
| ICC CGID | 0.35 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.041 / 0.379 | ||
## Results
lev <- c(lev, "Team")
mod <- c(mod, "Model 2a")
icc <- c(icc, as.numeric(sjstats::icc(gm_2a)))
std <- c(std, as.numeric(attr(summary(gm_2a)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_2a)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_2a)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_2a)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_2a)[2,2]) ## Conditional R^2Model 2b (Clinical Variables + ICU)
- Clinical Variables (Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_2b <- glmer(NQF ~ SOFA +
FIRST_CAREUNIT +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_2b)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ SOFA + FIRST_CAREUNIT + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 11861.9 11913.1 -5923.9 11847.9 11152
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.5217 -0.6889 0.3914 0.5803 2.7499
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 1.913 1.383
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.947851 0.105952 8.946 < 2e-16 ***
## SOFA 0.028912 0.008497 3.403 0.000667 ***
## FIRST_CAREUNITCSRU -0.800718 0.168384 -4.755 1.98e-06 ***
## FIRST_CAREUNITMICU 0.068061 0.085470 0.796 0.425849
## FIRST_CAREUNITSICU -0.136505 0.122717 -1.112 0.265985
## FIRST_CAREUNITTSICU 0.308234 0.145968 2.112 0.034716 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) SOFA FIRST_CAREUNITC FIRST_CAREUNITM
## SOFA -0.266
## FIRST_CAREUNITC -0.321 -0.017
## FIRST_CAREUNITM -0.623 -0.006 0.358
## FIRST_CAREUNITS -0.501 0.060 0.352 0.557
## FIRST_CAREUNITT -0.438 0.051 0.324 0.467
## FIRST_CAREUNITS
## SOFA
## FIRST_CAREUNITC
## FIRST_CAREUNITM
## FIRST_CAREUNITS
## FIRST_CAREUNITT 0.596tab_model(gm_2b)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.58 | 2.10 – 3.18 | <0.001 |
| SOFA | 1.03 | 1.01 – 1.05 | 0.001 |
| FIRST_CAREUNITCSRU | 0.45 | 0.32 – 0.62 | <0.001 |
| FIRST_CAREUNITMICU | 1.07 | 0.91 – 1.27 | 0.426 |
| FIRST_CAREUNITSICU | 0.87 | 0.69 – 1.11 | 0.266 |
| FIRST_CAREUNITTSICU | 1.36 | 1.02 – 1.81 | 0.035 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 1.91 | ||
| ICC CGID | 0.37 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.007 / 0.372 | ||
## Results
lev <- c(lev, "Team")
mod <- c(mod, "Model 2b")
icc <- c(icc, as.numeric(sjstats::icc(gm_2b)))
std <- c(std, as.numeric(attr(summary(gm_2b)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_2b)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_2b)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_2b)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_2b)[2,2]) ## Conditional R^2Model 3 (Demographics + Clinical Variables + ICU)
- Demographics (as in Model 1a)
- SOFA Score (As Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_icu <- glmer(NQF ~
GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
SOFA +
FIRST_CAREUNIT +
(1 | CGID),
data = temp,
family = binomial,
control = glmerControl(optimizer = "bobyqa"),
nAGQ = 10)
## View
summary(gm_icu)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## SOFA + FIRST_CAREUNIT + (1 | CGID)
## Data: temp
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 11594.9 11704.7 -5782.4 11564.9 11144
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4725 -0.6227 0.3728 0.5760 3.4682
##
## Random effects:
## Groups Name Variance Std.Dev.
## CGID (Intercept) 1.808 1.345
## Number of obs: 11159, groups: CGID, 493
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.830118 0.145771 5.695 1.24e-08 ***
## GENDERM -0.375047 0.054874 -6.835 8.21e-12 ***
## AGE 0.268603 0.026323 10.204 < 2e-16 ***
## ETHNICITYOTHER 0.582375 0.128231 4.542 5.58e-06 ***
## ETHNICITYWHITE 0.049426 0.090794 0.544 0.58618
## LANGUAGEOTHER -0.128255 0.079515 -1.613 0.10675
## MARITAL_STATUSSINGLE 0.535471 0.071580 7.481 7.39e-14 ***
## MARITAL_STATUSUNKNOWN -0.382603 0.129097 -2.964 0.00304 **
## MARITAL_STATUSWIDOWED 0.184334 0.063793 2.890 0.00386 **
## SOFA 0.045565 0.008774 5.193 2.07e-07 ***
## FIRST_CAREUNITCSRU -0.784574 0.171129 -4.585 4.55e-06 ***
## FIRST_CAREUNITMICU 0.135794 0.086988 1.561 0.11851
## FIRST_CAREUNITSICU -0.180827 0.124228 -1.456 0.14550
## FIRST_CAREUNITTSICU 0.361765 0.148627 2.434 0.01493 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need ittab_model(gm_icu)| Â | NQF | ||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.29 | 1.72 – 3.05 | <0.001 |
| GENDERM | 0.69 | 0.62 – 0.77 | <0.001 |
| AGE | 1.31 | 1.24 – 1.38 | <0.001 |
| ETHNICITYOTHER | 1.79 | 1.39 – 2.30 | <0.001 |
| ETHNICITYWHITE | 1.05 | 0.88 – 1.26 | 0.586 |
| LANGUAGEOTHER | 0.88 | 0.75 – 1.03 | 0.107 |
| MARITAL_STATUSSINGLE | 1.71 | 1.48 – 1.97 | <0.001 |
| MARITAL_STATUSUNKNOWN | 0.68 | 0.53 – 0.88 | 0.003 |
| MARITAL_STATUSWIDOWED | 1.20 | 1.06 – 1.36 | 0.004 |
| SOFA | 1.05 | 1.03 – 1.06 | <0.001 |
| FIRST_CAREUNITCSRU | 0.46 | 0.33 – 0.64 | <0.001 |
| FIRST_CAREUNITMICU | 1.15 | 0.97 – 1.36 | 0.119 |
| FIRST_CAREUNITSICU | 0.83 | 0.65 – 1.06 | 0.146 |
| FIRST_CAREUNITTSICU | 1.44 | 1.07 – 1.92 | 0.015 |
| Random Effects | |||
| σ2 | 3.29 | ||
| Ï„00 CGID | 1.81 | ||
| ICC CGID | 0.35 | ||
| Observations | 11159 | ||
| Marginal R2 / Conditional R2 | 0.044 / 0.383 | ||
## Results
lev <- c(lev, "Team")
mod <- c(mod, "Model 3")
icc <- c(icc, as.numeric(sjstats::icc(gm_icu)))
std <- c(std, as.numeric(attr(summary(gm_icu)$varcor$CGID, "stddev")))
t_r2m <- c(t_r2m, MuMIn::r.squaredGLMM(gm_icu)[1,1])
t_R2c <- c(t_R2c, MuMIn::r.squaredGLMM(gm_icu)[1,2]) ## Conditional R^2
d_r2m <- c(d_r2m, MuMIn::r.squaredGLMM(gm_icu)[2,1]) ## Marginal R_glmm^2
d_R2c <- c(d_R2c, MuMIn::r.squaredGLMM(gm_icu)[2,2]) ## Conditional R^2Mixed Model ROC Curves (Team Level)
par(mfrow=c(2,3))
rocplot(as.numeric(predict(gm_initial, type="response")), temp$NQF, col="blue", main = "Initial Model")
rocplot(as.numeric(predict(gm_a, type="response")), temp$NQF, col="blue", main = "Model 1a (Demographics)")
rocplot(as.numeric(predict(gm_b, type="response")), temp$NQF, col="blue", main = "Model 1b (Clinical Variables)")
rocplot(as.numeric(predict(gm_2a, type="response")), temp$NQF, col="blue", main = "Model 2a (Demographics + ICU)")
rocplot(as.numeric(predict(gm_2b, type="response")), temp$NQF, col="blue", main = "Model 2b (Clinical Variables + ICU)")
rocplot(as.numeric(predict(gm_icu, type="response")), temp$NQF, col="blue", main = "Model 3 (All Covariates)")
Generalized Logistic Regression Models (Team Level)
Baseline Model (Provider Random Effects)
Not meaningful.
Model 1a (Patient Demographics)
- Gender
- Age
- Ethnicity
- Languagee
- Marrital Status
- Care Provider (Random)
gm_a <- glm(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_a)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS,
## family = binomial(link = "logit"), data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0450 -1.3494 0.7126 0.8503 1.3895
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.14611 0.08677 13.209 < 2e-16 ***
## GENDERM -0.32538 0.04645 -7.005 2.47e-12 ***
## AGE 0.25832 0.02236 11.554 < 2e-16 ***
## ETHNICITYOTHER 0.28433 0.10866 2.617 0.00888 **
## ETHNICITYWHITE -0.17346 0.07844 -2.211 0.02701 *
## LANGUAGEOTHER -0.03501 0.06831 -0.513 0.60826
## MARITAL_STATUSSINGLE 0.52681 0.06216 8.475 < 2e-16 ***
## MARITAL_STATUSUNKNOWN -0.68721 0.10841 -6.339 2.31e-10 ***
## MARITAL_STATUSWIDOWED 0.05823 0.05395 1.079 0.28049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13281 on 11158 degrees of freedom
## Residual deviance: 12914 on 11150 degrees of freedom
## AIC: 12932
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 3.1459451 2.6572706 3.7341953
## GENDERM 0.7222558 0.6593555 0.7910434
## AGE 1.2947534 1.2393381 1.3528639
## ETHNICITYOTHER 1.3288740 1.0740039 1.6445827
## ETHNICITYWHITE 0.8407523 0.7198356 0.9790896
## LANGUAGEOTHER 0.9655946 0.8452159 1.1047981
## MARITAL_STATUSSINGLE 1.6935260 1.5001033 1.9140858
## MARITAL_STATUSUNKNOWN 0.5029796 0.4068938 0.6225242
## MARITAL_STATUSWIDOWED 1.0599546 0.9536594 1.1782824Model 1b (Clinical Variables)
- SOFA
- Care Provider (Random)
gm_b <- glm(NQF ~ SOFA,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_b)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ SOFA, family = binomial(link = "logit"),
## data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7808 -1.5470 0.8029 0.8254 0.8717
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.826230 0.031302 26.395 < 2e-16 ***
## SOFA 0.032490 0.007181 4.525 6.05e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13281 on 11158 degrees of freedom
## Residual deviance: 13261 on 11157 degrees of freedom
## AIC: 13265
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 2.284689 2.149072 2.429646
## SOFA 1.033024 1.018642 1.047725Model 2a (Demographics + ICU)
- Demographics (Model 1a)
- First Intensive Care Unit
- Care Provider (Random)
gm_2a <- glm(NQF ~ GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
FIRST_CAREUNIT,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_2a)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## FIRST_CAREUNIT, family = binomial(link = "logit"), data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1913 -1.1323 0.6659 0.8157 1.5564
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.18385 0.10583 11.186 < 2e-16 ***
## GENDERM -0.31922 0.04750 -6.720 1.81e-11 ***
## AGE 0.27033 0.02279 11.861 < 2e-16 ***
## ETHNICITYOTHER 0.47616 0.11192 4.254 2.10e-05 ***
## ETHNICITYWHITE -0.13330 0.07985 -1.669 0.0950 .
## LANGUAGEOTHER -0.16040 0.06974 -2.300 0.0214 *
## MARITAL_STATUSSINGLE 0.50884 0.06306 8.069 7.10e-16 ***
## MARITAL_STATUSUNKNOWN -0.50575 0.11215 -4.510 6.49e-06 ***
## MARITAL_STATUSWIDOWED 0.11793 0.05527 2.134 0.0329 *
## FIRST_CAREUNITCSRU -1.19074 0.12032 -9.897 < 2e-16 ***
## FIRST_CAREUNITMICU 0.14864 0.06518 2.281 0.0226 *
## FIRST_CAREUNITSICU -0.84509 0.08508 -9.933 < 2e-16 ***
## FIRST_CAREUNITTSICU -0.66840 0.09909 -6.745 1.53e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13281 on 11158 degrees of freedom
## Residual deviance: 12540 on 11146 degrees of freedom
## AIC: 12566
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 3.2669367 2.6582230 4.0251932
## GENDERM 0.7267132 0.6620572 0.7975655
## AGE 1.3103975 1.2532525 1.3703799
## ETHNICITYOTHER 1.6098768 1.2929677 2.0052989
## ETHNICITYWHITE 0.8752003 0.7472729 1.0220439
## LANGUAGEOTHER 0.8517989 0.7434865 0.9772923
## MARITAL_STATUSSINGLE 1.6633622 1.4707744 1.8833109
## MARITAL_STATUSUNKNOWN 0.6030541 0.4843306 0.7519202
## MARITAL_STATUSWIDOWED 1.1251604 1.0097372 1.2540385
## FIRST_CAREUNITCSRU 0.3039954 0.2399291 0.3846012
## FIRST_CAREUNITMICU 1.1602561 1.0203280 1.3174134
## FIRST_CAREUNITSICU 0.4295194 0.3634068 0.5072978
## FIRST_CAREUNITTSICU 0.5125257 0.4220521 0.6224501Model 2b (Clinical Variables + ICU)
- Clinical Variables (Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_2b <- glm(NQF ~ SOFA +
FIRST_CAREUNIT,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_2b)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ SOFA + FIRST_CAREUNIT, family = binomial(link = "logit"),
## data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8415 -1.2707 0.7305 0.7639 1.3155
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.997956 0.061864 16.131 < 2e-16 ***
## SOFA 0.025667 0.007342 3.496 0.000473 ***
## FIRST_CAREUNITCSRU -1.299375 0.117808 -11.030 < 2e-16 ***
## FIRST_CAREUNITMICU 0.076010 0.063712 1.193 0.232859
## FIRST_CAREUNITSICU -0.800513 0.083025 -9.642 < 2e-16 ***
## FIRST_CAREUNITTSICU -0.738141 0.096241 -7.670 1.72e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13281 on 11158 degrees of freedom
## Residual deviance: 12892 on 11153 degrees of freedom
## AIC: 12904
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 2.7127318 2.4051091 3.0653430
## SOFA 1.0259989 1.0113926 1.0409289
## FIRST_CAREUNITCSRU 0.2727022 0.2162857 0.3433094
## FIRST_CAREUNITMICU 1.0789731 0.9515195 1.2215327
## FIRST_CAREUNITSICU 0.4490984 0.3815144 0.5282975
## FIRST_CAREUNITTSICU 0.4780018 0.3958215 0.5772770Model 3 (Demographics + Clinical Variables + ICU)
- Demographics (as in Model 1a)
- SOFA Score (As Model 1b)
- First Intensive Care Unit
- Care Provider (Random)
gm_icu <- glm(NQF ~
GENDER +
AGE +
ETHNICITY +
LANGUAGE +
MARITAL_STATUS +
SOFA +
FIRST_CAREUNIT,
data = temp,
family = binomial(link = "logit"))
## View
model_info(gm_icu)## Waiting for profiling to be done...## $`Model Summary`
##
## Call:
## glm(formula = NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS +
## SOFA + FIRST_CAREUNIT, family = binomial(link = "logit"),
## data = temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2894 -1.1498 0.6684 0.8166 1.5496
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.047484 0.108555 9.649 < 2e-16 ***
## GENDERM -0.349655 0.047894 -7.301 2.86e-13 ***
## AGE 0.281278 0.022939 12.262 < 2e-16 ***
## ETHNICITYOTHER 0.479055 0.112366 4.263 2.01e-05 ***
## ETHNICITYWHITE -0.119561 0.080072 -1.493 0.1354
## LANGUAGEOTHER -0.176088 0.070031 -2.514 0.0119 *
## MARITAL_STATUSSINGLE 0.521340 0.063216 8.247 < 2e-16 ***
## MARITAL_STATUSUNKNOWN -0.487552 0.112375 -4.339 1.43e-05 ***
## MARITAL_STATUSWIDOWED 0.126786 0.055320 2.292 0.0219 *
## SOFA 0.043068 0.007615 5.655 1.55e-08 ***
## FIRST_CAREUNITCSRU -1.223911 0.120647 -10.145 < 2e-16 ***
## FIRST_CAREUNITMICU 0.138801 0.065314 2.125 0.0336 *
## FIRST_CAREUNITSICU -0.828584 0.085265 -9.718 < 2e-16 ***
## FIRST_CAREUNITTSICU -0.646877 0.099174 -6.523 6.91e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13281 on 11158 degrees of freedom
## Residual deviance: 12508 on 11145 degrees of freedom
## AIC: 12536
##
## Number of Fisher Scoring iterations: 4
##
##
## $`OR Summary`
## OR 2.5 % 97.5 %
## (Intercept) 2.8504691 2.3067742 3.5305457
## GENDERM 0.7049316 0.6417118 0.7742484
## AGE 1.3248217 1.2666855 1.3858749
## ETHNICITYOTHER 1.6145473 1.2956016 2.0129012
## ETHNICITYWHITE 0.8873098 0.7572923 1.0366430
## LANGUAGEOTHER 0.8385444 0.7314905 0.9626274
## MARITAL_STATUSSINGLE 1.6842835 1.4888322 1.9075766
## MARITAL_STATUSUNKNOWN 0.6141281 0.4930084 0.7660711
## MARITAL_STATUSWIDOWED 1.1351745 1.0186282 1.2653224
## SOFA 1.0440092 1.0285974 1.0597692
## FIRST_CAREUNITCSRU 0.2940779 0.2319516 0.3722915
## FIRST_CAREUNITMICU 1.1488952 1.0100659 1.3048646
## FIRST_CAREUNITSICU 0.4366674 0.3693252 0.5159274
## FIRST_CAREUNITTSICU 0.5236788 0.4311703 0.6361001Logistic Regression ROC Curves (Team Level)
par(mfrow=c(2,3))
rocplot(as.numeric(predict(gm_a, type="response")), temp$NQF, col="blue", main = "Model 1a (Demographics)")
rocplot(as.numeric(predict(gm_b, type="response")), temp$NQF, col="blue", main = "Model 1b (Clinical Variables)")
rocplot(as.numeric(predict(gm_2a, type="response")), temp$NQF, col="blue", main = "Model 2a (Demographics + ICU)")
rocplot(as.numeric(predict(gm_2b, type="response")), temp$NQF, col="blue", main = "Model 2b (Clinical Variables + ICU)")
rocplot(as.numeric(predict(gm_icu, type="response")), temp$NQF, col="blue", main = "Model 3 (All Covariates)")
Results
tmp <- rbind(results, as.data.frame(cbind(mod, icc, std, lev, t_r2m, t_R2c, d_r2m, d_R2c)))
tmp$index <- c(0, 1, 1, 2, 2, 3, 0, 1, 1, 2, 2, 3)
## Correct variables
tmp$mod <- as.character(tmp$mod)
tmp$icc <- as.numeric(as.character(tmp$icc))
tmp$std <- as.numeric(as.character(tmp$std))
tmp$lev <- as.character(tmp$lev)
tmp$t_r2m <- as.numeric(as.character(tmp$t_r2m))
tmp$t_R2c <- as.numeric(as.character(tmp$t_R2c))
tmp$std <- round(tmp$std, 4)
tmp$icc <- round(tmp$icc, 4)
tmp$t_r2m <- round(tmp$t_r2m, 4)
tmp$t_R2c <- round(tmp$t_R2c, 4)
colnames(tmp) <- c("Model", "ICC", "Standard Deviation", "Level", colnames(tmp)[5:ncol(tmp)])
tmp$Letter <- c("Model: No Covariates", "Model: a (Demographics)",
"Model: b (Clinical Variables)", "Model: a (Demographics)",
"Model: b (Clinical Variables)", "Model: All Covariates", "Model: No Covariates", "Model: a (Demographics)", "Model: b (Clinical Variables)", "Model: a (Demographics)", "Model: b (Clinical Variables)", "Model: All Covariates")
colnames(tmp)[which(colnames(tmp) == "t_r2m")] <- "Marginal Coefficient of Determination"
colnames(tmp)[which(colnames(tmp) == "t_R2c")] <- "Conditional Coefficient of Determination"
print(tmp)## Model ICC Standard Deviation Level
## 1 Model 0 0.3938 1.4620 Provider
## 2 Model 1a 0.3911 1.4537 Provider
## 3 Model 1b 0.3938 1.4620 Provider
## 4 Model 2a 0.3987 1.4768 Provider
## 5 Model 2b 0.4054 1.4976 Provider
## 6 Model 3 0.3998 1.4802 Provider
## 7 Model 0 0.3794 1.4181 Team
## 8 Model 1a 0.3707 1.3920 Team
## 9 Model 1b 0.3787 1.4162 Team
## 10 Model 2a 0.3532 1.3403 Team
## 11 Model 2b 0.3677 1.3832 Team
## 12 Model 3 0.3547 1.3447 Team
## Marginal Coefficient of Determination
## 1 0.0000
## 2 0.0244
## 3 0.0000
## 4 0.0322
## 5 0.0076
## 6 0.0328
## 7 0.0000
## 8 0.0317
## 9 0.0013
## 10 0.0406
## 11 0.0074
## 12 0.0436
## Conditional Coefficient of Determination d_r2m
## 1 0.3938 0
## 2 0.4060 0.0216184415541336
## 3 0.3938 1.14290378300496e-07
## 4 0.4180 0.0285735799337614
## 5 0.4099 0.0067323447518925
## 6 0.4194 0.0290856883952417
## 7 0.3794 0
## 8 0.3906 0.025492554357851
## 9 0.3795 0.00104391524155094
## 10 0.3795 0.0325798018098079
## 11 0.3724 0.00595304054668518
## 12 0.3828 0.0349931785928889
## d_R2c index Letter
## 1 0.34752823595739 0 Model: No Covariates
## 2 0.359113489717591 1 Model: a (Demographics)
## 3 0.347527296071577 1 Model: b (Clinical Variables)
## 4 0.370611777313085 2 Model: a (Demographics)
## 5 0.362829255995969 2 Model: b (Clinical Variables)
## 6 0.37196177250928 3 Model: All Covariates
## 7 0.304086725763605 0 Model: No Covariates
## 8 0.314213109198477 1 Model: a (Demographics)
## 9 0.304241721258049 1 Model: b (Clinical Variables)
## 10 0.304173223878352 2 Model: a (Demographics)
## 11 0.297849819413253 2 Model: b (Clinical Variables)
## 12 0.307211798782246 3 Model: All CovariatesGraphs of Marginal and Conditional Coefficients of Determination
ggplot(tmp, aes(x=Model, y=`Marginal Coefficient of Determination`,
label = Level,
colour = Level)) +
geom_point(size=5, aes(col=Level), stat = "identity") +
geom_segment(aes(x=Model,
xend=Model,
y=min(`Marginal Coefficient of Determination`),
yend=max(`Marginal Coefficient of Determination`)+0.02),
linetype="dashed",
size=0.1) + # Draw dashed lines
labs(title="Dot Plot of Marginal Coefficient of Determination by Model", #,
subtitle="Care Team Level (Left), Provider Level (Right)") +#,
geom_text(size=4,aes(label=`Marginal Coefficient of Determination`, fontface = 2), position=position_dodge(width=0.9), hjust=0, vjust = -1,colour="black") +
coord_flip()
ggplot(tmp, aes(x=Model, y=`Conditional Coefficient of Determination`,
label = Level,
colour = Level)) +
geom_point(size=5, aes(col=Level), stat = "identity") +
geom_segment(aes(x=Model,
xend=Model,
y=min(`Conditional Coefficient of Determination`),
yend=max(`Conditional Coefficient of Determination`)+0.02),
linetype="dashed",
size=0.1) +
labs(title="Dot Plot of Conditional Coefficient of Determination by Model", #,
subtitle="Care Team Level (Left), Provider Level (Right)") +#,
geom_text(size=4,aes(label=`Conditional Coefficient of Determination`, fontface = 2),position=position_dodge(width=0.9), hjust=0, vjust = -1,colour="black") +
coord_flip()