Care Measure Variation

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.

It is important to note that 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.

Model 0 (No Covariates)

Care Provider (Random)

Model 1a (Patient Demographics)

Gender
Age
Ethnicity
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)

Libraries

library("ggplot2") ## Graphing
library("GGally") ## Extension to ggplot
library("reshape2") ## for melt()
library("lme4") ## for Hierachal Models
library("sjPlot") ## Lovely Presentation of Model Output
library("rcompanion") ## Pairwise nominal test
library("ROCR")
library("pROC")

Utility Functions

caremeasure_rate_check

caremeasure_rate_check will go through each hospital admission and keep only those admissions with attendings who have logged a patient note that was captured by MIMIC-III. In addition to checking for an attending physician, it will also check to see if the caremeasure was implimented.

caremeasure_rate_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)){
            ## Aggregate
            tmp_frame <- aggregate(cbind(CIM.machine,
                       AGE,
                       SOFA) ~ 
                     ETHNICITY +
                     GENDER +
                     MARITAL_STATUS +
                     FIRST_CAREUNIT +
                     HADM_ID +
                     SUBJECT_ID +
                     CGID,
                 data = tmp_frame, 
                 FUN = mean)
            
            res <- rbind(res, tmp_frame)
        }
    }
    ## Return control to outer level
    return(res)
}

caremeasure_check

caremeasure_check()

caremeasure_check <- function(dat){
    ## Create Caremeasure Variable
    dat$NQF <- 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 care providers are "Attending"
        if (any("Attending" %in% tmp_frame$CG_DESCRIPTION)){
            ## Check if caremeasure was implemented
            if (any(tmp_frame$CIM.machine == 1)){
                tmp_frame$NQF <- rep(1, each = nrow(tmp_frame))
            }
            ## add hospital admission to results
            res <- rbind(res, tmp_frame)
        }
    }
    ## Return control to outer level
    return(res)
}

plotDat

A simple plotting function.

plotDat <- function(dat, column, x_col, bs, mn, xl, yl){
  tmp <- as.matrix(table(dat[[column]], dat[[x_col]]))
  prop <- prop.table(tmp, margin = 2)#2 for column-wise proportions
  par(mar = c(5.0, 4.0, 4.0, 15), xpd = TRUE)
  barplot(prop, col = cm.colors(length(rownames(prop))), beside = bs, width = 2, main = mn, xlab = xl, ylab = yl)
  legend("topright", inset = c(-0.90,0), fill = cm.colors(length(rownames(prop))), legend=rownames(prop))
}

detach_packages

detach_packages will keep only base R packages, and will remove all other supplementary packages to avoid functional conflicts.

zero_range

# https://stackoverflow.com/questions/4752275/test-for-equality-among-all-elements-of-a-single-vector
# Determine if range of vector is FP 0.
zero_range <- function(x, tol = .Machine$double.eps ^ 0.5) {
    if (length(x) == 1) return(TRUE)
    x <- range(x) / mean(x)
    isTRUE(all.equal(x[1], x[2], tolerance = tol))
}

t_tests

t_tests large scale paired t-testing with p-value adjustment.

t_tests <- function(dat, tmp, var_name){
    #Create empty list 
    results <- list()
    
    for (i in 1:length(tmp)){
        # Subset data
        df <- dat[dat[[var_name]] %in% as.character(tmp[[i]]),]
        
        # Subset data for both cytokines
        tmp_one <- df[df[[var_name]] == levels(factor(df[[var_name]]))[1], ]$NQF
        tmp_two <- df[df[[var_name]] == levels(factor(df[[var_name]]))[1], ]$NQF
        
        # Check if the range is zero-- t.test will throw error if so
        if (zero_range(tmp_one) | zero_range(tmp_two)){
            # If zero, NA
            results[[i]] <- NA
        } else {
            results[[i]] <- t.test(df$NQF ~ df[[var_name]])
        }
    }
    
    # Rename list ,entries based on gene pairs  
    names(results) <- paste(matrix(unlist(tmp), ncol = 2, byrow = TRUE)[ ,1],
                            matrix(unlist(tmp), ncol = 2, byrow = TRUE)[ ,2], sep = " vs. ")
    
    # Convert to data frame
    results <- as.data.frame(unlist(results))
    
    # Add names as variable
    results$test_pair <- row.names(results)
    
    # Keep only p values
    results <- results[grepl("p.value", results$test_pair), ]
    row.names(results) <- 1:nrow(results)
    
    # Column names or orientation
    colnames(results) <- c("p.value", "test_pair")
    results <- results[ ,c("test_pair", "p.value")]
    
    # Convert p-value to numeric (from factor)
    results$p.value <- as.numeric(as.character(results$p.value))
    
    # Adjust p value using Benjamini & Hochberg (1995) "fdr" method
    results$adjusted.p <- p.adjust(results$p.value, method = "fdr", n = nrow(results))
    
    # Significance codes
    results$signif <- ifelse(results$adjusted.p <= 0.1, '.', ' ')
    results$signif <- ifelse(results$adjusted.p <= 0.05, '*', results$signif)
    results$signif <- ifelse(results$adjusted.p <= 0.01, "**", results$signif)
    results$signif <- ifelse(results$adjusted.p <= 0.001, "***", results$signif)
    
    return(results)
}

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.8, 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)
}

Load Data and Merge to Care Providers Table

The latest dataset with NeuroNER predictions will be used. Multiple clinicians are associated in the database with identical notes, and those notes will be re-introduced to the data set by performing an inner join between old and new data, which will apply labels to the old data (which contains duplicate notes associated with different care providers).

