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.

A mixed model is

\[\beta \sim \mathcal{N}(\mu, \sigma)\]

\[\mathbf{y} = \boldsymbol{X\beta} + \boldsymbol{Zu} + \boldsymbol{\varepsilon}\]

Statistics of Interest

  1. Intra Class Correlation

Calculated according to:

Wu S, Crespi CM, Wong WK. 2012. Comparison of methods for estimating the intraclass correlation coefficient for binary responses in cancer prevention cluster randomized trials. Contempory Clinical Trials 33: 869-880

  1. Standard Deviation of Random Effects

Calculated using:

Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48

  1. Marginal (variance explained by fixed effects) R^2

\[R_{GLMM(m)}² = \frac{(σ_f²)}{(σ_f² + σ_α² + σ_ε²)}\]

  1. Conditional (variance explained by the entire model) R^2

\[R_{GLMM(c)}² = \frac{(σ_f² + σ_α²)}{(σ_f² + σ_α² + σ_ε²)}\]

Calculated according to:

Nakagawa, S., Johnson, P.C.D., Schielzeth, H. (2017) The coefficient of determination R² and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. J. R. Soc. Interface 14: 20170213.

Model 0 (No Covariates)

  1. Care Provider (Random)

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Language
  5. Marrital Status
  6. Care Provider (Random)

Model 1b (Clinical Variables)

  1. Sequential Organ Failure Assessment (SOFA)
  2. Care Provider (Random)

Model 2a (Demographics + ICU)

  1. Demographics (Model 1a)
  2. First Intensive Care Unit
  3. Care Provider (Random)

Model 2b (Clinical Variables + ICU)

  1. SOFA (Model 1b)
  2. First Intensive Care Unit
  3. 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.

  1. Demographics (as in Model 1)
  2. First Intensive Care Unit
  3. SOFA Score (As Model 2)
  4. Care Provider (Random)

