Validation of ICD-9-CM/ICD-10-CM Codes for Automated Electronic Scoring of APACHE II, APACHE III, and SAPS II

Analysis for ASHP Poster

Patient Selection

• Patients meeting inclusion criteria: 16242
• Reasons for exclusion:
  • Pregnant: 203
  • Prisoner: 77
  • ICU stay < 12 hours: 1077
  • Diagnosis coding contained both ICD-9-CM and ICD-10-CM: 548
  • Missing data needed to calculate risk scores: 9165
  • Re-encounters (only included each patient once) : 1782
• Total number of eligible patients remaining to sample from: 3390

Patient Demographics

deceased <- "Deceased/Hospice"
home <- "Home"
transfer <- "Transferred"
dispos <- list("Deceased" = deceased,
               "Expired/Donor" = deceased, 
               "Hospice-Home" = deceased, 
               "Hospice-Medical Facility" = deceased,
               "Home Care with Home Health" = home,
               "Home or Self Care" = home,
               "Left Against Medical Advise" = home,
               "Skilled Nursing Facility" = transfer,
               "DC/DISC TO REHAB" = transfer,
               "Long Term Care" = transfer,
               "Intermediate Care" = transfer,
               "Acute Care" = transfer,
               "Federal Healthcare Facility" = transfer,
               "DC/TF To Psych Hosp" = transfer)
df <- data_demographics %>%
    dmap_at("disposition", str_replace_all, pattern = dispos) %>%
    dmap_at("race", factor) %>%
    dmap_at("race", fct_lump, n = 4) %>%
    select(-pie.id, Age = age, Sex = sex, Race = race, Disposition = disposition, `Hospital Length of Stay (days)` = length.stay, `ICU Location` = location) 
vars <- names(df)
tbl <- CreateTableOne(vars, data = df)
ptbl <- print(tbl, printToggle = FALSE)
rownames(ptbl) <- str_replace_all(rownames(ptbl), "   ", "- ")
knitr::kable(ptbl, caption = "Demographics")

	Overall
n	200
Age (mean (sd))	60.68 (16.54)
Sex = Male (%)	123 (61.5)
Race (%)
- African American	49 (28.0)
- Asian	8 ( 4.6)
- White/Caucasian	73 (41.7)
- Other	45 (25.7)
Disposition (%)
- Deceased/Hospice	37 (18.5)
- Home	109 (54.5)
- Transferred	54 (27.0)
Hospital Length of Stay (days) (mean (sd))	12.07 (12.04)
ICU Location (%)
- Cullen 2 E Medical Intensive Care Unit	60 (30.0)
- Hermann 3 Shock Trauma Intensive Care Unit	14 ( 7.0)
- Hermann 3 Transplant Surgical ICU	19 ( 9.5)
- HVI Cardiac Care Unit	20 (10.0)
- HVI Cardiovascular Intensive Care Unit	34 (17.0)
- Jones 7 J Elective Neuro ICU	53 (26.5)

Primary endpoint

APACHE II

• Mean and median score for each comorbidity set

scores <- ls(pattern = "score_apache2")
all_scores <- map(scores, ~mutate(get(.x), score = .x)) %>%
    map(~select(.x, pie.id, score, apache2)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "score_apache2_", replacement = "") 
