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.

Model 0 (No Covariates)

  1. Care Provider (Random)

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Marrital Status
  5. 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)
}
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.5,0), fill = cm.colors(length(rownames(prop))), legend=rownames(prop))
}
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 <- 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_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"

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 <- provider_caremeasure_check(attending_check(dat))

temp <- expire_check(temp)

Data Exploration

plotDat(temp, "NQF", "GENDER", F, "Gender (Provider-Level)", "Gender", "Frequency")

test <- table(temp$GENDER, temp$NQF)
test
##    
##        0    1
##   F 2187 3344
##   M 2768 2860
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 104.66, 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 1.45e-24    1.45e-24
plotDat(temp, "NQF", "ADMISSION_TYPE", F, "Admission Type (Provider-Level)", "Admission Type", "Frequency")

test <- table(temp$ADMISSION_TYPE, temp$NQF)
test
##            
##                0    1
##   ELECTIVE   294   28
##   EMERGENCY 4661 6176
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 293.48, 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 ELECTIVE : EMERGENCY 8.65e-66    8.65e-66
plotDat(temp, "NQF","ETHNICITY", F, "Ethnicity (Provider-Level)", "Ethnicity", "Frequency")

test <- table(temp$ETHNICITY, temp$NQF)
test
##        
##            0    1
##   BLACK  405  619
##   OTHER  475  636
##   WHITE 4075 4949
chisq.test(test)
## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 13.069, df = 2, p-value = 0.001452
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = F, gtest = F, chisq = T, method = "fdr")
##      Comparison  p.Chisq p.adj.Chisq
## 1 BLACK : OTHER 0.145000     0.14500
## 2 BLACK : WHITE 0.000706     0.00212
## 3 OTHER : WHITE 0.137000     0.14500
plotDat(temp, "NQF", "MARITAL_STATUS", F, "Marital Status (Provider-Level)", "Marital Status", "Frequency")

test <- table(temp$MARITAL_STATUS, temp$NQF)
test
##          
##              0    1
##   MARRIED 2315 2343
##   SINGLE   884 1469
##   UNKNOWN  227  187
##   WIDOWED 1529 2205
chisq.test(test)
## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 133.74, 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 7.70e-22    4.62e-21
## 2 MARRIED : UNKNOWN 5.12e-02    5.12e-02
## 3 MARRIED : WIDOWED 1.53e-15    4.59e-15
## 4  SINGLE : UNKNOWN 5.65e-11    1.13e-10
## 5  SINGLE : WIDOWED 9.39e-03    1.13e-02
## 6 UNKNOWN : WIDOWED 7.79e-08    1.17e-07
x <- barplot(table(mtcars$cyl), xaxt="n")
labs <- paste(names(table(mtcars$cyl)), "cylinders")
text(cex=1, x=x-.25, y=-1.25, labs, xpd=TRUE, srt=45)

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