## Latest Dataset of NeuroNER Predictions
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)

Load Severity of Illness Data

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)

Data Cleaning

## 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"

Check Attending and Caremeasure Implementation

temp <- caremeasure_rate_check(dat)

## Change caremeasure name
colnames(temp)[which(colnames(temp) == "CIM.machine")] <- "NQF"

cat(length(unique(dat$CGID)) - length(unique(temp$CGID)), "Clinicians dropped due to no attending data in notes.\n")

## 4 Clinicians dropped due to no attending data in notes.

cat(length(unique(dat$SUBJECT_ID)) - length(unique(temp$SUBJECT_ID)), "Patients dropped due to no attending data in notes.\n")

## 78 Patients dropped due to no attending data in notes.

cat(length(unique(dat$HADM_ID)) - length(unique(temp$HADM_ID)), "Hospital Admissions dropped due to no attending data in notes.\n")

## 95 Hospital Admissions dropped due to no attending data in notes.

Data Exploration

boxplot(temp$NQF ~ temp$FIRST_CAREUNIT)

test <- t_tests(temp, combn(unique(temp$FIRST_CAREUNIT), 2, simplify = F), "FIRST_CAREUNIT")

test

##                 test_pair      p.value   adjusted.p signif
## 1    MICU vs. CCU.p.value 4.363966e-02 4.848851e-02      *
## 2   MICU vs. CSRU.p.value 1.096719e-16 1.096719e-15    ***
## 3   MICU vs. SICU.p.value 1.306458e-09 4.354861e-09    ***
## 4  MICU vs. TSICU.p.value 6.116334e-07 1.529084e-06    ***
## 5    CCU vs. CSRU.p.value 3.433527e-11 1.716763e-10    ***
## 6    CCU vs. SICU.p.value 7.899781e-04 1.128540e-03     **
## 7   CCU vs. TSICU.p.value 3.112150e-03 3.890188e-03     **
## 8   CSRU vs. SICU.p.value 2.734587e-05 5.469174e-05    ***
## 9  CSRU vs. TSICU.p.value 1.132300e-04 1.887167e-04    ***
## 10 SICU vs. TSICU.p.value 9.736742e-01 9.736742e-01

boxplot(temp$NQF ~ temp$ETHNICITY)

test <- t_tests(temp, combn(unique(temp$ETHNICITY), 2, simplify = F), "ETHNICITY")

test

##                 test_pair    p.value adjusted.p signif
## 1 WHITE vs. OTHER.p.value 0.68970852 0.68970852       
## 2 WHITE vs. BLACK.p.value 0.02421985 0.07265955      .
## 3 OTHER vs. BLACK.p.value 0.14590895 0.21886343

boxplot(temp$NQF ~ temp$GENDER)

t.test(temp$NQF ~ temp$GENDER)

## 
##  Welch Two Sample t-test
## 
## data:  temp$NQF by temp$GENDER
## t = 6.3717, df = 4947.9, p-value = 2.039e-10
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.05959236 0.11256032
## sample estimates:
## mean in group F mean in group M 
##       0.5154297       0.4293534

test <- t_tests(temp, combn(unique(temp$GENDER), 2, simplify = F), "GENDER")

test

##         test_pair      p.value   adjusted.p signif
## 1 F vs. M.p.value 2.039191e-10 2.039191e-10    ***

boxplot(temp$NQF ~ temp$MARITAL_STATUS)

test <- t_tests(temp, combn(unique(temp$MARITAL_STATUS), 2, simplify = F), "MARITAL_STATUS")

test

##                     test_pair      p.value   adjusted.p signif
## 1 MARRIED vs. WIDOWED.p.value 6.407815e-09 3.844689e-08    ***
## 2  MARRIED vs. SINGLE.p.value 4.492615e-07 1.347785e-06    ***
## 3 MARRIED vs. UNKNOWN.p.value 9.614203e-01 9.614203e-01       
## 4  WIDOWED vs. SINGLE.p.value 9.368273e-01 9.614203e-01       
## 5 WIDOWED vs. UNKNOWN.p.value 1.257390e-02 2.112715e-02      *
## 6  SINGLE vs. UNKNOWN.p.value 1.408477e-02 2.112715e-02      *

print("Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1")

## [1] "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"

Data Standardization and Factoring (for modeling)

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

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

In this way, we can analyze the data in units of their deviation from the mean.

temp$AGE <- (temp$AGE - mean(temp$AGE))/sd(temp$AGE)
temp$SOFA <- (temp$SOFA - mean(temp$SOFA)/sd(temp$SOFA))

## Factor
temp <- within(temp, {
    ETHNICITY <- factor(ETHNICITY)
    GENDER <- factor(GENDER)
    MARITAL_STATUS <- factor(MARITAL_STATUS)
    FIRST_CAREUNIT <- factor(FIRST_CAREUNIT)
    SUBJECT_ID <- factor(SUBJECT_ID)
    CGID <- factor(CGID)
})

Baseline Model (Provider Random Effects)

m_initial <- lmer(NQF ~ (1 | CGID),
                   data = temp)

## View
summary(m_initial)

## Linear mixed model fit by REML ['lmerMod']
## Formula: NQF ~ (1 | CGID)
##    Data: temp
## 
## REML criterion at convergence: 6540.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7054 -0.9195 -0.2895  0.9609  1.7612 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CGID     (Intercept) 0.0254   0.1594  
##  Residual             0.2053   0.4531  
## Number of obs: 4950, groups:  CGID, 493
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.47842    0.01092    43.8

