Load Data

dat <- read.csv("~/nqf_caregivers/data/nqf_dat.csv", header = T, stringsAsFactors = F)

## Edit language
dat$LANGUAGE <- ifelse(dat$LANGUAGE == "ENGL", "ENGL", "OTHER")

team <- team_caremeasure_check(dat)

team <- expire_check(team)

## Load ICU stay details
details <- read.csv("~/nqf_caregivers/data/icustay_detail.csv", header = T, stringsAsFactors = F)
colnames(details) <- toupper(colnames(details))

## Clean details
details <- details[,c("SUBJECT_ID", "HADM_ID", "ICUSTAY_ID", "LOS_HOSPITAL", "LOS_ICU")]

## Merge details
team <- merge(team, details, by = c("SUBJECT_ID", "HADM_ID", "ICUSTAY_ID"))
## remove details
rm(details)

## Load ventdurations.csv
vent <- read.csv("~/nqf_caregivers/data/ventdurations.csv", header = T, stringsAsFactors = F)
colnames(vent) <- toupper(colnames(vent))

## Change colnames
colnames(vent) <- c("ICUSTAY_ID", "VENTNUM", "VENT_STARTTIME", "VENT_ENDTIME", "VENT_DURATION_HOURS")

## Clean vent for only those stays within our cohort
vent <- vent[(vent$ICUSTAY_ID %in% team$ICUSTAY_ID),]

## Find length until VENT_STARTTIME from their admission
vent <- merge(vent, team[,c("ICUSTAY_ID", "ADMITTIME")], by = "ICUSTAY_ID")

## Drop VENTNUM as it is unimportant
vent$VENTNUM <- NULL

## Remove duplicates from join
vent <- vent[!duplicated(vent),]

## Adjust dates
vent$VENT_STARTTIME <- as.numeric(as.Date(vent$VENT_STARTTIME, "%Y-%m-%d %H:%M:%S"))
vent$VENT_ENDTIME <- as.numeric(as.Date(vent$VENT_ENDTIME, "%Y-%m-%d %H:%M:%S"))

## Time until ventilation
vent$TIME_UNTIL_VENTILATION <- vent$VENT_STARTTIME - vent$ADMITTIME

## If time until ventilation is less than or equal to 2 (two days), give it a 1, else 0
vent$VENT_FIRST_48 <- ifelse(vent$TIME_UNTIL_VENTILATION <= 2, 1, 0)

## Get unique icustay ids with ventilation data
## ids <- unique(team[team$ICUSTAY_ID %in% vent$ICUSTAY_ID,]$ICUSTAY_ID)
ids <- unique(team$ICUSTAY_ID)
## Create a vector to hold results
MECHVENT <- vector()
for (id in ids){
    if (any(vent[(vent$ICUSTAY_ID == id), ]$VENT_FIRST_48 == 1)){
        MECHVENT <- c(MECHVENT, rep(1, each = nrow(team[(team$ICUSTAY_ID == id), ])))
    } else {
        MECHVENT <- c(MECHVENT, rep(0, each = nrow(team[(team$ICUSTAY_ID == id), ])))
    }
}

## Add to data frame
team$MECHVENT <- MECHVENT

## Clean up data
rm(MECHVENT)

Hospital Admission Level Analysis (Hospital Admission Level)

temp <- aggregate(cbind(NQF,
                        MECHVENT,
                        AGE,
                        SOFA,
                        LOS_HOSPITAL,
                        LOS_ICU) ~ 
                     SUBJECT_ID +
                     ETHNICITY +
                     GENDER +
                     MARITAL_STATUS +
                     FIRST_CAREUNIT +
                     HADM_ID +
                     HOSPITAL_EXPIRE_FLAG +
                     ADMISSION_TYPE +
                     EXPIRE_FLAG +
                     LANGUAGE,
                 data = team, 
                 FUN = mean)

nrow(temp)
## [1] 1350

Table 1 Cohort characteristics by the five ICUs

  1. In table 1 I think we are missing data for SOFA score and Age stratified by ICU.
  2. Is it possible to add ICU length of stay (in days) and if on mechanical ventilation at time of documentation (or within 48 hours of ICU admission)
  3. I like the set up for table two. I added a ghost section on Table 2 for the remainder of the analysis. Basically just comparing each covariate based on performance of NQF measure.
  4. In the logistic regression: Can you change the reference group to MICU and white race. Can you add mechanical ventilation at time of NQF documentation as a covariate? Can you add a cox regression for mortality. I have made a ghost table with the covariates in the google doc.
