-
Notifications
You must be signed in to change notification settings - Fork 0
/
misccode.R
80 lines (60 loc) · 2.08 KB
/
misccode.R
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
# misc code
# ggplot for lasso
ggplot(tidied_cv, aes(lambda, estimate)) + geom_line(color = "red") +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2) +
scale_x_log10() +
geom_vline(xintercept = glance_cv$lambda.min) +
geom_vline(xintercept = glance_cv$lambda.1se, lty = 2)
## Backward Logistic
```{r hccBackwardLog, eval=F}
hccbackward <- list()
for (i in 1:length(cohorts)) {
hccmodeldata <- train %>% filter(cohort==cohorts[[i]])
glmhccfull = glm(hccfm, data=hccmodeldata, family="binomial")
hccbackward[[i]] <- step(glmhccfull, direction = "backward")
}
```
## Forward Logistic
```{r hccForwardLog, eval=F}
library(foreach)
library(doMC)
registerDoMC(cores=8)
hccforward <- vector('list', 1)
for(i in seq_along(hccforward)){
a <- replicate(4, list())
hccforward[[i]] <- a
}
hccforward <- list(c(0))
for (i in 1:2) {
hccmodeldata <- train %>% filter(cohort==cohorts[[i]]) %>% ungroup()
glmhccfull = glm(hccfm, data=hccmodeldata, family="binomial")
glmhccempty = glm(isreadmit30dc ~ 1, data=hccmodeldata, family="binomial")
for (j in 1:2) {
hccforward[[i]] <- list(hccforward[[i]], step(glmhccempty, scope = list(
lower = formula(glmhccempty),
upper = formula(glmhccfull)),
direction = "forward",
steps=j))
}
}
t <- step(glmhccempty, scope = list(
lower = formula(glmhccempty),
upper = formula(glmhccfull)),
direction = "forward",
steps=1)[c("aic", "call")]
t <- hccmodeldata %>% ungroup %>% select(isreadmit30dc, starts_with("hcc_"))
glmhccfull = glm(hccfm, data=hccmodeldata, family="binomial")
glmhccempty = glm(isreadmit30dc ~ 1, data=hccmodeldata, family="binomial")
for (i in 1:10) {
hccforward[[i]] <- step(glmhccempty, scope = list(
lower = formula(glmhccempty),
upper = formula(glmhccfull)),
direction = "forward",
steps = i)
}
bestglm(Xy = t, family=binomial, IC = "AIC", method="forward")
?bestglm
# Generate plot of AIC vs. Number of predictors for each cohort
hccmodeldata <- train %>% filter(cohort==cohorts[[1]])
hccf <- regsubsets(hccfm, data=hccmodeldata, method="backward", nvmax = 49)
```