## Intraclass Correlation
sjstats::icc(m_initial)

## 
## Linear mixed model
##  Family: gaussian (identity)
## Formula: NQF ~ (1 | CGID)
## 
##   ICC (CGID): 0.110119

Model 1a (Patient Demographics)

Gender
Age
Ethnicity
Marrital Status
Care Provider (Random)

m_a <- lmer(NQF ~ GENDER +
                     AGE +
                     ETHNICITY +
                     MARITAL_STATUS +
                     (1 | CGID), 
                   data = temp)

## View
summary(m_a)

## Linear mixed model fit by REML ['lmerMod']
## Formula: NQF ~ GENDER + AGE + ETHNICITY + MARITAL_STATUS + (1 | CGID)
##    Data: temp
## 
## REML criterion at convergence: 6420.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7949 -0.9044 -0.2545  0.9598  1.9602 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CGID     (Intercept) 0.02349  0.1533  
##  Residual             0.19926  0.4464  
## Number of obs: 4950, groups:  CGID, 493
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            0.509208   0.027080  18.804
## GENDERM               -0.063470   0.014112  -4.498
## AGE                    0.069035   0.006649  10.383
## ETHNICITYOTHER        -0.004560   0.030375  -0.150
## ETHNICITYWHITE        -0.023341   0.023133  -1.009
## MARITAL_STATUSSINGLE   0.060012   0.018332   3.274
## MARITAL_STATUSUNKNOWN -0.017329   0.033903  -0.511
## MARITAL_STATUSWIDOWED  0.028246   0.016573   1.704
## 
## Correlation of Fixed Effects:
##                 (Intr) GENDER AGE    ETHNICITYO ETHNICITYW MARITAL_STATUSS
## GENDERM         -0.322                                                    
## AGE              0.115 -0.013                                             
## ETHNICITYOT     -0.607 -0.052 -0.054                                      
## ETHNICITYWH     -0.783 -0.068 -0.089  0.689                               
## MARITAL_STATUSS -0.362  0.257  0.007  0.085      0.073                    
## MARITAL_STATUSU -0.151  0.089 -0.074 -0.100      0.026      0.174         
## MARITAL_STATUSW -0.455  0.369 -0.177  0.104      0.096      0.437         
##                 MARITAL_STATUSU
## GENDERM                        
## AGE                            
## ETHNICITYOT                    
## ETHNICITYWH                    
## MARITAL_STATUSS                
## MARITAL_STATUSU                
## MARITAL_STATUSW  0.217

sjstats::icc(m_a)

## 
## Linear mixed model
##  Family: gaussian (identity)
## Formula: NQF ~ GENDER + AGE + ETHNICITY + MARITAL_STATUS + (1 | CGID)
## 
##   ICC (CGID): 0.105475

Model 1b (Clinical Variables)

SOFA
Care Provider (Random)

m_b <- lmer(NQF ~ SOFA +
                     (1 | CGID), 
                   data = temp)

## View
summary(m_b)

## Linear mixed model fit by REML ['lmerMod']
## Formula: NQF ~ SOFA + (1 | CGID)
##    Data: temp
## 
## REML criterion at convergence: 6550.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7096 -0.9168 -0.2940  0.9600  1.7594 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CGID     (Intercept) 0.02544  0.1595  
##  Residual             0.20528  0.4531  
## Number of obs: 4950, groups:  CGID, 493
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.481713   0.013064   36.87
## SOFA        -0.001031   0.002242   -0.46
## 
## Correlation of Fixed Effects:
##      (Intr)
## SOFA -0.548

sjstats::icc(m_b)

## 
## Linear mixed model
##  Family: gaussian (identity)
## Formula: NQF ~ SOFA + (1 | CGID)
## 
##   ICC (CGID): 0.110254

Model 2a (Demographics + ICU)

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

m_2a <- lmer(NQF ~ GENDER +
                      AGE +
                      ETHNICITY +
                      MARITAL_STATUS +
                      FIRST_CAREUNIT +
                      (1 | CGID), 
                   data = temp)