test <- table(temp$FIRST_CAREUNIT, temp$NQF)
test
##        
##            0    1
##   CCU    696  905
##   CSRU   273  111
##   MICU  3031 4312
##   SICU   593  546
##   TSICU  362  330
chisq.test(test)
## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 185.04, 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.13e-22    2.06e-21
## 2    CCU : MICU 1.13e-01    1.26e-01
## 3    CCU : SICU 1.08e-05    1.54e-05
## 4   CCU : TSICU 1.17e-04    1.46e-04
## 5   CSRU : MICU 2.10e-30    2.10e-29
## 6   CSRU : SICU 1.10e-10    2.75e-10
## 7  CSRU : TSICU 2.91e-09    5.82e-09
## 8   MICU : SICU 9.46e-12    3.15e-11
## 9  MICU : TSICU 2.43e-08    4.05e-08
## 10 SICU : TSICU 9.56e-01    9.56e-01
boxplot(temp$AGE ~ temp$NQF, 
        main = "Caremeasure Implementation\n by Age (Provider-Level)",
        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 = -12.58, df = 10825, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.425367 -1.041058
## sample estimates:
## mean in group 0 mean in group 1 
##        83.40638        84.63959
boxplot(temp$SOFA ~ temp$NQF, 
        main = "Caremeasure Implementation\n by SOFA Score (Provider-Level)",
        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 = 0.11386, df = 10337, p-value = 0.9094
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1062057  0.1193046
## sample estimates:
## mean in group 0 mean in group 1 
##        5.037336        5.030787

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)

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

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Marrital Status
  5. 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 
##  13472.3  13538.1  -6727.1  13454.3    11150 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4389 -0.7877  0.3476  0.7307  4.1859 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CGID   (Intercept) 2.112    1.453   
## Number of obs: 11159, groups:  CGID, 493
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.22100    0.11608   1.904  0.05692 .  
## GENDERM               -0.33510    0.04934  -6.792 1.11e-11 ***
## AGE                    0.25396    0.02348  10.817  < 2e-16 ***
## ETHNICITYOTHER         0.03976    0.10762   0.369  0.71182    
## ETHNICITYWHITE        -0.02416    0.08113  -0.298  0.76585    
## MARITAL_STATUSSINGLE   0.38954    0.06319   6.164 7.08e-10 ***
## MARITAL_STATUSUNKNOWN -0.14376    0.12465  -1.153  0.24878    
## MARITAL_STATUSWIDOWED  0.17954    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 MARITAL_STATUSS
## GENDERM         -0.266                                                    
## AGE              0.089 -0.054                                             
## ETHNICITYOT     -0.495 -0.037 -0.051                                      
## ETHNICITYWH     -0.647 -0.064 -0.064  0.675                               
## MARITAL_STATUSS -0.302  0.248  0.020  0.093      0.095                    
## MARITAL_STATUSU -0.130  0.066 -0.077 -0.068      0.048      0.161         
## MARITAL_STATUSW -0.369  0.368 -0.177  0.101      0.096      0.432         
##                 MARITAL_STATUSU
## GENDERM                        
## AGE                            
## ETHNICITYOT                    
## ETHNICITYWH                    
## MARITAL_STATUSS                
## MARITAL_STATUSU                
## MARITAL_STATUSW  0.205
sjt.glmer(gm_a)
    NQF
    Odds Ratio CI p
Fixed Parts
(Intercept)   1.25 0.99 – 1.57 .057
GENDER (M)   0.72 0.65 – 0.79 <.001
AGE   1.29 1.23 – 1.35 <.001
ETHNICITY (OTHER)   1.04 0.84 – 1.28 .712
ETHNICITY (WHITE)   0.98 0.83 – 1.14 .766
MARITAL_STATUS (SINGLE)   1.48 1.30 – 1.67 <.001
MARITAL_STATUS (UNKNOWN)   0.87 0.68 – 1.11 .249
MARITAL_STATUS (WIDOWED)   1.20 1.07 – 1.34 .002
Random Parts
τ00, CGID   2.112
NCGID   493
ICCCGID   0.391
Observations   11159
Deviance   12244.947

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
sjt.glmer(gm_b)
    NQF
    Odds Ratio CI p
Fixed Parts
(Intercept)   1.18 1.01 – 1.38 .042
SOFA   1.00 0.99 – 1.02 .971
Random Parts
τ00, CGID   2.138
NCGID   493
ICCCGID   0.394
Observations   11159
Deviance   12489.152

Model 2a (Demographics + ICU)

  1. Demographics (Model 1a)
  2. First Intensive Care Unit
  3. 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 
##  13425.1  13520.2  -6699.5  13399.1    11146 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3984 -0.7883  0.3440  0.7329  4.7920 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CGID   (Intercept) 2.184    1.478   
## Number of obs: 11159, groups:  CGID, 493
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -0.089391   0.136182  -0.656 0.511562    
## GENDERM               -0.317199   0.049573  -6.399 1.57e-10 ***
## AGE                    0.259135   0.023635  10.964  < 2e-16 ***
## ETHNICITYOTHER         0.084372   0.108248   0.779 0.435723    
## ETHNICITYWHITE         0.006944   0.081463   0.085 0.932066    
## MARITAL_STATUSSINGLE   0.402595   0.063371   6.353 2.11e-10 ***
## MARITAL_STATUSUNKNOWN -0.189578   0.125578  -1.510 0.131135    
## MARITAL_STATUSWIDOWED  0.199230   0.058503   3.405 0.000660 ***
## FIRST_CAREUNITCSRU    -0.453809   0.171432  -2.647 0.008117 ** 
## FIRST_CAREUNITMICU     0.301187   0.078861   3.819 0.000134 ***
## FIRST_CAREUNITSICU     0.347727   0.116786   2.977 0.002906 ** 
## FIRST_CAREUNITTSICU    0.764958   0.142851   5.355 8.56e-08 ***
## ---
## 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.242                                                    
## AGE              0.044 -0.050                                             
## ETHNICITYOT     -0.455 -0.033 -0.045                                      
## ETHNICITYWH     -0.577 -0.061 -0.062  0.673                               
## MARITAL_STATUSS -0.272  0.246  0.022  0.096      0.097                    
## MARITAL_STATUSU -0.083  0.064 -0.078 -0.069      0.043      0.158         
## MARITAL_STATUSW -0.310  0.367 -0.176  0.103      0.098      0.431         
## FIRST_CAREUNITC -0.205 -0.031  0.002  0.011      0.014      0.000         
## FIRST_CAREUNITM -0.475  0.031  0.083  0.070      0.035      0.027         
## FIRST_CAREUNITS -0.386  0.033  0.014  0.026      0.043      0.015         
## FIRST_CAREUNITT -0.359  0.009  0.011  0.045      0.066      0.027         
##                 MARITAL_STATUSU MARITAL_STATUSW FIRST_CAREUNITC
## GENDERM                                                        
## AGE                                                            
## ETHNICITYOT                                                    
## ETHNICITYWH                                                    
## MARITAL_STATUSS                                                
## MARITAL_STATUSU                                                
## MARITAL_STATUSW  0.202                                         
## FIRST_CAREUNITC -0.031          -0.049                         
## FIRST_CAREUNITM -0.033          -0.010           0.300         
## FIRST_CAREUNITS -0.041          -0.018           0.290         
## FIRST_CAREUNITT -0.076           0.021           0.266         
##                 FIRST_CAREUNITM FIRST_CAREUNITS
## GENDERM                                        
## AGE                                            
## ETHNICITYOT                                    
## ETHNICITYWH                                    
## MARITAL_STATUSS                                
## MARITAL_STATUSU                                
## MARITAL_STATUSW                                
## FIRST_CAREUNITC                                
## FIRST_CAREUNITM                                
## FIRST_CAREUNITS  0.540                         
## FIRST_CAREUNITT  0.441           0.556
sjt.glmer(gm_2a)
    NQF
    Odds Ratio CI p
Fixed Parts
(Intercept)   0.91 0.70 – 1.19 .512
GENDER (M)   0.73 0.66 – 0.80 <.001
AGE   1.30 1.24 – 1.36 <.001
ETHNICITY (OTHER)   1.09 0.88 – 1.35 .436
ETHNICITY (WHITE)   1.01 0.86 – 1.18 .932
MARITAL_STATUS (SINGLE)   1.50 1.32 – 1.69 <.001
MARITAL_STATUS (UNKNOWN)   0.83 0.65 – 1.06 .131
MARITAL_STATUS (WIDOWED)   1.22 1.09 – 1.37 <.001
FIRST_CAREUNIT (CSRU)   0.64 0.45 – 0.89 .008
FIRST_CAREUNIT (MICU)   1.35 1.16 – 1.58 <.001
FIRST_CAREUNIT (SICU)   1.42 1.13 – 1.78 .003
FIRST_CAREUNIT (TSICU)   2.15 1.62 – 2.84 <.001
Random Parts
τ00, CGID   2.184
NCGID   493
ICCCGID   0.399
Observations   11159
Deviance   12180.260

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
sjt.glmer(gm_2b)
    NQF
    Odds Ratio CI p
Fixed Parts
(Intercept)   0.93 0.75 – 1.13 .455
SOFA   1.00 0.99 – 1.02 .676
FIRST_CAREUNIT (CSRU)   0.63 0.45 – 0.87 .006
FIRST_CAREUNIT (MICU)   1.27 1.09 – 1.48 .002
FIRST_CAREUNIT (SICU)   1.47 1.17 – 1.84 <.001
FIRST_CAREUNIT (TSICU)   2.09 1.58 – 2.75 <.001
Random Parts
τ00, CGID   2.243
NCGID   493
ICCCGID   0.405
Observations   11159
Deviance   12421.574

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 +
                   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 
##  13420.5  13523.0  -6696.3  13392.5    11145 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2806 -0.7926  0.3433  0.7313  4.9852 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CGID   (Intercept) 2.194    1.481   
## Number of obs: 11159, groups:  CGID, 493
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -0.158075   0.138990  -1.137 0.255410    
## GENDERM               -0.331264   0.049905  -6.638 3.18e-11 ***
## AGE                    0.263075   0.023705  11.098  < 2e-16 ***
## ETHNICITYOTHER         0.080285   0.108411   0.741 0.458960    
## ETHNICITYWHITE         0.014362   0.081614   0.176 0.860314    
## MARITAL_STATUSSINGLE   0.408120   0.063430   6.434 1.24e-10 ***
## MARITAL_STATUSUNKNOWN -0.179696   0.125682  -1.430 0.152785    
## MARITAL_STATUSWIDOWED  0.203522   0.058556   3.476 0.000510 ***
## SOFA                   0.020113   0.007848   2.563 0.010383 *  
## FIRST_CAREUNITCSRU    -0.462606   0.171437  -2.698 0.006967 ** 
## FIRST_CAREUNITMICU     0.297664   0.078919   3.772 0.000162 ***
## FIRST_CAREUNITSICU     0.362566   0.117035   3.098 0.001949 ** 
## FIRST_CAREUNITTSICU    0.782931   0.143072   5.472 4.44e-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
sjt.glmer(gm_icu)
    NQF
    Odds Ratio CI p
Fixed Parts
(Intercept)   0.85 0.65 – 1.12 .255
GENDER (M)   0.72 0.65 – 0.79 <.001
AGE   1.30 1.24 – 1.36 <.001
ETHNICITY (OTHER)   1.08 0.88 – 1.34 .459
ETHNICITY (WHITE)   1.01 0.86 – 1.19 .860
MARITAL_STATUS (SINGLE)   1.50 1.33 – 1.70 <.001
MARITAL_STATUS (UNKNOWN)   0.84 0.65 – 1.07 .153
MARITAL_STATUS (WIDOWED)   1.23 1.09 – 1.37 <.001
SOFA   1.02 1.00 – 1.04 .010
FIRST_CAREUNIT (CSRU)   0.63 0.45 – 0.88 .007
FIRST_CAREUNIT (MICU)   1.35 1.15 – 1.57 <.001
FIRST_CAREUNIT (SICU)   1.44 1.14 – 1.81 .002
FIRST_CAREUNIT (TSICU)   2.19 1.65 – 2.90 <.001
Random Parts
τ00, CGID   2.194
NCGID   493
ICCCGID   0.400
Observations   11159
Deviance   12172.161
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 Linear Regression Models

Baseline Model (Provider Random Effects)

Not meaningful.

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Marrital Status
  5. Care Provider (Random)
gm_a <- glm(NQF ~ GENDER +
                     AGE +
                     ETHNICITY +
                     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 + MARITAL_STATUS, 
##     family = binomial(link = "logit"), data = temp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6955  -1.2224   0.8935   1.0738   1.5704  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.41451    0.07676   5.400 6.67e-08 ***
## GENDERM               -0.31305    0.04191  -7.470 8.04e-14 ***
## AGE                    0.25069    0.02009  12.476  < 2e-16 ***
## ETHNICITYOTHER        -0.00460    0.09117  -0.050 0.959760    
## ETHNICITYWHITE        -0.16072    0.06906  -2.327 0.019959 *  
## MARITAL_STATUSSINGLE   0.40458    0.05420   7.464 8.37e-14 ***
## MARITAL_STATUSUNKNOWN -0.39663    0.10611  -3.738 0.000185 ***
## MARITAL_STATUSWIDOWED  0.10679    0.04917   2.172 0.029868 *  
## ---
## 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: 14980  on 11151  degrees of freedom
## AIC: 14996
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## $`OR Summary`
##                              OR     2.5 %    97.5 %
## (Intercept)           1.5136245 1.3026783 1.7601383
## GENDERM               0.7312139 0.6735309 0.7937934
## AGE                   1.2849074 1.2353570 1.3366010
## ETHNICITYOTHER        0.9954107 0.8324493 1.1901204
## ETHNICITYWHITE        0.8515325 0.7433614 0.9745490
## MARITAL_STATUSSINGLE  1.4986755 1.3478250 1.6669219
## MARITAL_STATUSUNKNOWN 0.6725828 0.5458913 0.8276527
## MARITAL_STATUSWIDOWED 1.1127017 1.0104664 1.2252801

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 +
                      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 + MARITAL_STATUS + 
##     FIRST_CAREUNIT, family = binomial(link = "logit"), data = temp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7632  -1.2201   0.8364   1.0560   1.7705  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.36198    0.09268   3.906 9.39e-05 ***
## GENDERM               -0.29432    0.04240  -6.941 3.90e-12 ***
## AGE                    0.25485    0.02026  12.577  < 2e-16 ***
## ETHNICITYOTHER         0.07258    0.09240   0.785  0.43217    
## ETHNICITYWHITE        -0.13145    0.06950  -1.891  0.05858 .  
## MARITAL_STATUSSINGLE   0.39380    0.05451   7.225 5.01e-13 ***
## MARITAL_STATUSUNKNOWN -0.29419    0.10809  -2.722  0.00649 ** 
## MARITAL_STATUSWIDOWED  0.15119    0.04982   3.034  0.00241 ** 
## FIRST_CAREUNITCSRU    -1.08414    0.12547  -8.641  < 2e-16 ***
## FIRST_CAREUNITMICU     0.13121    0.05697   2.303  0.02127 *  
## FIRST_CAREUNITSICU    -0.36046    0.07956  -4.531 5.88e-06 ***
## FIRST_CAREUNITTSICU   -0.28322    0.09372  -3.022  0.00251 ** 
## ---
## 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: 14804  on 11147  degrees of freedom
## AIC: 14828
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## $`OR Summary`
##                              OR     2.5 %    97.5 %
## (Intercept)           1.4361771 1.1980433 1.7229352
## GENDERM               0.7450390 0.6856020 0.8095905
## AGE                   1.2902742 1.2401053 1.3426299
## ETHNICITYOTHER        1.0752767 0.8970970 1.2887367
## ETHNICITYWHITE        0.8768212 0.7647764 1.0043495
## MARITAL_STATUSSINGLE  1.4826089 1.3325846 1.6500387
## MARITAL_STATUSUNKNOWN 0.7451372 0.6024862 0.9205865
## MARITAL_STATUSWIDOWED 1.1632135 1.0550064 1.2825623
## FIRST_CAREUNITCSRU    0.3381933 0.2636760 0.4313495
## FIRST_CAREUNITMICU    1.1402075 1.0195868 1.2747458
## FIRST_CAREUNITSICU    0.6973538 0.5965642 0.8149293
## FIRST_CAREUNITTSICU   0.7533563 0.6268141 0.9051693

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 +
              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 + MARITAL_STATUS + 
##     SOFA + FIRST_CAREUNIT, family = binomial(link = "logit"), 
##     data = temp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7742  -1.2190   0.8377   1.0578   1.7829  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.323744   0.095265   3.398 0.000678 ***
## GENDERM               -0.302932   0.042703  -7.094 1.30e-12 ***
## AGE                    0.257538   0.020330  12.668  < 2e-16 ***
## ETHNICITYOTHER         0.070420   0.092457   0.762 0.446265    
## ETHNICITYWHITE        -0.126631   0.069594  -1.820 0.068826 .  
## MARITAL_STATUSSINGLE   0.396581   0.054539   7.272 3.55e-13 ***
## MARITAL_STATUSUNKNOWN -0.289381   0.108153  -2.676 0.007458 ** 
## MARITAL_STATUSWIDOWED  0.153350   0.049840   3.077 0.002092 ** 
## SOFA                   0.011589   0.006655   1.741 0.081608 .  
## FIRST_CAREUNITCSRU    -1.091802   0.125550  -8.696  < 2e-16 ***
## FIRST_CAREUNITMICU     0.128070   0.057006   2.247 0.024664 *  
## FIRST_CAREUNITSICU    -0.355293   0.079631  -4.462 8.13e-06 ***
## FIRST_CAREUNITTSICU   -0.276171   0.093780  -2.945 0.003231 ** 
## ---
## 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: 14801  on 11146  degrees of freedom
## AIC: 14827
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## $`OR Summary`
##                              OR     2.5 %    97.5 %
## (Intercept)           1.3822940 1.1472187 1.6666589
## GENDERM               0.7386491 0.6793220 0.8031131
## AGE                   1.2937404 1.2432756 1.3464150
## ETHNICITYOTHER        1.0729589 0.8950591 1.2861032
## ETHNICITYWHITE        0.8810590 0.7683379 1.0093910
## MARITAL_STATUSSINGLE  1.4867326 1.3362079 1.6547356
## MARITAL_STATUSUNKNOWN 0.7487268 0.6053150 0.9251427
## MARITAL_STATUSWIDOWED 1.1657328 1.0572550 1.2853853
## SOFA                  1.0116560 0.9985574 1.0249514
## FIRST_CAREUNITCSRU    0.3356111 0.2616212 0.4281265
## FIRST_CAREUNITMICU    1.1366331 1.0163182 1.2708378
## FIRST_CAREUNITSICU    0.7009680 0.5995752 0.8192669
## FIRST_CAREUNITTSICU   0.7586832 0.6311746 0.9116708
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_caremeasure_check(attending_check(dat))

temp <- expire_check(temp)

Data Exploration

plotDat(temp, "NQF", "GENDER", F, "Gender (Care Team Level)", "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", "ADMISSION_TYPE", F, "Admission Type (Care Team Level)", "Admission Type", "Frequency")

test <- table(temp$ADMISSION_TYPE, temp$NQF)
test
##            
##                0    1
##   ELECTIVE   275   47
##   EMERGENCY 2875 7962
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 532.1, 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 ELECTIVE : EMERGENCY 9.89e-118   9.89e-118
plotDat(temp, "NQF","ETHNICITY", F, "Ethnicity (Care Team Level)", "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 (Care Team Level)", "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 (Care Team Level)", "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\n by Age (Care Team Level)",
        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\n by SOFA Score (Care Team Level)",
        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

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)

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

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Marrital Status
  5. 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
NCGID   493
ICCCGID   0.371
Observations   11159
Deviance   10536.307

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
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
NCGID   493
ICCCGID   0.379
Observations   11159
Deviance   10761.908

Model 2a (Demographics + ICU)

  1. Demographics (Model 1a)
  2. First Intensive Care Unit
  3. 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
NCGID   493
ICCCGID   0.354
Observations   11159
Deviance   10513.905

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
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
NCGID   493
ICCCGID   0.368
Observations   11159
Deviance   10736.210

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 +
                   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
NCGID   493
ICCCGID   0.355
Observations   11159
Deviance   10485.611
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 Linear Regression Models

Baseline Model (Provider Random Effects)

Not meaningful.

Model 1a (Patient Demographics)

  1. Gender
  2. Age
  3. Ethnicity
  4. Marrital Status
  5. Care Provider (Random)
gm_a <- glm(NQF ~ GENDER +
                     AGE +
                     ETHNICITY +
                     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 + MARITAL_STATUS, 
##     family = binomial(link = "logit"), data = temp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0524  -1.3558   0.7131   0.8514   1.3900  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.14273    0.08652  13.208  < 2e-16 ***
## GENDERM               -0.32540    0.04645  -7.006 2.46e-12 ***
## AGE                    0.25791    0.02235  11.542  < 2e-16 ***
## ETHNICITYOTHER         0.27144    0.10569   2.568   0.0102 *  
## ETHNICITYWHITE        -0.17328    0.07843  -2.209   0.0272 *  
## MARITAL_STATUSSINGLE   0.52593    0.06214   8.463  < 2e-16 ***
## MARITAL_STATUSUNKNOWN -0.68581    0.10836  -6.329 2.47e-10 ***
## MARITAL_STATUSWIDOWED  0.05847    0.05395   1.084   0.2785    
## ---
## 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 11151  degrees of freedom
## AIC: 12930
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## $`OR Summary`
##                              OR     2.5 %    97.5 %
## (Intercept)           3.1353184 2.6496285 3.7197749
## GENDERM               0.7222378 0.6593390 0.7910238
## AGE                   1.2942246 1.2388622 1.3522800
## ETHNICITYOTHER        1.3118546 1.0664176 1.6140752
## ETHNICITYWHITE        0.8409056 0.7199808 0.9792482
## MARITAL_STATUSSINGLE  1.6920322 1.4988336 1.9123353
## MARITAL_STATUSUNKNOWN 0.5036813 0.4075010 0.6233374
## MARITAL_STATUSWIDOWED 1.0602117 0.9538886 1.1785708

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 +
                      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 + MARITAL_STATUS + 
##     FIRST_CAREUNIT, family = binomial(link = "logit"), data = temp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1567  -1.1270   0.6710   0.8155   1.5558  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.17362    0.10565  11.108  < 2e-16 ***
## GENDERM               -0.31896    0.04750  -6.716 1.87e-11 ***
## AGE                    0.26826    0.02277  11.779  < 2e-16 ***
## ETHNICITYOTHER         0.41483    0.10855   3.822 0.000133 ***
## ETHNICITYWHITE        -0.13291    0.07975  -1.667 0.095600 .  
## MARITAL_STATUSSINGLE   0.50508    0.06306   8.010 1.15e-15 ***
## MARITAL_STATUSUNKNOWN -0.50087    0.11204  -4.470 7.81e-06 ***
## MARITAL_STATUSWIDOWED  0.11838    0.05528   2.141 0.032235 *  
## FIRST_CAREUNITCSRU    -1.18681    0.12034  -9.863  < 2e-16 ***
## FIRST_CAREUNITMICU     0.13976    0.06504   2.149 0.031650 *  
## FIRST_CAREUNITSICU    -0.84086    0.08503  -9.889  < 2e-16 ***
## FIRST_CAREUNITTSICU   -0.66598    0.09906  -6.723 1.78e-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: 12545  on 11147  degrees of freedom
## AIC: 12569
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## $`OR Summary`
##                              OR     2.5 %    97.5 %
## (Intercept)           3.2336696 2.6320563 3.9828340
## GENDERM               0.7269040 0.6622378 0.7977673
## AGE                   1.3076887 1.2507043 1.3675004
## ETHNICITYOTHER        1.5141176 1.2240603 1.8735154
## ETHNICITYWHITE        0.8755423 0.7477099 1.0222409
## MARITAL_STATUSSINGLE  1.6571136 1.4652722 1.8762073
## MARITAL_STATUSUNKNOWN 0.6060026 0.4868062 0.7554407
## MARITAL_STATUSWIDOWED 1.1256749 1.0101810 1.2546353
## FIRST_CAREUNITCSRU    0.3051929 0.2408667 0.3861290
## FIRST_CAREUNITMICU    1.1500031 1.0115715 1.3054206
## FIRST_CAREUNITSICU    0.4313385 0.3649816 0.5093978
## FIRST_CAREUNITTSICU   0.5137711 0.4231037 0.6239222

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 +
              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 + MARITAL_STATUS + 
##     SOFA + FIRST_CAREUNIT, family = binomial(link = "logit"), 
##     data = temp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2506  -1.1409   0.6676   0.8177   1.5492  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.037661   0.108429   9.570  < 2e-16 ***
## GENDERM               -0.348734   0.047881  -7.283 3.26e-13 ***
## AGE                    0.278813   0.022917  12.166  < 2e-16 ***
## ETHNICITYOTHER         0.411646   0.108942   3.779 0.000158 ***
## ETHNICITYWHITE        -0.118724   0.079960  -1.485 0.137599    
## MARITAL_STATUSSINGLE   0.517150   0.063203   8.182 2.78e-16 ***
## MARITAL_STATUSUNKNOWN -0.482863   0.112254  -4.302 1.70e-05 ***
## MARITAL_STATUSWIDOWED  0.127479   0.055333   2.304 0.021231 *  
## SOFA                   0.042354   0.007607   5.568 2.58e-08 ***
## FIRST_CAREUNITCSRU    -1.218994   0.120658 -10.103  < 2e-16 ***
## FIRST_CAREUNITMICU     0.129467   0.065185   1.986 0.047016 *  
## FIRST_CAREUNITSICU    -0.824054   0.085209  -9.671  < 2e-16 ***
## FIRST_CAREUNITTSICU   -0.644605   0.099133  -6.502 7.90e-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: 12514  on 11146  degrees of freedom
## AIC: 12540
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## $`OR Summary`
##                              OR     2.5 %    97.5 %
## (Intercept)           2.8226063 2.2847918 3.4951818
## GENDERM               0.7055805 0.6423184 0.7749432
## AGE                   1.3215604 1.2636228 1.3824017
## ETHNICITYOTHER        1.5093001 1.2192309 1.8690038
## ETHNICITYWHITE        0.8880526 0.7580906 1.0372807
## MARITAL_STATUSSINGLE  1.6772412 1.4826441 1.8995524
## MARITAL_STATUSUNKNOWN 0.6170144 0.4954494 0.7694974
## MARITAL_STATUSWIDOWED 1.1359613 1.0193093 1.2662317
## SOFA                  1.0432640 1.0278801 1.0589950
## FIRST_CAREUNITCSRU    0.2955274 0.2330900 0.3741354
## FIRST_CAREUNITMICU    1.1382211 1.0009305 1.2924047
## FIRST_CAREUNITSICU    0.4386497 0.3710425 0.5182143
## FIRST_CAREUNITTSICU   0.5248698 0.4321841 0.6374953
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)")

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.