Attending check

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 <- function(dat){
    ## Create Caremeasure Variable
    dat$NQF <- rep(0, each = nrow(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$CIM.machine == 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 <- function(dat){
    ## Create Caremeasure Variable
    dat$NQF <- rep(0, each = nrow(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$CIM.machine == 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 <- 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 <- 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 <- 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 & Cleansing

## Latest Dataset of NeuroNER Predictions
#dat <- read.csv("~/nqf_dat.csv", header = T, stringsAsFactors = F)
dat <- read.csv("~/nqf_caregivers/data/20180607_EOL_data_ICU.csv", header = T, stringsAsFactors = F)
dim(dat)
## [1] 10250    57
## Note: Notes had been logged by multiple Care providers, we will reintroduce those annotations

## Load Labeled Note Data for NQF Caremeasure Cohort (From NOTEEVENTS table)
tmp <- read.csv("~/nqf_caregivers/data/note_labels_over75.csv", header = T, stringsAsFactors = F)
dim(tmp)
## [1] 11575    25
## Keep only TEXT and ROW_ID from tmp
tmp <- tmp[ ,c("ROW_ID", "TEXT")]

## Inner join
dat <- merge(tmp, dat, by = "TEXT")

## Clean tmp
rm(tmp)

## Check column names
colnames(dat)
##  [1] "TEXT"                 "ROW_ID.x"             "SUBJECT_ID"          
##  [4] "HADM_ID"              "ROW_ID.y"             "CHARTDATE"           
##  [7] "CHARTTIME"            "STORETIME"            "CATEGORY"            
## [10] "DESCRIPTION"          "CGID"                 "ISERROR"             
## [13] "ADMITTIME"            "DISCHTIME"            "DEATHTIME"           
## [16] "ADMISSION_TYPE"       "ADMISSION_LOCATION"   "DISCHARGE_LOCATION"  
## [19] "INSURANCE"            "LANGUAGE"             "RELIGION"            
## [22] "MARITAL_STATUS"       "ETHNICITY"            "EDREGTIME"           
## [25] "EDOUTTIME"            "DIAGNOSIS"            "HOSPITAL_EXPIRE_FLAG"
## [28] "HAS_CHARTEVENTS_DATA" "GENDER"               "DOB"                 
## [31] "DOD"                  "DOD_HOSP"             "DOD_SSN"             
## [34] "EXPIRE_FLAG"          "ICUSTAY_ID"           "DBSOURCE"            
## [37] "FIRST_CAREUNIT"       "LAST_CAREUNIT"        "FIRST_WARDID"        
## [40] "LAST_WARDID"          "INTIME"               "OUTTIME"             
## [43] "LOS"                  "AGE"                  "ADMISSION_NUMBER"    
## [46] "DAYS_UNTIL_DEATH"     "TIME_SINCE_ADMIT"     "CGID.1"              
## [49] "HADM_ID.1"            "FAM.machine"          "CIM.machine"         
## [52] "LIM.machine"          "CAR.machine"          "COD.machine"         
## [55] "check.CGID"           "check.dadm_id"        "CIM.or.FAM"          
## [58] "Died.in.Hospital"
## What is HADM_ID.1?
head(table(dat$HADM_ID.1))
## 
##   #N/A 100102 100153 100347 100391 100525 
##     32      8      3     12     15      2
## What is HADM_ID?
head(table(dat$HADM_ID))
## 
## 100102 100153 100347 100391 100525 100575 
##      8      3     12     15      2     11
## #N/A? Clean HADM_ID.1
dat$HADM_ID.1 <- NULL

## What is CGID.1?
head(table(dat$CGID.1))
## 
##  #N/A 14010 14022 14037 14045 14056 
##    32    22     1    93    37    14
## What is CGID
head(table(dat$CGID))
## 
## 14010 14022 14037 14045 14056 14080 
##    22     1    93    37    14     6
## #N/A? Clean CGID.1
dat$CGID.1 <- NULL

## What is check.CGID
head(table(dat$check.CGID))
##     0 
## 11575
## Clean it
dat$check.CGID <- NULL

## What is check.dadm_id?
head(table(dat$check.dadm_id))
##     0 
## 11575
## Clean it
dat$check.dadm_id <- NULL

## Clean column names
dat$ROW_ID.y <- NULL

colnames(dat)[which(colnames(dat) == "ROW_ID.x")] <- "ROW_ID"


## Load CAREGIVERS Table for join on CGID
cg <- read.csv("~/nqf_caregivers/data/mimic/CAREGIVERS.csv", 
               header = T, stringsAsFactors = F)

## Change column name of "NOTEEVENTS.DESCRIPTION" to explicitly mention that it describes the note
colnames(dat)[which(colnames(dat) == "DESCRIPTION")] <- "NOTE_DESCRIPTION"

## Change column name of "CAREGIVERS. DESCRIPTION" to explicitly mention that it describes the careprovider
colnames(cg)[which(colnames(cg) == "DESCRIPTION")] <- "CG_DESCRIPTION"

## Remove ROW_ID from CG
cg$ROW_ID <- NULL

## Remove TEXT
dat$TEXT <- NULL

## Merge to caregivers
dat <- merge(dat, cg, by = "CGID")
dim(dat)
## [1] 11575    54
## Clean CG
rm(cg)

sofa <- read.csv("~/nqf_caregivers/data/sofa.csv", header = T, stringsAsFactors = F)
#oasis <- read.csv("~/nqf_caregivers/data/oasis.csv", header = T, stringsAsFactors = F)
#saps <- read.csv("~/nqf_caregivers/data/saps.csv", header = T, stringsAsFactors = F)

colnames(sofa) <- toupper(colnames(sofa))

dat <- merge(dat, sofa, by = c("SUBJECT_ID", "HADM_ID", "ICUSTAY_ID"))
dim(dat)
## [1] 11575    61
## Clean environment
rm(sofa)

## Clean ethnicity to Black/White/Other
dat[(grepl("WHITE|PORTUGUESE", dat$ETHNICITY)),]$ETHNICITY <- "WHITE" 
dat[(grepl("ASIAN", dat$ETHNICITY)),]$ETHNICITY <- "OTHER" 
dat[(grepl("BLACK", dat$ETHNICITY)),]$ETHNICITY <- "BLACK" 
dat[(grepl("HISPANIC", dat$ETHNICITY)),]$ETHNICITY <- "OTHER"
dat[(grepl("MIDDLE|NATIVE|MULTI|DECLINED|UNABLE|OTHER|NOT", dat$ETHNICITY)),]$ETHNICITY <- "OTHER"

## Clean Marital Status to Married, Single, Widowed, Unknown
dat$MARITAL_STATUS[dat$MARITAL_STATUS == ""] <- "UNKNOWN (DEFAULT)"
dat$MARITAL_STATUS[dat$MARITAL_STATUS == "UNKNOWN (DEFAULT)"] <- "UNKNOWN"
dat$MARITAL_STATUS[dat$MARITAL_STATUS == "SEPARATED"] <- "SINGLE"
dat$MARITAL_STATUS[dat$MARITAL_STATUS == "DIVORCED"] <- "SINGLE"

## Clean admission data
dat$ADMISSION_TYPE[dat$ADMISSION_TYPE == "URGENT"] <- "EMERGENCY"

## Clean Languages
dat$LANGUAGE <- ifelse(dat$LANGUAGE == "ENGL", "ENGL", "OTHER")

#write.csv(dat, file = "~/nqf_caregivers/data/nqf_dat.csv", row.names = F)

NQF Care Measure Check

prov <- provider_caremeasure_check(attending_check(dat))

prov <- expire_check(prov)

team <- team_caremeasure_check(attending_check(dat))

team <- expire_check(team)

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.

temp <- prov

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

Data Standardization and Factoring (for modeling)

Continuous data, Age and SOFA, will be standardized in the form:

\[z = \frac{x - \mu}{\sigma}\]

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)

##  as.numeric(attr(summary(gm_initial)$varcor$CGID, "stddev"))

## 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 ' ' 1
sjstats::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:
plot_model(gm_initial, type = "diag")
## $CGID

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Language
  5. Marrital Status
  6. 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.205
tab_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^2

Model 1b (Clinical Variables)

  1. SOFA
  2. 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.311
tab_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^2

Model 2a (Demographics + ICU)

  1. Demographics (Model 1a)
  2. First Intensive Care Unit
  3. 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 it
tab_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^2

Model 2b (Clinical Variables + ICU)

  1. Clinical Variables (Model 1b)
  2. First Intensive Care Unit
  3. 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.557
tab_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^2

Model 3 (Demographics + Clinical Variables + ICU)

  1. Demographics (as in Model 1a)
  2. SOFA Score (As Model 1b)
  3. First Intensive Care Unit
  4. 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 it
tab_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.37196177250928

ROC Curves

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

Baseline Model (Provider Random Effects)

Not meaningful.

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Language
  5. Marrital Status
  6. 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.2247760

Model 1b (Clinical Variables)

  1. SOFA
  2. 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.011798

Model 2a (Demographics + ICU)

  1. Demographics (Model 1a)
  2. First Intensive Care Unit
  3. 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.9037148

Model 2b (Clinical Variables + ICU)

  1. Clinical Variables (Model 1b)
  2. First Intensive Care Unit
  3. 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.8366620

Model 3 (Demographics + Clinical Variables + ICU)

  1. Demographics (as in Model 1a)
  2. SOFA Score (As Model 1b)
  3. First Intensive Care Unit
  4. 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.9104643
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.

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

Data Standardization and Factoring (for modeling)

Continuous data, Age and SOFA, will be standardized in the form:

\[z = \frac{x - \mu}{\sigma}\]

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)

##  as.numeric(attr(summary(gm_initial)$varcor$CGID, "stddev"))

## 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 ' ' 1
sjstats::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)

  1. Gender
  2. Age
  3. Ethnicity
  4. Language
  5. Marrital Status
  6. 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.202
tab_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^2

Model 1b (Clinical Variables)

  1. SOFA
  2. 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.334
tab_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^2

Model 2a (Demographics + ICU)

  1. Demographics (Model 1a)
  2. First Intensive Care Unit
  3. 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 it
tab_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^2

Model 2b (Clinical Variables + ICU)

  1. Clinical Variables (Model 1b)
  2. First Intensive Care Unit
  3. 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.596
tab_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^2

Model 3 (Demographics + Clinical Variables + ICU)

  1. Demographics (as in Model 1a)
  2. SOFA Score (As Model 1b)
  3. First Intensive Care Unit
  4. 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 it
tab_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^2

ROC Curves

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

Baseline Model (Provider Random Effects)

Not meaningful.

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Languagee
  5. Marrital Status
  6. 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.1782824

Model 1b (Clinical Variables)

  1. SOFA
  2. 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.047725

Model 2a (Demographics + ICU)

  1. Demographics (Model 1a)
  2. First Intensive Care Unit
  3. 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.6224501

Model 2b (Clinical Variables + ICU)

  1. Clinical Variables (Model 1b)
  2. First Intensive Care Unit
  3. 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.5772770

Model 3 (Demographics + Clinical Variables + ICU)

  1. Demographics (as in Model 1a)
  2. SOFA Score (As Model 1b)
  3. First Intensive Care Unit
  4. 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.6361001
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 Covariates
##theme_set(theme_classic())

ggplot(tmp, aes(x=Model, y=`Marginal Coefficient of Determination`, 
  label = Level, 
  colour = Level)) + 
  geom_point(size=5, aes(col=Level), stat = "identity") +   # Draw points
  #geom_text(color="black", size=3) + 
  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)") +#, 
       #caption="source: mpg")# +  
  # Change fontface. Allowed values : 1(normal),
  # 2(bold), 3(italic), 4(bold.italic)
    geom_text(size=4,aes(label=`Marginal Coefficient of Determination`, fontface = 2), position=position_dodge(width=0.9), hjust=0, vjust = -1,colour="black") +
    #geom_text(size=2,aes(label=Level, fontface = 2),position=position_dodge(width=0.9), hjust=-0.45, 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") +   # Draw points
  #geom_text(color="black", size=3) + 
  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) +   # Draw dashed lines
  labs(title="Dot Plot of Conditional Coefficient of Determination by Model",  #, 
       subtitle="Care Team Level (Left), Provider Level (Right)") +#, 
       #caption="source: mpg")# +  
  # Change fontface. Allowed values : 1(normal),
  # 2(bold), 3(italic), 4(bold.italic)
    geom_text(size=4,aes(label=`Conditional Coefficient of Determination`, fontface = 2),position=position_dodge(width=0.9), hjust=0, vjust = -1,colour="black") +
    #geom_text(size=2,aes(label=Level, fontface = 2),position=position_dodge(width=0.9), hjust=-0.45, vjust = 1,colour="black") +
  coord_flip()

References

  1. Vincent JL, de Mendonça A, Cantraine F, Moreno R, Takala J, Suter PM, Sprung CL, Colardyn F, Blecher S. Use of the SOFA score to assess the incidence of organ dysfunction/failure in intensive care units: results of a multicenter, prospective study. Working group on “sepsis-related problems” of the European Society of Intensive Care Medicine. Crit Care Med 1998 Nov;26(11):1793-800. PMID 9824069.

  2. Moreno R, Vincent JL, Matos R, Mendonça A, Cantraine F, Thijs L, Takala J, Sprung C, Antonelli M, Bruining H, Willatts S. The use of maximum SOFA score to quantify organ dysfunction/failure in intensive care. Results of a prospective, multicentre study. Working Group on Sepsis related Problems of the ESICM. Intensive Care Med 1999 Jul;25(7):686-96. PMID 10470572.

  3. Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B 57, 289–300.