## View
summary(m_2a)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## NQF ~ GENDER + AGE + ETHNICITY + MARITAL_STATUS + FIRST_CAREUNIT +  
##     (1 | CGID)
##    Data: temp
## 
## REML criterion at convergence: 6405.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8291 -0.9152 -0.2343  0.9633  1.9654 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CGID     (Intercept) 0.0176   0.1326  
##  Residual             0.2002   0.4474  
## Number of obs: 4950, groups:  CGID, 493
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            0.482593   0.032097  15.036
## GENDERM               -0.060634   0.014144  -4.287
## AGE                    0.070705   0.006659  10.617
## ETHNICITYOTHER         0.005973   0.030468   0.196
## ETHNICITYWHITE        -0.021393   0.023140  -0.924
## MARITAL_STATUSSINGLE   0.059841   0.018331   3.264
## MARITAL_STATUSUNKNOWN -0.011231   0.034005  -0.330
## MARITAL_STATUSWIDOWED  0.031766   0.016596   1.914
## FIRST_CAREUNITCSRU    -0.146203   0.039847  -3.669
## FIRST_CAREUNITMICU     0.057552   0.020878   2.757
## FIRST_CAREUNITSICU    -0.029194   0.027747  -1.052
## FIRST_CAREUNITTSICU   -0.011955   0.032465  -0.368
## 
## Correlation of Fixed Effects:
##                 (Intr) GENDER AGE    ETHNICITYO ETHNICITYW MARITAL_STATUSS
## GENDERM         -0.289                                                    
## AGE              0.066 -0.013                                             
## ETHNICITYOT     -0.540 -0.050 -0.050                                      
## ETHNICITYWH     -0.677 -0.066 -0.087  0.690                               
## MARITAL_STATUSS -0.314  0.256  0.006  0.086      0.074                    
## MARITAL_STATUSU -0.101  0.086 -0.074 -0.100      0.024      0.171         
## MARITAL_STATUSW -0.373  0.369 -0.179  0.105      0.098      0.436         
## FIRST_CAREUNITC -0.245 -0.035  0.029  0.003     -0.003      0.016         
## FIRST_CAREUNITM -0.531  0.026  0.061  0.065      0.033      0.008         
## FIRST_CAREUNITS -0.402  0.044  0.031 -0.002      0.001      0.014         
## FIRST_CAREUNITT -0.377  0.035  0.012  0.031      0.028      0.031         
##                 MARITAL_STATUSU MARITAL_STATUSW FIRST_CAREUNITC
## GENDERM                                                        
## AGE                                                            
## ETHNICITYOT                                                    
## ETHNICITYWH                                                    
## MARITAL_STATUSS                                                
## MARITAL_STATUSU                                                
## MARITAL_STATUSW  0.216                                         
## FIRST_CAREUNITC -0.040          -0.051                         
## FIRST_CAREUNITM -0.024          -0.017           0.387         
## FIRST_CAREUNITS -0.054          -0.029           0.345         
## FIRST_CAREUNITT -0.084           0.015           0.304         
##                 FIRST_CAREUNITM FIRST_CAREUNITS
## GENDERM                                        
## AGE                                            
## ETHNICITYOT                                    
## ETHNICITYWH                                    
## MARITAL_STATUSS                                
## MARITAL_STATUSU                                
## MARITAL_STATUSW                                
## FIRST_CAREUNITC                                
## FIRST_CAREUNITM                                
## FIRST_CAREUNITS  0.596                         
## FIRST_CAREUNITT  0.505           0.494

sjstats::icc(m_2a)

## 
## Linear mixed model
##  Family: gaussian (identity)
## Formula: NQF ~ GENDER + AGE + ETHNICITY + MARITAL_STATUS + FIRST_CAREUNIT + (1 | CGID)
## 
##   ICC (CGID): 0.080788

Model 2b (Clinical Variables + ICU)

Clinical Variables (Model 1b)
First Intensive Care Unit
Care Provider (Random)

m_2b <- lmer(NQF ~ SOFA +
                  FIRST_CAREUNIT +
                      (1 | CGID), 
                   data = temp)

## View
summary(m_2b)

## Linear mixed model fit by REML ['lmerMod']
## Formula: NQF ~ SOFA + FIRST_CAREUNIT + (1 | CGID)
##    Data: temp
## 
## REML criterion at convergence: 6539.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6999 -0.9118 -0.2877  0.9650  1.8780 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CGID     (Intercept) 0.01944  0.1394  
##  Residual             0.20634  0.4542  
## Number of obs: 4950, groups:  CGID, 493
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)          0.466809   0.021207  22.012
## SOFA                -0.001245   0.002254  -0.552
## FIRST_CAREUNITCSRU  -0.156965   0.040487  -3.877
## FIRST_CAREUNITMICU   0.048269   0.021171   2.280
## FIRST_CAREUNITSICU  -0.026038   0.028227  -0.922
## FIRST_CAREUNITTSICU -0.012234   0.032995  -0.371
## 
## Correlation of Fixed Effects:
##                 (Intr) SOFA   FIRST_CAREUNITC FIRST_CAREUNITM
## SOFA            -0.342                                       
## FIRST_CAREUNITC -0.395 -0.029                                
## FIRST_CAREUNITM -0.765 -0.020  0.387                         
## FIRST_CAREUNITS -0.628  0.054  0.342           0.594         
## FIRST_CAREUNITT -0.551  0.055  0.303           0.503         
##                 FIRST_CAREUNITS
## SOFA                           
## FIRST_CAREUNITC                
## FIRST_CAREUNITM                
## FIRST_CAREUNITS                
## FIRST_CAREUNITT  0.496

sjstats::icc(m_2b)

## 
## Linear mixed model
##  Family: gaussian (identity)
## Formula: NQF ~ SOFA + FIRST_CAREUNIT + (1 | CGID)
## 
##   ICC (CGID): 0.086102

Model 3 (Demographics + Clinical Variables + ICU)

Demographics (as in Model 1a)
SOFA Score (As Model 1b)
First Intensive Care Unit
Care Provider (Random)

m_icu <- lmer(NQF ~ 
                   GENDER +
                   AGE +
                   ETHNICITY +
                   MARITAL_STATUS +
                   SOFA +
                   FIRST_CAREUNIT +
                   (1 | CGID), 
                data = temp)