table(temp$GENDER)
## 
##   F   M 
## 677 673
table(temp$GENDER, temp$FIRST_CAREUNIT)
##    
##     CCU CSRU MICU SICU TSICU
##   F  82   27  385  116    67
##   M 101   47  372   89    64
summary(temp$AGE)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.01   79.97   83.93   84.17   88.40   91.40
sd(temp$AGE)
## [1] 5.177834
summary(temp$AGE[temp$FIRST_CAREUNIT == "CCU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.13   80.88   84.77   84.62   88.75   91.40
sd(temp$AGE[temp$FIRST_CAREUNIT == "CCU"])
## [1] 4.871286
summary(temp$AGE[temp$FIRST_CAREUNIT == "CSRU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.27   79.12   82.35   83.06   86.02   91.40
sd(temp$AGE[temp$FIRST_CAREUNIT == "CSRU"])
## [1] 5.069464
summary(temp$AGE[temp$FIRST_CAREUNIT == "MICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.04   79.87   84.12   84.11   88.38   91.40
sd(temp$AGE[temp$FIRST_CAREUNIT == "MICU"])
## [1] 5.217173
summary(temp$AGE[temp$FIRST_CAREUNIT == "SICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.02   79.70   83.64   84.29   88.99   91.40
sd(temp$AGE[temp$FIRST_CAREUNIT == "SICU"])
## [1] 5.349049
summary(temp$AGE[temp$FIRST_CAREUNIT == "TSICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.01   80.24   83.56   84.33   88.43   91.40
sd(temp$AGE[temp$FIRST_CAREUNIT == "TSICU"])
## [1] 5.129848
table(temp$MARITAL_STATUS)
## 
## MARRIED  SINGLE UNKNOWN WIDOWED 
##     561     270      57     462
table(temp$MARITAL_STATUS, temp$FIRST_CAREUNIT)
##          
##           CCU CSRU MICU SICU TSICU
##   MARRIED  85   32  308   81    55
##   SINGLE   39   11  162   36    22
##   UNKNOWN   4    4   21   14    14
##   WIDOWED  55   27  266   74    40
table(temp$ETHNICITY)
## 
## BLACK OTHER WHITE 
##   110   135  1105
table(temp$ETHNICITY, temp$FIRST_CAREUNIT)
##        
##         CCU CSRU MICU SICU TSICU
##   BLACK   9    4   77   11     9
##   OTHER  23   10   59   28    15
##   WHITE 151   60  621  166   107
table(temp$LANGUAGE)
## 
##  ENGL OTHER 
##  1167   183
table(temp$LANGUAGE, temp$FIRST_CAREUNIT)
##        
##         CCU CSRU MICU SICU TSICU
##   ENGL  159   68  637  185   118
##   OTHER  24    6  120   20    13
table(temp$ADMISSION_TYPE)
## 
##  ELECTIVE EMERGENCY 
##        60      1290
table(temp$ADMISSION_TYPE, temp$FIRST_CAREUNIT)
##            
##             CCU CSRU MICU SICU TSICU
##   ELECTIVE    1   21   17   13     8
##   EMERGENCY 182   53  740  192   123
summary(temp$SOFA)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    3.00    4.00    4.73    6.00   18.00
sd(temp$SOFA)
## [1] 2.913163
summary(temp$SOFA[temp$FIRST_CAREUNIT == "CCU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   4.000   4.705   6.000  18.000
sd(temp$SOFA[temp$FIRST_CAREUNIT == "CCU"])
## [1] 3.130904
summary(temp$SOFA[temp$FIRST_CAREUNIT == "CSRU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   3.000   5.000   5.324   7.000  15.000
sd(temp$SOFA[temp$FIRST_CAREUNIT == "CSRU"])
## [1] 2.819249
summary(temp$SOFA[temp$FIRST_CAREUNIT == "MICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   3.000   5.000   4.978   6.000  18.000
sd(temp$SOFA[temp$FIRST_CAREUNIT == "MICU"])
## [1] 2.91743
summary(temp$SOFA[temp$FIRST_CAREUNIT == "SICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   4.000   4.015   5.000  18.000
sd(temp$SOFA[temp$FIRST_CAREUNIT == "SICU"])
## [1] 2.599921
summary(temp$SOFA[temp$FIRST_CAREUNIT == "TSICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   4.000   4.115   5.000  15.000
sd(temp$SOFA[temp$FIRST_CAREUNIT == "TSICU"])
## [1] 2.832887
summary(temp$LOS_ICU)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0245  1.2696  2.1644  3.8971  4.1884 60.9390
sd(temp$LOS_ICU)
## [1] 5.161508
summary(temp$LOS_ICU[temp$FIRST_CAREUNIT == "CCU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4224  1.4459  2.4031  3.9177  4.7260 30.9020
sd(temp$LOS_ICU[temp$FIRST_CAREUNIT == "CCU"])
## [1] 4.202086
summary(temp$LOS_ICU[temp$FIRST_CAREUNIT == "CSRU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6406  1.7817  3.0964  5.2287  5.4189 43.1362
sd(temp$LOS_ICU[temp$FIRST_CAREUNIT == "CSRU"])
## [1] 6.863161
summary(temp$LOS_ICU[temp$FIRST_CAREUNIT == "MICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0453  1.1796  2.0308  3.7672  3.9230 60.9390
sd(temp$LOS_ICU[temp$FIRST_CAREUNIT == "MICU"])
## [1] 5.493012
summary(temp$LOS_ICU[temp$FIRST_CAREUNIT == "SICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5636  1.4198  2.2541  3.9900  4.9275 30.9994
sd(temp$LOS_ICU[temp$FIRST_CAREUNIT == "SICU"])
## [1] 4.681092
summary(temp$LOS_ICU[temp$FIRST_CAREUNIT == "TSICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0245  1.4347  2.3408  3.7212  4.4143 20.8696
sd(temp$LOS_ICU[temp$FIRST_CAREUNIT == "TSICU"])
## [1] 3.778244
summary(temp$LOS_HOSPITAL)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##   0.4479   3.8476   5.9858   8.2009   9.7984 105.5799
sd(temp$LOS_HOSPITAL)
## [1] 7.47117
summary(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "CCU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.800   3.515   5.521   7.143   8.633  34.727
sd(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "CCU"])
## [1] 5.893634
summary(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "CSRU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.760   4.979   7.182   8.950  10.711  43.302
sd(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "CSRU"])
## [1] 6.860786
summary(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "MICU"])
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##   0.4479   3.8417   5.9326   8.2398   9.7604 105.5799
sd(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "MICU"])
## [1] 7.99587
summary(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "SICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8444  3.9924  6.3958  8.3527 11.2104 32.8618
sd(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "SICU"])
## [1] 6.31487
summary(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "TSICU"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8792  4.0983  6.0507  8.7921 10.4639 56.7653
sd(temp$LOS_HOSPITAL[temp$FIRST_CAREUNIT == "TSICU"])
## [1] 8.21543
table(temp$NQF)
## 
##   0   1 
## 477 873
table(temp$NQF, temp$FIRST_CAREUNIT)
##    
##     CCU CSRU MICU SICU TSICU
##   0  53   52  190  113    69
##   1 130   22  567   92    62
table(temp$MECHVENT)
## 
##   0   1 
## 870 480
table(temp$MECHVENT, temp$FIRST_CAREUNIT)
##    
##     CCU CSRU MICU SICU TSICU
##   0 124   29  536  109    72
##   1  59   45  221   96    59
table(temp$HOSPITAL_EXPIRE_FLAG)
## 
##    0    1 
## 1075  275
table(temp$HOSPITAL_EXPIRE_FLAG, temp$FIRST_CAREUNIT)
##    
##     CCU CSRU MICU SICU TSICU
##   0 154   55  618  153    95
##   1  29   19  139   52    36

Table 2 Unadjusted analysis

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

test <- table(temp$GENDER, temp$NQF)
test
##    
##       0   1
##   F 207 470
##   M 270 403
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 13.037, df = 1, p-value = 0.0003055
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##   Comparison p.Fisher p.adj.Fisher
## 1      F : M 0.000267     0.000267
plotDat(temp, "NQF","ETHNICITY", F, "Ethnicity", "Ethnicity", "Frequency")

test <- table(temp$ETHNICITY, temp$NQF)
test
##        
##           0   1
##   BLACK  35  75
##   OTHER  41  94
##   WHITE 401 704
chisq.test(test)
## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 2.4924, df = 2, p-value = 0.2876
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##      Comparison p.Fisher p.adj.Fisher
## 1 BLACK : OTHER    0.890        0.890
## 2 BLACK : WHITE    0.405        0.608
## 3 OTHER : WHITE    0.184        0.552
plotDat(temp, "NQF","LANGUAGE", F, "Language", "Language", "Frequency")

test <- table(temp$LANGUAGE, temp$NQF)
test
##        
##           0   1
##   ENGL  418 749
##   OTHER  59 124
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 0.73662, df = 1, p-value = 0.3907
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##     Comparison p.Fisher p.adj.Fisher
## 1 ENGL : OTHER    0.361        0.361
plotDat(temp, "NQF", "MARITAL_STATUS", F, "Marital Status", "Marital Status", "Frequency")

test <- table(temp$MARITAL_STATUS, temp$NQF)
test
##          
##             0   1
##   MARRIED 227 334
##   SINGLE   80 190
##   UNKNOWN  23  34
##   WIDOWED 147 315
chisq.test(test)
## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 13.433, df = 3, p-value = 0.003789
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##          Comparison p.Fisher p.adj.Fisher
## 1  MARRIED : SINGLE  0.00273        0.015
## 2 MARRIED : UNKNOWN  1.00000        1.000
## 3 MARRIED : WIDOWED  0.00500        0.015
## 4  SINGLE : UNKNOWN  0.11900        0.238
## 5  SINGLE : WIDOWED  0.56300        0.676
## 6 UNKNOWN : WIDOWED  0.23100        0.347
plotDat(temp, "NQF", "FIRST_CAREUNIT", F, "First Careunit", "First Careunit", "Frequency")

test <- table(temp$FIRST_CAREUNIT, temp$NQF)
test
##        
##           0   1
##   CCU    53 130
##   CSRU   52  22
##   MICU  190 567
##   SICU  113  92
##   TSICU  69  62
chisq.test(test)
## 
##  Pearson's Chi-squared test
## 
## data:  test
## X-squared = 129.85, df = 4, p-value < 2.2e-16
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##      Comparison p.Fisher p.adj.Fisher
## 1    CCU : CSRU 1.70e-09     4.25e-09
## 2    CCU : MICU 3.01e-01     3.34e-01
## 3    CCU : SICU 2.38e-07     4.76e-07
## 4   CCU : TSICU 2.40e-05     4.00e-05
## 5   CSRU : MICU 1.55e-14     7.75e-14
## 6   CSRU : SICU 2.72e-02     3.40e-02
## 7  CSRU : TSICU 1.78e-02     2.54e-02
## 8   MICU : SICU 1.77e-15     1.77e-14
## 9  MICU : TSICU 1.25e-09     4.17e-09
## 10 SICU : TSICU 7.36e-01     7.36e-01
plotDat(temp, "NQF", "ADMISSION_TYPE", F, "Admission Type", "Admission Type", "Frequency")

test <- table(temp$ADMISSION_TYPE, temp$NQF)
test
##            
##               0   1
##   ELECTIVE   54   6
##   EMERGENCY 423 867
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 79.64, df = 1, p-value < 2.2e-16
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##             Comparison p.Fisher p.adj.Fisher
## 1 ELECTIVE : EMERGENCY 2.54e-19     2.54e-19
plotDat(temp, "NQF", "MECHVENT", F, "Mechanical Ventilation", "Mechanical Ventilation", "Frequency")

test <- table(temp$MECHVENT, temp$NQF)
test
##    
##       0   1
##   0 281 589
##   1 196 284
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 9.4909, df = 1, p-value = 0.002065
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##   Comparison p.Fisher p.adj.Fisher
## 1      0 : 1  0.00197      0.00197
plotDat(temp, "NQF", "HOSPITAL_EXPIRE_FLAG", F, "In-Hospital Mortality", "In-Hospital Mortality", "Frequency")

test <- table(temp$HOSPITAL_EXPIRE_FLAG, temp$NQF)
test
##    
##       0   1
##   0 397 678
##   1  80 195
chisq.test(test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  test
## X-squared = 5.5517, df = 1, p-value = 0.01846
pairwiseNominalIndependence(
  as.matrix(test), 
  fisher = T, gtest = F, chisq = F, method = "fdr")
##   Comparison p.Fisher p.adj.Fisher
## 1      0 : 1   0.0162       0.0162
boxplot(temp$AGE ~ temp$NQF, 
        main = "Caremeasure Implementation by Age",
        xlab = "Implementation (1 == Yes)",
        ylab = "Age (Years)")

summary(temp$AGE[temp$NQF == 1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.04   80.35   84.76   84.70   91.40   91.40
sd(temp$AGE[temp$NQF == 1])
## [1] 5.244684
summary(temp$AGE[temp$NQF == 0])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.01   79.30   83.03   83.21   86.63   91.40
sd(temp$AGE[temp$NQF == 0])
## [1] 4.914834
## Parametric
t.test(temp$AGE ~ temp$NQF)
## 
##  Welch Two Sample t-test
## 
## data:  temp$AGE by temp$NQF
## t = -5.1907, df = 1034.1, p-value = 2.521e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.0501541 -0.9253222
## sample estimates:
## mean in group 0 mean in group 1 
##        83.21094        84.69867
## Non-parametric
wilcox.test(temp$AGE ~ temp$NQF)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  temp$AGE by temp$NQF
## W = 174540, p-value = 7.548e-07
## alternative hypothesis: true location shift is not equal to 0
boxplot(temp$SOFA ~ temp$NQF, 
        main = "Caremeasure Implementation by SOFA Score",
        xlab = "Implementation (1 == Yes)",
        ylab = "SOFA Score")

summary(temp$SOFA[temp$NQF == 1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   3.000   4.000   4.881   6.000  18.000
sd(temp$SOFA[temp$NQF == 1])
## [1] 2.962226
summary(temp$SOFA[temp$NQF == 0])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   4.000   4.453   6.000  18.000
sd(temp$SOFA[temp$NQF == 0])
## [1] 2.80313
## Parametric
t.test(temp$SOFA ~ temp$NQF)
## 
##  Welch Two Sample t-test
## 
## data:  temp$SOFA by temp$NQF
## t = -2.6282, df = 1025.7, p-value = 0.008711
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.7476217 -0.1084590
## sample estimates:
## mean in group 0 mean in group 1 
##        4.452830        4.880871
## Non-parametric
wilcox.test(temp$SOFA ~ temp$NQF)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  temp$SOFA by temp$NQF
## W = 189480, p-value = 0.005871
## alternative hypothesis: true location shift is not equal to 0
boxplot(temp$LOS_ICU ~ temp$NQF, 
        main = "Caremeasure Implementation by\nICU Length of Stay (Days)",
        xlab = "Implementation (1 == Yes)",
        ylab = "Length of Stay (Days)")

summary(temp$LOS_ICU[temp$NQF == 1])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0453  1.2880  2.1152  3.6999  3.9584 43.4559
sd(temp$LOS_ICU[temp$NQF == 1])
## [1] 4.795096
summary(temp$LOS_ICU[temp$NQF == 0])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0245  1.2382  2.2827  4.2581  4.9731 60.9390
sd(temp$LOS_ICU[temp$NQF == 0])
## [1] 5.760079
## Parametric
t.test(temp$LOS_ICU ~ temp$NQF)
## 
##  Welch Two Sample t-test
## 
## data:  temp$LOS_ICU by temp$NQF
## t = 1.8026, df = 839.06, p-value = 0.07181
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.04959986  1.16603126
## sample estimates:
## mean in group 0 mean in group 1 
##        4.258062        3.699846
## Non-parametric
wilcox.test(temp$LOS_ICU ~ temp$NQF)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  temp$LOS_ICU by temp$NQF
## W = 219180, p-value = 0.1092
## alternative hypothesis: true location shift is not equal to 0
boxplot(temp$LOS_HOSPITAL ~ temp$NQF, 
        main = "Caremeasure Implementation by\nHOSPITAL Length of Stay (Days)",
        xlab = "Implementation (1 == Yes)",
        ylab = "Length of Stay (Days)")

summary(temp$LOS_HOSPITAL[temp$NQF == 1])
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##   0.4479   3.6472   5.7000   7.8009   9.1396 105.5799
sd(temp$LOS_HOSPITAL[temp$NQF == 1])
## [1] 7.527101
summary(temp$LOS_HOSPITAL[temp$NQF == 0])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7479  4.4063  6.9174  8.9328 11.2104 60.9250
sd(temp$LOS_HOSPITAL[temp$NQF == 0])
## [1] 7.319067
## Parametric
t.test(temp$LOS_HOSPITAL ~ temp$NQF)
## 
##  Welch Two Sample t-test
## 
## data:  temp$LOS_HOSPITAL by temp$NQF
## t = 2.6888, df = 1002.4, p-value = 0.00729
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.305802 1.957909
## sample estimates:
## mean in group 0 mean in group 1 
##        8.932782        7.800926
## Non-parametric
wilcox.test(temp$LOS_HOSPITAL ~ temp$NQF)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  temp$LOS_HOSPITAL by temp$NQF
## W = 238280, p-value = 1.128e-05
## alternative hypothesis: true location shift is not equal to 0

Generalized Linear Model

Table 3 multivariable model (above)

## Change names of reference variables to bring them first in alphabetical order
temp$ETHNICITY <- ifelse(temp$ETHNICITY == "WHITE", "aWHITE", temp$ETHNICITY)
temp$LANGUAGE <- ifelse(temp$LANGUAGE == "ENGL", "aENGL", temp$LANGUAGE)
temp$FIRST_CAREUNIT <- ifelse(temp$FIRST_CAREUNIT == "MICU", "aMICU", temp$FIRST_CAREUNIT)


## Factor Categorical Variables
temp$ETHNICITY <- factor(temp$ETHNICITY)
temp$GENDER <- factor(temp$GENDER)
temp$LANGUAGE <- factor(temp$LANGUAGE)
temp$MARITAL_STATUS <- factor(temp$MARITAL_STATUS)
temp$FIRST_CAREUNIT <- factor(temp$FIRST_CAREUNIT)
temp$ADMISSION_TYPE <- factor(temp$ADMISSION_TYPE)
temp$NQF <- factor(temp$NQF)
temp$MECHVENT <- factor(temp$MECHVENT)

gm_icu <- glm(NQF ~ 
              GENDER +
              AGE +
              ETHNICITY +
              LANGUAGE +
              MARITAL_STATUS +
              SOFA +
              FIRST_CAREUNIT +
              MECHVENT +
              ADMISSION_TYPE, 
                   data = temp,
                   family = binomial(link = "logit"))

(tmp <- model_info(gm_icu))
## Waiting for profiling to be done...
## $`Model Summary`
## 
## Call:
## glm(formula = NQF ~ GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS + 
##     SOFA + FIRST_CAREUNIT + MECHVENT + ADMISSION_TYPE, family = binomial(link = "logit"), 
##     data = temp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3358  -1.0437   0.6387   0.8439   2.5488  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -6.05956    1.15490  -5.247 1.55e-07 ***
## GENDERM                 -0.40635    0.13584  -2.991  0.00278 ** 
## AGE                      0.05495    0.01267   4.336 1.45e-05 ***
## ETHNICITYBLACK          -0.04150    0.23410  -0.177  0.85930    
## ETHNICITYOTHER           0.61838    0.23696   2.610  0.00906 ** 
## LANGUAGEOTHER           -0.28190    0.19263  -1.463  0.14335    
## MARITAL_STATUSSINGLE     0.35225    0.17511   2.012  0.04426 *  
## MARITAL_STATUSUNKNOWN   -0.03530    0.31265  -0.113  0.91012    
## MARITAL_STATUSWIDOWED    0.19182    0.15749   1.218  0.22323    
## SOFA                     0.05650    0.02367   2.387  0.01699 *  
## FIRST_CAREUNITCCU       -0.26935    0.18963  -1.420  0.15550    
## FIRST_CAREUNITCSRU      -1.60540    0.29557  -5.431 5.59e-08 ***
## FIRST_CAREUNITSICU      -1.32409    0.17712  -7.476 7.68e-14 ***
## FIRST_CAREUNITTSICU     -1.18051    0.20748  -5.690 1.27e-08 ***
## MECHVENT1               -0.14468    0.14171  -1.021  0.30728    
## ADMISSION_TYPEEMERGENCY  2.46688    0.45391   5.435 5.49e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1753.6  on 1349  degrees of freedom
## Residual deviance: 1523.5  on 1334  degrees of freedom
## AIC: 1555.5
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## $`OR Summary`
##                                   OR        2.5 %      97.5 %
## (Intercept)              0.002335427 0.0002358047  0.02193287
## GENDERM                  0.666077875 0.5099168894  0.86875235
## AGE                      1.056490137 1.0306957843  1.08322099
## ETHNICITYBLACK           0.959353367 0.6109065420  1.53300326
## ETHNICITYOTHER           1.855926194 1.1754635478  2.98099826
## LANGUAGEOTHER            0.754346922 0.5185676664  1.10476311
## MARITAL_STATUSSINGLE     1.422266518 1.0114206969  2.01054045
## MARITAL_STATUSUNKNOWN    0.965319519 0.5252878163  1.79747967
## MARITAL_STATUSWIDOWED    1.211455236 0.8901372979  1.65101809
## SOFA                     1.058121371 1.0105167797  1.10886513
## FIRST_CAREUNITCCU        0.763879320 0.5289786688  1.11373959
## FIRST_CAREUNITCSRU       0.200809517 0.1110656321  0.35561249
## FIRST_CAREUNITSICU       0.266045981 0.1876220007  0.37592397
## FIRST_CAREUNITTSICU      0.307122856 0.2041447848  0.46101980
## MECHVENT1                0.865301580 0.6559162887  1.14353031
## ADMISSION_TYPEEMERGENCY 11.785596059 5.2020875142 31.76556955
tab_model(gm_icu)
## Warning: package 'bindrcpp' was built under R version 3.4.4
  NQF
Predictors Odds Ratios CI p
(Intercept) 0.00 0.00 – 0.02 <0.001
GENDERM 0.67 0.51 – 0.87 0.003
AGE 1.06 1.03 – 1.08 <0.001
ETHNICITYBLACK 0.96 0.61 – 1.53 0.859
ETHNICITYOTHER 1.86 1.18 – 2.98 0.009
LANGUAGEOTHER 0.75 0.52 – 1.10 0.143
MARITAL_STATUSSINGLE 1.42 1.01 – 2.01 0.044
MARITAL_STATUSUNKNOWN 0.97 0.53 – 1.80 0.910
MARITAL_STATUSWIDOWED 1.21 0.89 – 1.65 0.223
SOFA 1.06 1.01 – 1.11 0.017
FIRST_CAREUNITCCU 0.76 0.53 – 1.11 0.156
FIRST_CAREUNITCSRU 0.20 0.11 – 0.36 <0.001
FIRST_CAREUNITSICU 0.27 0.19 – 0.38 <0.001
FIRST_CAREUNITTSICU 0.31 0.20 – 0.46 <0.001
MECHVENT 0 0.87 0.66 – 1.14 0.307
ADMISSION_TYPEEMERGENCY 11.79 5.20 – 31.77 <0.001
Observations 1350
Cox & Snell’s R2 / Nagelkerke’s R2 0.157 / 0.216 
#write.csv(cbind(tmp$`OR Summary`, tmp$`Model Summary`$coefficients), file = "~/Desktop/temp.csv")

Cox Regression and Survival Curves

## SURVIVAL
## if no death, time is length of data collection, 
## else it is length from first admission to death
## Aggregate NQF by admission, admittime, data of death, and keep expiration info

tmp <- aggregate(cbind(NQF,
                       AGE,
                       SOFA,
                       MECHVENT) ~ 
                     SUBJECT_ID +
                     ETHNICITY +
                     GENDER +
                     MARITAL_STATUS +
                     FIRST_CAREUNIT +
                     HADM_ID +
                     HOSPITAL_EXPIRE_FLAG +
                     ADMISSION_TYPE +
                     EXPIRE_FLAG +
                     LANGUAGE +
                     HADM_ID +
                     DOD +
                     ADMITTIME +
                     EXPIRE_FLAG +
                     HOSPITAL_EXPIRE_FLAG,
                 data = team, 
                 FUN = mean)

## Change names of reference variables to bring them first in alphabetical order
tmp$ETHNICITY <- ifelse(tmp$ETHNICITY == "WHITE", "aWHITE", tmp$ETHNICITY)
tmp$LANGUAGE <- ifelse(tmp$LANGUAGE == "ENGL", "aENGL", tmp$LANGUAGE)
tmp$FIRST_CAREUNIT <- ifelse(tmp$FIRST_CAREUNIT == "MICU", "aMICU", tmp$FIRST_CAREUNIT)


tmp$TIME <- ifelse(is.na(tmp$DOD), NA, (tmp$DOD - tmp$ADMITTIME)/365)

# Data are right censored
nqf_cox <- coxph(Surv(TIME, HOSPITAL_EXPIRE_FLAG == 1, type = "right") ~ 
              NQF +
              GENDER +
              AGE +
              ETHNICITY +
              LANGUAGE +
              MARITAL_STATUS +
              SOFA +
              FIRST_CAREUNIT +
              MECHVENT +
              ADMISSION_TYPE, 
                   data = tmp)

(temp <- summary(nqf_cox))
## Call:
## coxph(formula = Surv(TIME, HOSPITAL_EXPIRE_FLAG == 1, type = "right") ~ 
##     NQF + GENDER + AGE + ETHNICITY + LANGUAGE + MARITAL_STATUS + 
##         SOFA + FIRST_CAREUNIT + MECHVENT + ADMISSION_TYPE, data = tmp)
## 
##   n= 1350, number of events= 275 
## 
##                             coef exp(coef) se(coef)      z Pr(>|z|)    
## NQF                      0.59797   1.81843  0.14469  4.133 3.58e-05 ***
## GENDERM                 -0.09022   0.91373  0.13423 -0.672 0.501511    
## AGE                     -0.03000   0.97045  0.01232 -2.435 0.014899 *  
## ETHNICITYBLACK          -0.31029   0.73323  0.25117 -1.235 0.216679    
## ETHNICITYOTHER          -0.24336   0.78399  0.20788 -1.171 0.241722    
## LANGUAGEOTHER           -0.48293   0.61698  0.20584 -2.346 0.018968 *  
## MARITAL_STATUSSINGLE    -0.17921   0.83593  0.17572 -1.020 0.307805    
## MARITAL_STATUSUNKNOWN    0.85120   2.34246  0.25092  3.392 0.000693 ***
## MARITAL_STATUSWIDOWED   -0.08313   0.92023  0.16015 -0.519 0.603708    
## SOFA                     0.19326   1.21319  0.01844 10.482  < 2e-16 ***
## FIRST_CAREUNITCCU       -0.15436   0.85696  0.20588 -0.750 0.453407    
## FIRST_CAREUNITCSRU       0.70363   2.02108  0.26673  2.638 0.008341 ** 
## FIRST_CAREUNITSICU       0.59848   1.81934  0.17183  3.483 0.000496 ***
## FIRST_CAREUNITTSICU      0.67106   1.95631  0.19652  3.415 0.000639 ***
## MECHVENT                 0.72790   2.07074  0.13821  5.267 1.39e-07 ***
## ADMISSION_TYPEEMERGENCY  1.35229   3.86628  0.47879  2.824 0.004737 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## NQF                        1.8184     0.5499    1.3694    2.4146
## GENDERM                    0.9137     1.0944    0.7024    1.1887
## AGE                        0.9704     1.0305    0.9473    0.9942
## ETHNICITYBLACK             0.7332     1.3638    0.4482    1.1996
## ETHNICITYOTHER             0.7840     1.2755    0.5216    1.1783
## LANGUAGEOTHER              0.6170     1.6208    0.4122    0.9236
## MARITAL_STATUSSINGLE       0.8359     1.1963    0.5924    1.1796
## MARITAL_STATUSUNKNOWN      2.3425     0.4269    1.4325    3.8305
## MARITAL_STATUSWIDOWED      0.9202     1.0867    0.6723    1.2596
## SOFA                       1.2132     0.8243    1.1701    1.2578
## FIRST_CAREUNITCCU          0.8570     1.1669    0.5724    1.2829
## FIRST_CAREUNITCSRU         2.0211     0.4948    1.1982    3.4090
## FIRST_CAREUNITSICU         1.8193     0.5496    1.2991    2.5479
## FIRST_CAREUNITTSICU        1.9563     0.5112    1.3309    2.8755
## MECHVENT                   2.0707     0.4829    1.5794    2.7150
## ADMISSION_TYPEEMERGENCY    3.8663     0.2586    1.5127    9.8818
## 
## Concordance= 0.77  (se = 0.018 )
## Rsquare= 0.173   (max possible= 0.944 )
## Likelihood ratio test= 256.6  on 16 df,   p=<2e-16
## Wald test            = 278.2  on 16 df,   p=<2e-16
## Score (logrank) test = 306  on 16 df,   p=<2e-16
write.csv(cbind(temp$conf.int, temp$coefficients), file = "~/Desktop/temp.csv")

Table 4/or Figure 2 - Kaplan meier plot of 1-year survival stratified by NQF measure (mortality).

We also need to somehow show days in the ICU until death or 1-year (burden), number of readmissions (burden)

All-Cause Mortality (1 Year)

# Data are right censored
fit <- survfit(Surv(TIME, EXPIRE_FLAG == 1, type = "right") ~ NQF,
               data = tmp)

print(fit)
## Call: survfit(formula = Surv(TIME, EXPIRE_FLAG == 1, type = "right") ~ 
##     NQF, data = tmp)
## 
##         n events median 0.95LCL 0.95UCL
## NQF=0 477    477  0.537   0.438   0.707
## NQF=1 873    873  0.238   0.195   0.326
ggsurvplot(fit, data = tmp,
 surv.median.line = "hv", # Add medians survival

 # Change legends: title & labels
 legend.title = "NQF Care Measure",
 legend.labs = c("Care Preferences\n Not Documented", "Care Preferences\n Documented"),
 # Add p-value and intervals
 pval =  "             p < 0.0001",
 conf.int = TRUE,
 
 #Limit (xlim, 1 year)
 xlim = c(0,1),
 xlab = "Time (Years)",
 #break.time.by = 120,
 
 # Add risk table
 risk.table = TRUE,
 tables.height = 0.2,
 tables.theme = theme_cleantable(),

 # Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
 # or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
 palette = cm.colors(2),#c("#E7B800", "#2E9FDF"),
 ggtheme = theme_bw(), break.time.by = 0.25 # Change ggplot2 theme
)

In-Hospital Mortality (1 Year)

#Data are right censored
fit <- survfit(Surv(TIME, HOSPITAL_EXPIRE_FLAG == 1, type = "right") ~ NQF,
               data = tmp)

print(fit)
## Call: survfit(formula = Surv(TIME, HOSPITAL_EXPIRE_FLAG == 1, type = "right") ~ 
##     NQF, data = tmp)
## 
##         n events median 0.95LCL 0.95UCL
## NQF=0 477     80     NA      NA      NA
## NQF=1 873    195     NA      NA      NA
ggsurvplot(fit, data = tmp,
 surv.median.line = "hv", # Add medians survival

 # Change legends: title & labels
 legend.title = "NQF Care Measure",
 legend.labs = c("Care Preferences\n Not Documented", "Care Preferences\n Documented"),
 # Add p-value and intervals
 pval =  "             p < 0.0001",
 conf.int = TRUE,
 
 #Limit (xlim, 1 year)
 xlim = c(0,1),
 xlab = "Time (Years)",
 #break.time.by = 120,
 
 # Add risk table
 risk.table = TRUE,
 tables.height = 0.2,
 tables.theme = theme_cleantable(),

 # Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
 # or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
 palette = cm.colors(2),#c("#E7B800", "#2E9FDF"),
 ggtheme = theme_bw(), break.time.by = 0.25 # Change ggplot2 theme
)
## Warning in .add_surv_median(p, fit, type = surv.median.line, fun = fun, :
## Median survival not reached.