man <- filter(all_scores, score == "manual")$apache2
result <- all_scores %>%
    select(-pie.id) %>%
    group_by(score) %>%
    summarize_all(funs(mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(result)

score	mean	sd	median	IQR_low	IQR_high
ahrq	20.20	6.858220	20	15	24.00
ahrq_drg	19.79	6.832225	19	15	24.00
ek	20.05	6.909865	20	15	24.00
ek_drg	19.65	6.866824	19	15	24.00
elixhauser	20.18	6.927304	20	15	24.00
elixhauser_drg	19.77	6.900375	19	15	24.00
manual	19.87	7.020651	19	15	24.25
quan	20.18	6.927304	20	15	24.00
quan_drg	19.78	6.891959	19	15	24.00

ggplot(all_scores, aes(x = score, y = apache2)) +
    geom_boxplot() 

    # geom_jitter(color = "light blue", alpha = 0.6)

ggplot(all_scores, aes(x = apache2)) +
    geom_histogram(binwidth = 0.5) +
    facet_wrap(~ score)

df <- all_scores %>%
    spread(score, apache2) %>%
    gather(score, apache2, -pie.id, -manual)
# ggplot(df, aes(x = apache2, y = manual)) +
#     geom_jitter(aes(color = score), alpha = 0.6) 
ggplot(df, aes(x = apache2, y = manual)) +
    geom_point() +
    # geom_smooth() +
    facet_wrap(~ score)

• Comparison of scores
  • Results run using both the t-test and Wilcoxon signed rank test (assumes non-parametric data, is more difficult to achieve significance)
  • To account for performing multiple comparisons, a BH (Benjamini and Hochberg) correction is applied to the p-values; use the bh.adjust value as the p-value
    • Reference: http://ies.ed.gov/ncee/pubs/20084018/app_b.asp

t <- all_scores[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(t.test(.$apache2, man, paired = TRUE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(t, digits = 4)

score	estimate	statistic	p.value	bh.adjust	parameter	conf.low	conf.high	method	alternative
ahrq	0.33	2.0591	0.0408	0.1306	199	0.0140	0.6460	Paired t-test	two.sided
ahrq_drg	-0.08	-0.4981	0.6190	0.6190	199	-0.3967	0.2367	Paired t-test	two.sided
ek	0.18	1.1170	0.2653	0.4245	199	-0.1378	0.4978	Paired t-test	two.sided
ek_drg	-0.22	-1.3754	0.1706	0.3411	199	-0.5354	0.0954	Paired t-test	two.sided
elixhauser	0.31	1.9808	0.0490	0.1306	199	0.0014	0.6186	Paired t-test	two.sided
elixhauser_drg	-0.10	-0.6386	0.5238	0.6190	199	-0.4088	0.2088	Paired t-test	two.sided
quan	0.31	1.9808	0.0490	0.1306	199	0.0014	0.6186	Paired t-test	two.sided
quan_drg	-0.09	-0.5734	0.5670	0.6190	199	-0.3995	0.2195	Paired t-test	two.sided

mw <- all_scores[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(wilcox.test(.$apache2, man, paired = TRUE, conf.int = TRUE, exact = FALSE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(mw, digits = 4)

score	estimate	statistic	p.value	bh.adjust	conf.low	conf.high	method	alternative
ahrq	1.5	1644.0	0.0563	0.1676	-0.0001	2.0000	Wilcoxon signed rank test with continuity correction	two.sided
ahrq_drg	0.0	1070.0	0.5179	0.5179	-1.5000	0.0000	Wilcoxon signed rank test with continuity correction	two.sided
ek	0.0	1448.5	0.3147	0.5036	-0.0001	2.0000	Wilcoxon signed rank test with continuity correction	two.sided
ek_drg	0.0	936.0	0.1357	0.2715	-1.5000	0.0000	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser	1.5	1551.0	0.0628	0.1676	-0.0001	1.9999	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser_drg	0.0	990.0	0.4483	0.5179	-1.5000	0.0001	Wilcoxon signed rank test with continuity correction	two.sided
quan	1.5	1551.0	0.0628	0.1676	-0.0001	1.9999	Wilcoxon signed rank test with continuity correction	two.sided
quan_drg	0.0	1030.5	0.4863	0.5179	-1.5000	0.0000	Wilcoxon signed rank test with continuity correction	two.sided

options(scipen = 999)
mw_pair <- all_scores %>%
    dmap_at("score", factor) %>%
    dmap_at("score", ~fct_relevel(.x, c("manual", "ek", "ek_drg"))) %>%
    do(tidy(pairwise.wilcox.test(.$apache2, .$score, p.adj = "BH", paired = TRUE, exact = FALSE)))

Linear Regression Model

lm <- df %>%
    group_by(score) %>%
    do(fit = lm(manual ~ apache2, data = .))
knitr::kable(glance(lm, fit), digits = 4)

score	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
ahrq	0.8967	0.8962	2.2624	1718.366	0	2	-446.0658	898.1315	908.0265	1013.4322	198
ahrq_drg	0.8960	0.8955	2.2695	1706.410	0	2	-446.6916	899.3832	909.2782	1019.7944	198
ek	0.8960	0.8955	2.2693	1706.705	0	2	-446.6761	899.3522	909.2472	1019.6365	198
ek_drg	0.8971	0.8966	2.2574	1726.889	0	2	-445.6220	897.2440	907.1390	1008.9449	198
elixhauser	0.9020	0.9015	2.2035	1822.143	0	2	-440.7920	887.5840	897.4790	961.3710	198
elixhauser_drg	0.9016	0.9011	2.2078	1814.341	0	2	-441.1790	888.3579	898.2529	965.0983	198
quan	0.9020	0.9015	2.2035	1822.143	0	2	-440.7920	887.5840	897.4790	961.3710	198
quan_drg	0.9011	0.9006	2.2135	1804.009	0	2	-441.6937	889.3874	899.2824	970.0789	198

augment(lm, fit) %>%
    ggplot(aes(x = .fitted, y = .resid)) + 
    geom_point() +
    # geom_smooth(se = FALSE) +
    facet_wrap(~ score)

mod_df <- df %>%
    filter(score == "ek_drg")
mod <- lm(manual ~ apache2, data = mod_df)
par(mfrow = c(2, 2))
plot(mod)

• Compare differences in score from comorbidity sets vs. manual
  • Column num_diff indicates the number of patients with a score that differed from manual
  • Mean and median indicate the difference from the manual score
    • These were only calculated for those patients whose score was different from manual
    • If all patients were included, the number is skewed towards 0 by all the patients with the same scores

diff <- all_scores %>%
    spread(score, apache2) %>%
    mutate_at(vars(-pie.id, -manual), funs(abs(. - manual))) %>%
    gather(score, apache2, -pie.id, -manual) %>%
    select(-pie.id, -manual) %>%
    filter(apache2 > 0) %>%
    group_by(score) %>%
    summarize_all(funs(num_diff = length, mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(diff, caption = "Differences in APACHE II Scores")

score	num_diff	mean	sd	median	IQR_low	IQR_high
ahrq	72	3.500000	1.510526	3.5	2	5
ahrq_drg	68	3.588235	1.508536	5.0	2	5
ek	71	3.521127	1.510526	5.0	2	5
ek_drg	68	3.588235	1.508536	5.0	2	5
elixhauser	70	3.457143	1.510214	2.0	2	5
elixhauser_drg	66	3.545454	1.510800	5.0	2	5
quan	70	3.457143	1.510214	2.0	2	5
quan_drg	67	3.522388	1.511153	5.0	2	5

APACHE III

scores <- ls(pattern = "score_apache3")
all_scores <- map(scores, ~mutate(get(.x), score = .x)) %>%
    map(~select(.x, pie.id, score, apache3)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "score_apache3_", replacement = "") 
man <- filter(all_scores, score == "manual")$apache3
result <- all_scores %>%
    select(-pie.id) %>%
    group_by(score) %>%
    summarize_all(funs(mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(result)

score	mean	sd	median	IQR_low	IQR_high
ahrq	85.715	29.29976	82.0	65.75	103.25
ahrq_drg	85.555	29.25932	82.0	65.00	103.25
ek	85.715	29.29976	82.0	65.75	103.25
ek_drg	85.555	29.25932	82.0	65.00	103.25
elixhauser	85.715	29.29976	82.0	65.75	103.25
elixhauser_drg	85.555	29.25932	82.0	65.00	103.25
manual	84.805	28.99847	81.5	63.75	103.00
quan	85.715	29.29976	82.0	65.75	103.25
quan_drg	85.555	29.25932	82.0	65.00	103.25

ggplot(all_scores, aes(x = score, y = apache3)) +
    geom_boxplot()

df <- all_scores %>%
    spread(score, apache3) %>%
    gather(score, apache3, -pie.id, -manual)
ggplot(df, aes(x = apache3, y = manual)) +
    geom_point() +
    # geom_smooth() +
    facet_wrap(~ score)

• Comparison of scores

t <- all_scores[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(t.test(.$apache3, man, paired = TRUE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(t, digits = 4)

score	estimate	statistic	p.value	bh.adjust	parameter	conf.low	conf.high	method	alternative
ahrq	0.91	2.9503	0.0036	0.0071	199	0.3018	1.5182	Paired t-test	two.sided
ahrq_drg	0.75	2.5453	0.0117	0.0117	199	0.1690	1.3310	Paired t-test	two.sided
ek	0.91	2.9503	0.0036	0.0071	199	0.3018	1.5182	Paired t-test	two.sided
ek_drg	0.75	2.5453	0.0117	0.0117	199	0.1690	1.3310	Paired t-test	two.sided
elixhauser	0.91	2.9503	0.0036	0.0071	199	0.3018	1.5182	Paired t-test	two.sided
elixhauser_drg	0.75	2.5453	0.0117	0.0117	199	0.1690	1.3310	Paired t-test	two.sided
quan	0.91	2.9503	0.0036	0.0071	199	0.3018	1.5182	Paired t-test	two.sided
quan_drg	0.75	2.5453	0.0117	0.0117	199	0.1690	1.3310	Paired t-test	two.sided

mw <- all_scores[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(wilcox.test(.$apache3, man, paired = TRUE, conf.int = TRUE, exact = FALSE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(mw, digits = 4)

score	estimate	statistic	p.value	bh.adjust	conf.low	conf.high	method	alternative
ahrq	9.0000	283.5	0.0062	0.0123	2.9999	13.0	Wilcoxon signed rank test with continuity correction	two.sided
ahrq_drg	7.0001	251.0	0.0176	0.0176	1.5000	11.5	Wilcoxon signed rank test with continuity correction	two.sided
ek	9.0000	283.5	0.0062	0.0123	2.9999	13.0	Wilcoxon signed rank test with continuity correction	two.sided
ek_drg	7.0001	251.0	0.0176	0.0176	1.5000	11.5	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser	9.0000	283.5	0.0062	0.0123	2.9999	13.0	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser_drg	7.0001	251.0	0.0176	0.0176	1.5000	11.5	Wilcoxon signed rank test with continuity correction	two.sided
quan	9.0000	283.5	0.0062	0.0123	2.9999	13.0	Wilcoxon signed rank test with continuity correction	two.sided
quan_drg	7.0001	251.0	0.0176	0.0176	1.5000	11.5	Wilcoxon signed rank test with continuity correction	two.sided

Linear Regression Model

lm <- df %>%
    group_by(score) %>%
    do(fit = lm(manual ~ apache3, data = .))
knitr::kable(glance(lm, fit), digits = 4)

score	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
ahrq	0.9778	0.9777	4.3280	8735.812	0	2	-575.8020	1157.604	1167.499	3708.786	198
ahrq_drg	0.9797	0.9796	4.1402	9564.716	0	2	-566.9292	1139.858	1149.753	3393.891	198
ek	0.9778	0.9777	4.3280	8735.812	0	2	-575.8020	1157.604	1167.499	3708.786	198
ek_drg	0.9797	0.9796	4.1402	9564.716	0	2	-566.9292	1139.858	1149.753	3393.891	198
elixhauser	0.9778	0.9777	4.3280	8735.812	0	2	-575.8020	1157.604	1167.499	3708.786	198
elixhauser_drg	0.9797	0.9796	4.1402	9564.716	0	2	-566.9292	1139.858	1149.753	3393.891	198
quan	0.9778	0.9777	4.3280	8735.812	0	2	-575.8020	1157.604	1167.499	3708.786	198
quan_drg	0.9797	0.9796	4.1402	9564.716	0	2	-566.9292	1139.858	1149.753	3393.891	198

augment(lm, fit) %>%
    ggplot(aes(x = .fitted, y = .resid)) + 
    geom_point() +
    facet_wrap(~ score)

• Differences in score

diff <- all_scores %>%
    spread(score, apache3) %>%
    mutate_at(vars(-pie.id, -manual), funs(abs(. - manual))) %>%
    gather(score, apache3, -pie.id, -manual) %>%
    select(-pie.id, -manual) %>%
    filter(apache3 > 0) %>%
    group_by(score) %>%
    summarize_all(funs(num_diff = length, mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(diff, caption = "Difference in APACHE III Scores")

score	num_diff	mean	sd	median	IQR_low	IQR_high
ahrq	26	11.38462	4.825573	12	10	16
ahrq_drg	25	10.88000	5.035871	12	6	16
ek	26	11.38462	4.825573	12	10	16
ek_drg	25	10.88000	5.035871	12	6	16
elixhauser	26	11.38462	4.825573	12	10	16
elixhauser_drg	25	10.88000	5.035871	12	6	16
quan	26	11.38462	4.825573	12	10	16
quan_drg	25	10.88000	5.035871	12	6	16

SAPS II

scores <- ls(pattern = "score_saps2")
all_scores <- map(scores, ~mutate(get(.x), score = .x)) %>%
    map(~select(.x, pie.id, score, saps2)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "score_saps2_", replacement = "") 
man <- filter(all_scores, score == "manual")$saps2
result <- all_scores %>%
    select(-pie.id) %>%
    group_by(score) %>%
    summarize_all(funs(mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(result)

score	mean	sd	median	IQR_low	IQR_high
ahrq	37.12	14.56386	36.0	26	47.00
ahrq_drg	37.12	14.56386	36.0	26	47.00
ek	37.12	14.56386	36.0	26	47.00
ek_drg	37.12	14.56386	36.0	26	47.00
elixhauser	37.12	14.56386	36.0	26	47.00
elixhauser_drg	37.12	14.56386	36.0	26	47.00
manual	36.94	14.44888	34.5	26	46.25
quan	37.12	14.56386	36.0	26	47.00
quan_drg	37.12	14.56386	36.0	26	47.00

ggplot(all_scores, aes(x = score, y = saps2)) +
    geom_boxplot()

df <- all_scores %>%
    spread(score, saps2) %>%
    gather(score, saps2, -pie.id, -manual)
ggplot(df, aes(x = saps2, y = manual)) +
    geom_point() +
    # geom_smooth() +
    facet_wrap(~ score)

• Comparison of scores

t <- all_scores[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(t.test(.$saps2, man, paired = TRUE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(t, digits = 4)

score	estimate	statistic	p.value	bh.adjust	parameter	conf.low	conf.high	method	alternative
ahrq	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided
ahrq_drg	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided
ek	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided
ek_drg	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided
elixhauser	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided
elixhauser_drg	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided
quan	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided
quan_drg	0.18	1.3406	0.1816	0.1816	199	-0.0848	0.4448	Paired t-test	two.sided

• Confidence interval was unable to be calculated for Wilcoxon analysis here

mw <- all_scores[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(wilcox.test(.$saps2, man, paired = TRUE, exact = FALSE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(mw, digits = 4)

score	statistic	p.value	bh.adjust	method	alternative
ahrq	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided
ahrq_drg	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided
ek	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided
ek_drg	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser_drg	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided
quan	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided
quan_drg	23	0.5182	0.5182	Wilcoxon signed rank test with continuity correction	two.sided

Linear Regression Model

lm <- df %>%
    group_by(score) %>%
    do(fit = lm(manual ~ saps2, data = .))
knitr::kable(glance(lm, fit), digits = 4)

score	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual
ahrq	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198
ahrq_drg	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198
ek	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198
ek_drg	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198
elixhauser	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198
elixhauser_drg	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198
quan	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198
quan_drg	0.983	0.9829	1.8886	11449.91	0	2	-409.9484	825.8968	835.7917	706.2182	198

augment(lm, fit) %>%
    ggplot(aes(x = .fitted, y = .resid)) + 
    geom_point() +
    facet_wrap(~ score)

• Differences in score

diff <- all_scores %>%
    spread(score, saps2) %>%
    mutate_at(vars(-pie.id, -manual), funs(abs(. - manual))) %>%
    gather(score, saps2, -pie.id, -manual) %>%
    select(-pie.id, -manual) %>%
    filter(saps2 > 0) %>%
    group_by(score) %>%
    summarize_all(funs(num_diff = length, mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(diff, caption = "Difference in SAPS II Scores")

score	num_diff	mean	sd	median	IQR_low	IQR_high
ahrq	8	9.5	0.5345225	9.5	9	10
ahrq_drg	8	9.5	0.5345225	9.5	9	10
ek	8	9.5	0.5345225	9.5	9	10
ek_drg	8	9.5	0.5345225	9.5	9	10
elixhauser	8	9.5	0.5345225	9.5	9	10
elixhauser_drg	8	9.5	0.5345225	9.5	9	10
quan	8	9.5	0.5345225	9.5	9	10
quan_drg	8	9.5	0.5345225	9.5	9	10

Secondary endpoints

Predicted mortality

SAPS II

• Mean and median risk of mortality (as a percent) for each comorbidity set

comorbid <- ls(pattern = "score_saps2_")
all_df <- map(comorbid, ~mutate(get(.x), score = .x)) %>%
    map(~select(.x, pie.id, score, saps2)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "score_saps2_", replacement = "") %>%
    distinct(pie.id, score, .keep_all = TRUE) %>%
    mutate(logit = -7.763 + 0.0737 * saps2 + 0.9971 * log(saps2 + 1),
           mortality = exp(logit) / (1 + exp(logit))) %>%
    select(pie.id, score, mortality) 
man <- filter(all_df, score == "manual")$mortality
result <- all_df %>%
    select(-pie.id) %>%
    group_by(score) %>%
    summarize_all(funs(mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(result, digits = 4)

score	mean	sd	median	IQR_low	IQR_high
ahrq	0.2513	0.2299	0.1810	0.0717	0.3919
ahrq_drg	0.2513	0.2299	0.1810	0.0717	0.3919
ek	0.2513	0.2299	0.1810	0.0717	0.3919
ek_drg	0.2513	0.2299	0.1810	0.0717	0.3919
elixhauser	0.2513	0.2299	0.1810	0.0717	0.3919
elixhauser_drg	0.2513	0.2299	0.1810	0.0717	0.3919
manual	0.2476	0.2280	0.1597	0.0717	0.3752
quan	0.2513	0.2299	0.1810	0.0717	0.3919
quan_drg	0.2513	0.2299	0.1810	0.0717	0.3919

t <- all_df[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>% 
    do(tidy(t.test(.$mortality, man, paired = TRUE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(t, digits = 4)

score	estimate	statistic	p.value	bh.adjust	parameter	conf.low	conf.high	method	alternative
ahrq	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided
ahrq_drg	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided
ek	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided
ek_drg	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided
elixhauser	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided
elixhauser_drg	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided
quan	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided
quan_drg	0.0037	1.8305	0.0687	0.0687	199	-0.0003	0.0077	Paired t-test	two.sided

mw <- all_df[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(wilcox.test(.$mortality, man, paired = TRUE, conf.int = TRUE, exact = FALSE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(mw, digits = 4)

score	estimate	statistic	p.value	bh.adjust	conf.low	conf.high	method	alternative
ahrq	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided
ahrq_drg	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided
ek	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided
ek_drg	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser_drg	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided
quan	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided
quan_drg	0.1146	31	0.0801	0.0801	-0.0088	0.1865	Wilcoxon signed rank test with continuity correction	two.sided

• Difference in risk of mortality
  • Mean and median calculated using only those patients whose predicted mortality score was different from manual

diff <- all_df %>%
    spread(score, mortality) %>%
    mutate_at(vars(-pie.id, -manual), funs(abs(. - manual))) %>%
    gather(score, saps2, -pie.id, -manual) %>%
    select(-pie.id, -manual) %>%
    filter(saps2 > 0) %>%
    group_by(score) %>%
    summarize_all(funs(num_diff = length, mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(diff, caption = "Difference in Predicted Mortality using SAPS II")

score	num_diff	mean	sd	median	IQR_low	IQR_high
ahrq	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449
ahrq_drg	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449
ek	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449
ek_drg	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449
elixhauser	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449
elixhauser_drg	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449
quan	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449
quan_drg	8	0.1353385	0.0544702	0.1366591	0.0862261	0.1691449

Number of comorbidities

APACHE II

• Mean represents the average number of comorbidities per patient

comorbid <- ls(pattern = "data_apache2_")
all_df <- map(comorbid, ~mutate(get(.x), score = .x)) %>%
    map(~distinct(.x, pie.id, score, liver, cardiovasc, respiratory, renal, immunocomp, comorbidity)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "data_apache2_", replacement = "") %>%
    mutate_if(is.logical, as.numeric) %>%
    by_row(function(x) sum(x[, 3:7], na.rm = TRUE), .collate = "rows", .to = "num_comorbid") %>%
    select(pie.id, score, num_comorbid) 
man <- filter(all_df, score == "manual")$num_comorbid
result <- all_df %>%
    select(-pie.id) %>%
    group_by(score) %>%
    summarize_all(funs(mean, sd, median, IQR_low = quantile(., 0.25), IQR_high = quantile(., 0.75)))
knitr::kable(result)

score	mean	sd	median	IQR_low	IQR_high
ahrq	0.760	0.7314885	1	0	1
ahrq_drg	0.595	0.6732319	0	0	1
ek	0.690	0.7186671	1	0	1
ek_drg	0.530	0.6412143	0	0	1
elixhauser	0.745	0.7366739	1	0	1
elixhauser_drg	0.575	0.6607252	0	0	1
manual	0.605	0.7222070	0	0	1
quan	0.750	0.7417892	1	0	1
quan_drg	0.585	0.6670832	0	0	1

t <- all_df[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(t.test(.$num_comorbid, man, paired = TRUE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(t, digits = 4)

score	estimate	statistic	p.value	bh.adjust	parameter	conf.low	conf.high	method	alternative
ahrq	0.155	2.5932	0.0102	0.0513	199	0.0371	0.2729	Paired t-test	two.sided
ahrq_drg	-0.010	-0.1737	0.8623	0.8623	199	-0.1236	0.1036	Paired t-test	two.sided
ek	0.085	1.4253	0.1556	0.2997	199	-0.0326	0.2026	Paired t-test	two.sided
ek_drg	-0.075	-1.3232	0.1873	0.2997	199	-0.1868	0.0368	Paired t-test	two.sided
elixhauser	0.140	2.3598	0.0193	0.0513	199	0.0230	0.2570	Paired t-test	two.sided
elixhauser_drg	-0.030	-0.5336	0.5942	0.7923	199	-0.1409	0.0809	Paired t-test	two.sided
quan	0.145	2.4379	0.0157	0.0513	199	0.0277	0.2623	Paired t-test	two.sided
quan_drg	-0.020	-0.3528	0.7246	0.8281	199	-0.1318	0.0918	Paired t-test	two.sided

mw <- all_df[, -1] %>%
    filter(score != "manual") %>%
    group_by(score) %>%
    do(tidy(wilcox.test(.$num_comorbid, man, paired = TRUE, conf.int = TRUE, exact = FALSE))) %>%
    ungroup() %>%
    mutate(bh.adjust = p.adjust(p.value, "BH")) %>%
    select(score:p.value, bh.adjust, everything())
knitr::kable(mw, digits = 4)

score	estimate	statistic	p.value	bh.adjust	conf.low	conf.high	method	alternative
ahrq	0	3059.5	0.0077	0.0521	0.0001	1.0000	Wilcoxon signed rank test with continuity correction	two.sided
ahrq_drg	0	1933.0	0.9118	0.9118	0.0000	0.0000	Wilcoxon signed rank test with continuity correction	two.sided
ek	0	2553.5	0.1276	0.2552	0.0000	0.5000	Wilcoxon signed rank test with continuity correction	two.sided
ek_drg	0	1560.5	0.2053	0.3285	-0.5000	0.0000	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser	0	2752.5	0.0196	0.0521	0.0000	0.5000	Wilcoxon signed rank test with continuity correction	two.sided
elixhauser_drg	0	1672.0	0.5878	0.7837	0.0000	0.0001	Wilcoxon signed rank test with continuity correction	two.sided
quan	0	2826.0	0.0161	0.0521	0.0000	0.5000	Wilcoxon signed rank test with continuity correction	two.sided
quan_drg	0	1791.0	0.7129	0.8148	0.0000	0.0000	Wilcoxon signed rank test with continuity correction	two.sided

Proportion with differing comorbidity

APACHE II

• Proportion of patients where a comorbidity was identified by ICD or by manual, but not both

comorbid <- ls(pattern = "data_apache2_")
all_df <- map(comorbid, ~mutate(get(.x), score = .x)) %>%
    map(~distinct(.x, pie.id, score, comorbidity)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "data_apache2_", replacement = "") %>%
    spread(score, comorbidity) %>%
    mutate_at(vars(-pie.id, -manual), funs(. != manual)) %>%
    mutate_if(is.logical, as.numeric) %>%
    ungroup() %>%
    gather(score, comorbidity, -pie.id, -manual) %>%
    select(-pie.id, -manual) %>%
    group_by(score) %>%
    summarize_all(funs(num_diff = sum(. > 0), proportion = mean))
knitr::kable(all_df)

score	num_diff	proportion
ahrq	72	0.360
ahrq_drg	68	0.340
ek	71	0.355
ek_drg	68	0.340
elixhauser	70	0.350
elixhauser_drg	66	0.330
quan	70	0.350
quan_drg	67	0.335

APACHE III

• Proportion of patients where the primary comorbidity identified by ICD varied from primary comorbidity identified by manual

comorbid <- ls(pattern = "data_apache3_")
all_df <- map(comorbid, ~mutate(get(.x), score = .x)) %>%
    map(~distinct(.x, pie.id, score, comorbidity)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "data_apache3_", replacement = "") %>%
    spread(score, comorbidity, fill = "") %>%
    dmap(str_replace_all, pattern = "none", replacement = "") %>%
    mutate_at(vars(-pie.id, -manual), funs(. != manual)) %>%
    mutate_if(is.logical, as.numeric) %>%
    gather(score, comorbidity, -pie.id, -manual) %>%
    select(-pie.id, -manual) %>%
    group_by(score) %>%
    summarize_all(funs(num_diff = sum(. > 0), proportion = mean))
knitr::kable(all_df)

score	num_diff	proportion
ahrq	32	0.160
ahrq_drg	31	0.155
ek	32	0.160
ek_drg	31	0.155
elixhauser	32	0.160
elixhauser_drg	31	0.155
quan	32	0.160
quan_drg	31	0.155

SAPS II

• Proportion of patients where the primary comorbidity identified by ICD varied from primary comorbidity identified by manual

comorbid <- ls(pattern = "data_saps2_")
all_df <- map(comorbid, ~mutate(get(.x), score = .x)) %>%
    map(~distinct(.x, pie.id, score, comorbidity)) %>%
    map_df(~semi_join(.x, manual_patients, by = "pie.id")) %>%
    arrange(pie.id, score) %>%
    dmap_at("score", str_replace_all, pattern = "data_saps2_", replacement = "") %>%
    spread(score, comorbidity, fill = "") %>%
    dmap(str_replace_all, pattern = "none", replacement = "") %>%
    mutate_at(vars(-pie.id, -manual), funs(. != manual)) %>%
    mutate_if(is.logical, as.numeric) %>%
    gather(score, comorbidity, -pie.id, -manual) %>%
    select(-pie.id, -manual) %>%
    group_by(score) %>%
    summarize_all(funs(num_diff = sum(. > 0), proportion = mean))
knitr::kable(all_df)

score	num_diff	proportion
ahrq	8	0.04
ahrq_drg	8	0.04
ek	8	0.04
ek_drg	8	0.04
elixhauser	8	0.04
elixhauser_drg	8	0.04
quan	8	0.04
quan_drg	8	0.04