## View
summary(m_icu)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## NQF ~ GENDER + AGE + ETHNICITY + MARITAL_STATUS + SOFA + FIRST_CAREUNIT +  
##     (1 | CGID)
##    Data: temp
## 
## REML criterion at convergence: 6414.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8376 -0.9138 -0.2331  0.9600  1.9782 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  CGID     (Intercept) 0.01761  0.1327  
##  Residual             0.20020  0.4474  
## Number of obs: 4950, groups:  CGID, 493
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            0.475932   0.032876  14.477
## GENDERM               -0.062025   0.014222  -4.361
## AGE                    0.071088   0.006672  10.655
## ETHNICITYOTHER         0.005144   0.030481   0.169
## ETHNICITYWHITE        -0.020753   0.023150  -0.896
## MARITAL_STATUSSINGLE   0.060192   0.018335   3.283
## MARITAL_STATUSUNKNOWN -0.009927   0.034033  -0.292
## MARITAL_STATUSWIDOWED  0.032067   0.016599   1.932
## SOFA                   0.002102   0.002249   0.935
## FIRST_CAREUNITCSRU    -0.147188   0.039864  -3.692
## FIRST_CAREUNITMICU     0.057072   0.020885   2.733
## FIRST_CAREUNITSICU    -0.027924   0.027781  -1.005
## FIRST_CAREUNITTSICU   -0.010396   0.032508  -0.320

## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it

sjstats::icc(m_icu)

## 
## Linear mixed model
##  Family: gaussian (identity)
## Formula: NQF ~ GENDER + AGE + ETHNICITY + MARITAL_STATUS + SOFA + FIRST_CAREUNIT + (1 | CGID)
## 
##   ICC (CGID): 0.080869

Generalized Linear Models

Note: aggregation occurs as the hospital admission level. If any of the careproviders documented preferences, it was considered a positive implementation of the care measure.

## Extra libs
library("DescTools")

## Warning: package 'DescTools' was built under R version 3.4.4

library("pscl")
library("PresenceAbsence")

temp <- caremeasure_check(dat)

cat(length(unique(dat$CGID)) - length(unique(temp$CGID)), "Clinicians dropped due to no attending data in notes.\n")

## 4 Clinicians dropped due to no attending data in notes.

cat(length(unique(dat$SUBJECT_ID)) - length(unique(temp$SUBJECT_ID)), "Patients dropped due to no attending data in notes.\n")

## 78 Patients dropped due to no attending data in notes.

cat(length(unique(dat$HADM_ID)) - length(unique(temp$HADM_ID)), "Hospital Admissions dropped due to no attending data in notes.\n")

## 95 Hospital Admissions dropped due to no attending data in notes.

plotDat(temp, "NQF", "GENDER", F, "Gender", "Gender", "Frequency")

test <- table(temp$GENDER, temp$NQF)
test

##    
##        0    1
##   F 1338 4193
##   M 1812 3816

chisq.test(test)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 87.841, df = 1, p-value < 2.2e-16

pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = F, gtest = F, chisq = T, method = "fdr")

##   Comparison  p.Chisq p.adj.Chisq
## 1      F : M 7.09e-21    7.09e-21

plotDat(temp, "NQF","ETHNICITY", F, "Ethnicity", "Ethnicity", "Frequency")

test <- table(temp$ETHNICITY, temp$NQF)
test

##        
##            0    1
##   BLACK  250  774
##   OTHER  255  856
##   WHITE 2645 6379

chisq.test(test)

## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 27.836, df = 2, p-value = 9.024e-07

pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = F, gtest = F, chisq = T, method = "fdr")

##      Comparison  p.Chisq p.adj.Chisq
## 1 BLACK : OTHER 4.57e-01    4.57e-01
## 2 BLACK : WHITE 1.19e-03    1.78e-03
## 3 OTHER : WHITE 1.14e-05    3.42e-05

plotDat(temp, "NQF", "MARITAL_STATUS", F, "Marital Status", "Marital Status", "Frequency")

test <- table(temp$MARITAL_STATUS, temp$NQF)
test

##          
##              0    1
##   MARRIED 1509 3149
##   SINGLE   490 1863
##   UNKNOWN  175  239
##   WIDOWED  976 2758

chisq.test(test)

## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 151.94, df = 3, p-value < 2.2e-16

pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = F, gtest = F, chisq = T, method = "fdr")

##          Comparison  p.Chisq p.adj.Chisq
## 1  MARRIED : SINGLE 5.20e-24    3.12e-23
## 2 MARRIED : UNKNOWN 5.48e-05    5.48e-05
## 3 MARRIED : WIDOWED 5.10e-10    7.65e-10
## 4  SINGLE : UNKNOWN 8.35e-21    2.50e-20
## 5  SINGLE : WIDOWED 2.72e-06    3.26e-06
## 6 UNKNOWN : WIDOWED 5.29e-12    1.06e-11

plotDat(temp, "NQF", "FIRST_CAREUNIT", F, "First Careunit", "First Careunit", "Frequency")

test <- table(temp$FIRST_CAREUNIT, temp$NQF)
test

##        
##            0    1
##   CCU    406 1195
##   CSRU   211  173
##   MICU  1749 5594
##   SICU   494  645
##   TSICU  290  402

chisq.test(test)

## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 405.13, df = 4, p-value < 2.2e-16

pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = F, gtest = F, chisq = T, method = "fdr")

##      Comparison  p.Chisq p.adj.Chisq
## 1    CCU : CSRU 4.60e-29    1.53e-28
## 2    CCU : MICU 2.03e-01    2.26e-01
## 3    CCU : SICU 6.67e-23    1.33e-22
## 4   CCU : TSICU 3.79e-15    6.32e-15
## 5   CSRU : MICU 3.64e-42    1.82e-41
## 6   CSRU : SICU 1.06e-04    1.32e-04
## 7  CSRU : TSICU 5.24e-05    7.49e-05
## 8   MICU : SICU 7.77e-44    7.77e-43
## 9  MICU : TSICU 2.29e-25    5.72e-25
## 10 SICU : TSICU 5.72e-01    5.72e-01

