-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmodels.Rmd
121 lines (93 loc) · 3.36 KB
/
models.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
---
title: "Analysis"
output: html_notebook
---
```{r}
library(ggplot2)
library(dplyr)
library(purrr)
library(broom)
library(pROC)
library(data.table)
library(plotROC)
df <- read.csv("outpout_data_firsticu.csv") %>% as.data.table
setDT(df)
```
Define pesi
```{r}
df[, pesi:=(age>80) + cancer + lung_disease + (HeartRate_Mean >= 110) + (SysBP_Mean < 100) + (SpO2_Mean < 90)]
df[, auc(roc(DODin30~pesi))]
```
```{r}
pesi_vars <- c('age', 'cancer', 'lung_disease', 'HeartRate_Mean', 'SysBP_Mean', 'SpO2_Mean')
# map(pesi_vars, ~sum(is.na(df[[.x]])))
df <- df[df[, complete.cases(.SD), .SDcols=c(pesi_vars, "DODin30")],]
```
Marginals
```{r}
lr1 <- glm(DODin30~factor(pesi)+Afib, data = df, family = 'binomial')
summary(lr1)
lr2 <- glm(DODin30~age+cancer+lung_disease+(HeartRate_Mean >= 110)+(SysBP_Mean<100)+(SpO2_Mean),
data = df, family = 'binomial')
lr3 <- glm(DODin30~age+cancer+lung_disease+(HeartRate_Mean >= 110)+(SysBP_Mean<100)+(SpO2_Mean) + Afib,
data = df, family = 'binomial')
lr4 <- glm(DODin30~age+cancer + Afib,
data = df, family = 'binomial')
# anova(lr1, lr2, lr3, test="Chisq")
anova(lr2, lr3, test="Chisq")
anova(lr3, lr4, test="Chisq")
df[, auc(roc(DODin30 ~ predict(lr1)))]
df[, auc(roc(DODin30 ~ predict(lr2)))]
df[, auc(roc(DODin30 ~ predict(lr3)))]
d1 <- drop1(lr3)
var_df <- data.frame(variable=rownames(d1), AIC = d1$AIC)
var_df$AIC_diff <- (var_df$AIC - var_df[1, ]$AIC)
var_df <- arrange(var_df, AIC_diff)
library(magrittr)
var_df %>%
filter(variable!='<none>') %>%
ggplot(aes(variable, AIC_diff, fill=AIC_diff>0)) + geom_col() +
theme_minimal() + ggtitle("variable importance") + theme(legend.position = 'none')
attr(d1, 'heading')
df[, table(Afib, DODin30)]
df[, prop.table(table(Afib, DODin30), margin=1)]
df[, summary(table(Afib, DODin30))]
```
```{r}
df[, lp1:=predict(lr1)]
df[, pesi_refitted:=predict(lr2)]
df[, pesi_with_af:=predict(lr3)]
dfm <- melt.data.table(df[, .SD, .SDcols=c("DODin30", "pesi", "pesi_refitted", "pesi_with_af")], id.vars="DODin30", variable.name = 'model', value.name = 'prediction')
dfm %>%
filter(model=='pesi') %>%
ggplot(aes(m=prediction, d=DODin30, col=model)) + geom_roc(labels=F) +
theme_minimal() + theme(legend.position='bottom') +
ggtitle("AUC - pesi")
dfm %>%
filter(model %in% c('pesi', 'pesi_refitted')) %>%
ggplot(aes(m=prediction, d=DODin30, col=model)) + geom_roc(labels=F) +
theme_minimal() + theme(legend.position='bottom') +
ggtitle("AUC - pesi")
dfm %>%
filter(model %in% c('pesi', 'pesi_refitted')) %>%
ggplot(aes(m=prediction, d=DODin30, col=model)) + geom_roc(labels=F) +
theme_minimal() + theme(legend.position='bottom') +
ggtitle("AUC - pesi (refitted)")
dfm %>%
filter(model %in% c('pesi', 'pesi_refitted')) %>%
ggplot(aes(m=prediction, d=DODin30, col=model)) + geom_roc(labels=F) +
theme_minimal() + theme(legend.position='bottom') +
ggtitle("AUC - pesi (refitted)")
dfm %>%
ggplot(aes(m=prediction, d=DODin30, col=model)) + geom_roc(labels=F) +
theme_minimal() + theme(legend.position='bottom')
```
```{r}
library(rms)
fit2 <- lrm(DODin30~age+cancer+lung_disease+(HeartRate_Mean >= 110)+(SysBP_Mean<100)+(SpO2_Mean), data = df, x=T, y=T)
cal2 <- calibrate(fit2)
plot(cal2)
fit3 <- lrm(DODin30~age+cancer+lung_disease+(HeartRate_Mean >= 110)+(SysBP_Mean<100)+(SpO2_Mean) + Afib, data = df, x=T, y=T)
cal3 <- calibrate(fit3)
plot(cal3)
```