boxplot(temp$AGE ~ temp$NQF, 
        main = "Caremeasure Implementation by Age",
        xlab = "Implementation (1 == Yes)",
        ylab = "Age (Years)")

t.test(temp$AGE ~ temp$NQF)

## 
##  Welch Two Sample t-test
## 
## data:  temp$AGE by temp$NQF
## t = -11.503, df = 6320.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.4056701 -0.9963286
## sample estimates:
## mean in group 0 mean in group 1 
##        83.23003        84.43103

boxplot(temp$SOFA ~ temp$NQF, 
        main = "Caremeasure Implementation by SOFA Score",
        xlab = "Implementation (1 == Yes)",
        ylab = "SOFA Score")

t.test(temp$SOFA ~ temp$NQF)

## 
##  Welch Two Sample t-test
## 
## data:  temp$SOFA by temp$NQF
## t = -4.5397, df = 5777.5, p-value = 5.748e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4092173 -0.1623825
## sample estimates:
## mean in group 0 mean in group 1 
##        4.828571        5.114371

detach("package:PresenceAbsence")

Data Standardization and Factoring (for modeling)

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

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

In this way, we can analyze the data in units of their deviation from the mean.

temp$AGE <- (temp$AGE - mean(temp$AGE))/sd(temp$AGE)
temp$SOFA <- (temp$SOFA - mean(temp$SOFA)/sd(temp$SOFA))

## Factor
temp <- within(temp, {
    ETHNICITY <- factor(ETHNICITY)
    GENDER <- factor(GENDER)
    MARITAL_STATUS <- factor(MARITAL_STATUS)
    FIRST_CAREUNIT <- factor(FIRST_CAREUNIT)
    SUBJECT_ID <- factor(SUBJECT_ID)
    CGID <- factor(CGID)
})

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 
##  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.379369

# ROC
rocplot(as.numeric(predict(gm_initial, type="response")), temp$NQF, col="blue", main = "Initial Model")

Model 1a (Patient Demographics)

Gender
Age
Ethnicity
Marrital Status
Care Provider (Random)

gm_a <- glmer(NQF ~ GENDER +
                     AGE +
                     ETHNICITY +
                     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 + MARITAL_STATUS + (1 | CGID)
##    Data: temp
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  11663.7  11729.6  -5822.9  11645.7    11150 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3725 -0.6454  0.3750  0.5813  2.7621 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CGID   (Intercept) 1.937    1.392   
## Number of obs: 11159, groups:  CGID, 493
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.06543    0.12256   8.693  < 2e-16 ***
## GENDERM               -0.36430    0.05429  -6.710 1.95e-11 ***
## AGE                    0.25305    0.02607   9.706  < 2e-16 ***
## ETHNICITYOTHER         0.48539    0.12337   3.934 8.34e-05 ***
## ETHNICITYWHITE         0.01746    0.09065   0.193  0.84727    
## MARITAL_STATUSSINGLE   0.50199    0.07124   7.047 1.83e-12 ***
## MARITAL_STATUSUNKNOWN -0.38472    0.12847  -2.995  0.00275 ** 
## MARITAL_STATUSWIDOWED  0.14126    0.06340   2.228  0.02588 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                 (Intr) GENDER AGE    ETHNICITYO ETHNICITYW MARITAL_STATUSS
## GENDERM         -0.275                                                    
## AGE              0.108 -0.061                                             
## ETHNICITYOT     -0.507 -0.050 -0.030                                      
## ETHNICITYWH     -0.676 -0.078 -0.067  0.665                               
## MARITAL_STATUSS -0.289  0.226  0.034  0.093      0.087                    
## MARITAL_STATUSU -0.133  0.068 -0.078 -0.081      0.053      0.150         
## MARITAL_STATUSW -0.373  0.361 -0.176  0.102      0.095      0.390         
##                 MARITAL_STATUSU
## GENDERM                        
## AGE                            
## ETHNICITYOT                    
## ETHNICITYWH                    
## MARITAL_STATUSS                
## MARITAL_STATUSU                
## MARITAL_STATUSW  0.202

sjt.glmer(gm_a)

	NQF
	Odds Ratio	CI	p
Fixed Parts
(Intercept)	2.90	2.28 – 3.69	<.001
GENDER (M)	0.69	0.62 – 0.77	<.001
AGE	1.29	1.22 – 1.36	<.001
ETHNICITY (OTHER)	1.62	1.28 – 2.07	<.001
ETHNICITY (WHITE)	1.02	0.85 – 1.22	.847
MARITAL_STATUS (SINGLE)	1.65	1.44 – 1.90	<.001
MARITAL_STATUS (UNKNOWN)	0.68	0.53 – 0.88	.003
MARITAL_STATUS (WIDOWED)	1.15	1.02 – 1.30	.026
Random Parts
τ_{00, CGID}	1.937
N_CGID	493
ICC_CGID	0.371
Observations	11159
Deviance	10536.307

# ROC
rocplot(as.numeric(predict(gm_a, type="response")), temp$NQF, col="blue", main = "Model 1a (Demographics)")

Model 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.334

sjt.glmer(gm_b)

	NQF
	Odds Ratio	CI	p
Fixed Parts
(Intercept)	2.62	2.24 – 3.07	<.001
SOFA	1.03	1.01 – 1.05	.001
Random Parts
τ_{00, CGID}	2.006
N_CGID	493
ICC_CGID	0.379
Observations	11159
Deviance	10761.908

# ROC
rocplot(as.numeric(predict(gm_b, type="response")), temp$NQF, col="blue", main = "Model 1b (Clinical Variables)")

Model 2a (Demographics + ICU)

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

gm_2a <- glmer(NQF ~ GENDER +
                      AGE +
                      ETHNICITY +
                      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 + MARITAL_STATUS + FIRST_CAREUNIT +  
##     (1 | CGID)
##    Data: temp
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  11620.3  11715.5  -5797.1  11594.3    11146 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3780 -0.6342  0.3752  0.5795  3.1486 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CGID   (Intercept) 1.8      1.342   
## Number of obs: 11159, groups:  CGID, 493
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.97904    0.14253   6.869 6.47e-12 ***
## GENDERM               -0.34575    0.05451  -6.343 2.25e-10 ***
## AGE                    0.25758    0.02615   9.849  < 2e-16 ***
## ETHNICITYOTHER         0.53361    0.12383   4.309 1.64e-05 ***
## ETHNICITYWHITE         0.03406    0.09052   0.376  0.70675    
## MARITAL_STATUSSINGLE   0.51937    0.07140   7.275 3.48e-13 ***
## MARITAL_STATUSUNKNOWN -0.40106    0.12901  -3.109  0.00188 ** 
## MARITAL_STATUSWIDOWED  0.17425    0.06371   2.735  0.00624 ** 
## FIRST_CAREUNITCSRU    -0.76667    0.17108  -4.481 7.42e-06 ***
## FIRST_CAREUNITMICU     0.13239    0.08679   1.525  0.12718    
## FIRST_CAREUNITSICU    -0.21185    0.12372  -1.712  0.08683 .  
## FIRST_CAREUNITTSICU    0.32709    0.14833   2.205  0.02745 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                 (Intr) GENDER AGE    ETHNICITYO ETHNICITYW MARITAL_STATUSS
## GENDERM         -0.254                                                    
## AGE              0.071 -0.058                                             
## ETHNICITYOT     -0.466 -0.044 -0.020                                      
## ETHNICITYWH     -0.603 -0.076 -0.066  0.661                               
## MARITAL_STATUSS -0.262  0.226  0.036  0.099      0.091                    
## MARITAL_STATUSU -0.086  0.067 -0.078 -0.081      0.048      0.147         
## MARITAL_STATUSW -0.315  0.361 -0.176  0.105      0.098      0.390         
## FIRST_CAREUNITC -0.231 -0.043 -0.003  0.002      0.016     -0.010         
## FIRST_CAREUNITM -0.500  0.037  0.056  0.068      0.033      0.023         
## FIRST_CAREUNITS -0.391  0.034 -0.006  0.021      0.038      0.007         
## FIRST_CAREUNITT -0.364  0.009 -0.001  0.052      0.068      0.032         
##                 MARITAL_STATUSU MARITAL_STATUSW FIRST_CAREUNITC
## GENDERM                                                        
## AGE                                                            
## ETHNICITYOT                                                    
## ETHNICITYWH                                                    
## MARITAL_STATUSS                                                
## MARITAL_STATUSU                                                
## MARITAL_STATUSW  0.200                                         
## FIRST_CAREUNITC -0.037          -0.061                         
## FIRST_CAREUNITM -0.035          -0.009           0.351         
## FIRST_CAREUNITS -0.047          -0.025           0.350         
## FIRST_CAREUNITT -0.083           0.020           0.320         
##                 FIRST_CAREUNITM FIRST_CAREUNITS
## GENDERM                                        
## AGE                                            
## ETHNICITYOT                                    
## ETHNICITYWH                                    
## MARITAL_STATUSS                                
## MARITAL_STATUSU                                
## MARITAL_STATUSW                                
## FIRST_CAREUNITC                                
## FIRST_CAREUNITM                                
## FIRST_CAREUNITS  0.559                         
## FIRST_CAREUNITT  0.467           0.592

sjt.glmer(gm_2a)

	NQF
	Odds Ratio	CI	p
Fixed Parts
(Intercept)	2.66	2.01 – 3.52	<.001
GENDER (M)	0.71	0.64 – 0.79	<.001
AGE	1.29	1.23 – 1.36	<.001
ETHNICITY (OTHER)	1.71	1.34 – 2.17	<.001
ETHNICITY (WHITE)	1.03	0.87 – 1.24	.707
MARITAL_STATUS (SINGLE)	1.68	1.46 – 1.93	<.001
MARITAL_STATUS (UNKNOWN)	0.67	0.52 – 0.86	.002
MARITAL_STATUS (WIDOWED)	1.19	1.05 – 1.35	.006
FIRST_CAREUNIT (CSRU)	0.46	0.33 – 0.65	<.001
FIRST_CAREUNIT (MICU)	1.14	0.96 – 1.35	.127
FIRST_CAREUNIT (SICU)	0.81	0.63 – 1.03	.087
FIRST_CAREUNIT (TSICU)	1.39	1.04 – 1.85	.027
Random Parts
τ_{00, CGID}	1.800
N_CGID	493
ICC_CGID	0.354
Observations	11159
Deviance	10513.905

# ROC
rocplot(as.numeric(predict(gm_2a, type="response")), temp$NQF, col="blue", main = "Model 2a (Demographics + ICU)")

Model 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.596

sjt.glmer(gm_2b)

	NQF
	Odds Ratio	CI	p
Fixed Parts
(Intercept)	2.58	2.10 – 3.18	<.001
SOFA	1.03	1.01 – 1.05	<.001
FIRST_CAREUNIT (CSRU)	0.45	0.32 – 0.62	<.001
FIRST_CAREUNIT (MICU)	1.07	0.91 – 1.27	.426
FIRST_CAREUNIT (SICU)	0.87	0.69 – 1.11	.266
FIRST_CAREUNIT (TSICU)	1.36	1.02 – 1.81	.035
Random Parts
τ_{00, CGID}	1.913
N_CGID	493
ICC_CGID	0.368
Observations	11159
Deviance	10736.210

# ROC
rocplot(as.numeric(predict(gm_2b, type="response")), temp$NQF, col="blue", main = "Model 2b (Clinical Variables + ICU)")

Model 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 +
                   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 + MARITAL_STATUS + SOFA + FIRST_CAREUNIT +  
##     (1 | CGID)
##    Data: temp
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  11595.5  11698.0  -5783.7  11567.5    11145 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.4330 -0.6264  0.3722  0.5762  3.4517 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CGID   (Intercept) 1.812    1.346   
## Number of obs: 11159, groups:  CGID, 493
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.825834   0.145765   5.666 1.47e-08 ***
## GENDERM               -0.374802   0.054877  -6.830 8.50e-12 ***
## AGE                    0.265987   0.026273  10.124  < 2e-16 ***
## ETHNICITYOTHER         0.531797   0.124263   4.280 1.87e-05 ***
## ETHNICITYWHITE         0.048033   0.090726   0.529  0.59651    
## MARITAL_STATUSSINGLE   0.531716   0.071543   7.432 1.07e-13 ***
## MARITAL_STATUSUNKNOWN -0.377833   0.129057  -2.928  0.00342 ** 
## MARITAL_STATUSWIDOWED  0.184446   0.063817   2.890  0.00385 ** 
## SOFA                   0.045084   0.008769   5.141 2.73e-07 ***
## FIRST_CAREUNITCSRU    -0.782858   0.171328  -4.569 4.89e-06 ***
## FIRST_CAREUNITMICU     0.128886   0.086892   1.483  0.13800    
## FIRST_CAREUNITSICU    -0.179204   0.124222  -1.443  0.14913    
## FIRST_CAREUNITTSICU    0.361549   0.148659   2.432  0.01501 *  
## ---
## 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

sjt.glmer(gm_icu)

	NQF
	Odds Ratio	CI	p
Fixed Parts
(Intercept)	2.28	1.72 – 3.04	<.001
GENDER (M)	0.69	0.62 – 0.77	<.001
AGE	1.30	1.24 – 1.37	<.001
ETHNICITY (OTHER)	1.70	1.33 – 2.17	<.001
ETHNICITY (WHITE)	1.05	0.88 – 1.25	.597
MARITAL_STATUS (SINGLE)	1.70	1.48 – 1.96	<.001
MARITAL_STATUS (UNKNOWN)	0.69	0.53 – 0.88	.003
MARITAL_STATUS (WIDOWED)	1.20	1.06 – 1.36	.004
SOFA	1.05	1.03 – 1.06	<.001
FIRST_CAREUNIT (CSRU)	0.46	0.33 – 0.64	<.001
FIRST_CAREUNIT (MICU)	1.14	0.96 – 1.35	.138
FIRST_CAREUNIT (SICU)	0.84	0.66 – 1.07	.149
FIRST_CAREUNIT (TSICU)	1.44	1.07 – 1.92	.015
Random Parts
τ_{00, CGID}	1.812
N_CGID	493
ICC_CGID	0.355
Observations	11159
Deviance	10485.611

# ROC
rocplot(as.numeric(predict(gm_icu, type="response")), temp$NQF, col="blue", main = "Model 3 (All Covariates)")

References

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.
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.
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.

#layout(mat = matrix(c(1,1,2,2,3,3,
#                      0,4,4,5,5,0), nrow = 2, byrow = TRUE))

#layout(matrix(c(1,2,3,4,4,5), 2, 3, byrow = TRUE))
# ROC

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)")

Care Measure Variation

E.T. Moseley

7/16/2018

Background

Model 0 (No Covariates)

Model 1a (Patient Demographics)

Model 1b (Clinical Variables)

Model 2a (Demographics + ICU)

Model 2b (Clinical Variables + ICU)

Model 3 (Demographics + Clinical Variables + ICU)

Libraries

Utility Functions

caremeasure_rate_check

caremeasure_check

plotDat

detach_packages

zero_range

t_tests

rocplot

Load Data and Merge to Care Providers Table

Load Severity of Illness Data

Data Cleaning

Check Attending and Caremeasure Implementation

Data Exploration

Data Standardization and Factoring (for modeling)

Baseline Model (Provider Random Effects)

Model 1a (Patient Demographics)

Model 1b (Clinical Variables)

Model 2a (Demographics + ICU)

Model 2b (Clinical Variables + ICU)

Model 3 (Demographics + Clinical Variables + ICU)

Generalized Linear Models

Data Standardization and Factoring (for modeling)

Baseline Model (Provider Random Effects)

Model 1a (Patient Demographics)

Model 1b (Clinical Variables)

Model 2a (Demographics + ICU)

Model 2b (Clinical Variables + ICU)

Model 3 (Demographics + Clinical Variables + ICU)

References