-
Notifications
You must be signed in to change notification settings - Fork 1
/
CTIS_model_fit.Rmd
208 lines (164 loc) · 6.49 KB
/
CTIS_model_fit.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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
---
title: "CTIS_Model"
output: html_document
---
This document contains the models included in our paper.
```{r setup, include=FALSE}
knitr::opts_chunk$set(include = FALSE)
library(tidyverse)
library(party)
library(rpart)
library(caret)
library(rpart.plot)
library(mboost)
library(MLmetrics)
library(pROC)
library(corrplot)
library(olsrr)
library(glmnet)
library(effects)
library(arsenal)
library(mboost)
library(MASS)
library(effects)
library(sjPlot)
library(visreg)
library(glmnet)
library(plotmo)
library(lme4)
library(cAIC4)
library(splines)
```
### Setup for modelling
```{r include=FALSE}
CTIS_micro <- readRDS("data/protected_data/CTIS_microdata_cleanV3.RDS")
CTIS_macro <- readRDS("data/data_CTIS_policy.RDS")
data_CTIS_policy <- readRDS("data/data_CTIS_policy.RDS")
set.seed(2507)
CTIS_micro$anxious <- ifelse(CTIS_micro$D1 %in% levels(CTIS_micro$D1)[4:5], TRUE, FALSE )
CTIS_micro$depressed <- ifelse(CTIS_micro$D2 %in% levels(CTIS_micro$D1)[4:5], TRUE, FALSE )
CTIS_micro$V11[CTIS_micro$E3 %in% "Male" & is.na(CTIS_micro$V11)] <- "No"
CTIS_micro$E5 <- as.numeric(CTIS_micro$E5)
CTIS_micro$E7a <- as.numeric(CTIS_micro$E7a)
```
# Microdata logistic regression model
```{r}
# Set reference categories
CTIS_micro$E4 <- relevel(CTIS_micro$E4, ref = "45-54")
CTIS_micro$E3 <- relevel(CTIS_micro$E3, ref = "Male")
CTIS_micro$E8 <- relevel(CTIS_micro$E8, ref = "Secondary school complete")
CTIS_micro$E2 <- relevel(CTIS_micro$E2, ref = "Town")
# Fit model
# E4 = Age, E3 = Gender, E8 = Education, E2 = Where do you live
model_glm <- glm(anxious ~ E4 + E3 + E8 + E2,family=binomial(link='logit'), data=CTIS_micro)
summary(model_glm)
exp(model_glm$coefficients)
# Forestplot
all.models <- list()
all.models[[1]] <- model_glm
plot_models(all.models)
ggsave("Plots/Chris/models/model_glm.png", width = 12, height = 9)
```
## Elisa
```{r Regularisierung with policy data}
# Model with regularisierung - policy data
# only metric variables
## training (80%) und testdatensatz (20%)
## Lasso model included in paper
data_ctis_policy_model <- data_CTIS_policy %>% drop_na() %>% dplyr::select(anxious_7d, finance, food_security, worried_become_ill, retail_and_recreation, grocery_and_pharmacy, residential, transit_stations, parks, workplaces)
head(data_ctis_policy_model)
n_policy <- nrow(data_ctis_policy_model)
n_policy
p_policy <- ncol(data_ctis_policy_model) - 1
p_policy
set.seed(123)
ind_train_policy <- sample(x = 1:n_policy, size = ceiling(0.8 * n_policy))
set_train_policy <- data_ctis_policy_model[ind_train_policy,] ## dataset
ind_test_policy <- setdiff(x = 1:n_policy, ind_train_policy)
set_test_policy <- data_ctis_policy_model[ind_test_policy, ] ## dataset
## remove nas in set_train
set_train_policy <- set_train_policy %>% drop_na()
## unpenalized model: everything significant
model_unpenalized_policy <- lm(formula = anxious_7d ~ finance + food_security +
worried_become_ill + retail_and_recreation +
grocery_and_pharmacy + residential + transit_stations +
parks + workplaces, data = set_train_policy)
summary(model_unpenalized_policy)
# Ridge-Regularisierung:
model_ridge_policy <- glmnet(x = as.matrix(set_train_policy[, 2:10]), y = set_train_policy$anxious_7d,
alpha = 0)
summary(model_ridge_policy)
# LASSO-Regularisierung:
model_lasso_policy <- glmnet(x = as.matrix(set_train_policy[, 2:10]), y = set_train_policy$anxious_7d,
alpha = 1)
summary(model_lasso_policy)
# lasso coefficients path
plot(model_lasso_policy, label = T)
# Best λ:
lambda_ridge_policy <- cv.glmnet(x = data.matrix(set_train_policy[, 2:10]),
y = set_train_policy$anxious_7d, alpha = 0)$lambda.min
lambda_ridge_policy
lambda_lasso_policy <- cv.glmnet(x = data.matrix(set_train_policy[, 2:10]), y = set_train_policy$anxious_7d, alpha = 1)$lambda.min
lambda_lasso_policy
## plot mse for different lambdas
jpeg(file="Plots/Elisa/mse_lambda_lasso.jpeg",
width = 465, height = 225, units='mm', res = 300)
plot(cv.glmnet(x = data.matrix(set_train_policy[, 2:10]), y = set_train_policy$anxious_7d, alpha = 1), cex.lab = 1.5, cex.axis = 1.5,
cex.main = 1, cex.sub = 1.5)
# ggsave("Plots/Elisa/Model/mse_lambda_lasso.png", width = 7, height = 3)
dev.off()
# Anpassung der Modelle an Trainingsdaten:
y_train_policy <- set_train_policy$anxious_7d
predict_train_policy <- matrix(data = 0, nrow = nrow(set_train_policy), ncol = 3)
predict_train_policy[, 1] <- predict.glmnet(object = model_ridge_policy,
newx = data.matrix(set_train_policy[, 2:10]),
s = lambda_ridge_policy)
predict_train_policy[, 2] <- predict.glmnet(object = model_lasso_policy,
newx = data.matrix(set_train_policy[, 2:10]),
s = lambda_lasso_policy)
MSE_train_policy <- rep(x = 0, length.out = 2)
for (i in 1:2) {
MSE_train_policy[i] = mean((y_train_policy - predict_train_policy[, i])^2)
}
MSE_train_policy
format(MSE_train_policy, scientific = FALSE)
# model Lasso best one (lowest MSE)
## Prädiktion auf Testdaten:
# drop nas in set_test
set_test_policy <- set_test_policy %>% drop_na()
y_test_policy <- set_test_policy$anxious_7d
predict_test_policy <- matrix(data = 0, nrow = nrow(set_test_policy), ncol = 3)
predict_test_policy[, 1] <- predict.glmnet(object = model_ridge_policy,
newx = data.matrix(set_test_policy[, 2:10]),
s = lambda_ridge_policy)
predict_test_policy[, 2] <- predict.glmnet(object = model_lasso_policy,
newx = data.matrix(set_test_policy[, 2:10]), s = lambda_lasso_policy)
MSE_test_policy <- rep(x = 0, length.out = 2)
for (i in 1:2) {
MSE_test_policy[i] = mean((y_test_policy - predict_test_policy[, i])^2)
}
MSE_test_policy
format(MSE_test_policy, scientific = FALSE)
## lasso best one
## Model with lasso regularisierung
coef_lasso_policy <- model_lasso_policy$beta[, which(model_lasso_policy$lambda == lambda_lasso_policy)]
coef_lasso_policy
# Ridge-Regularisierung
# Coefficients path plot
plot_glmnet(x = model_ridge_policy, label = TRUE, xvar = "lambda")
title(main = "Ridge", line = 3)
# LASSO-Regularisierung
# Coefficients path plot
plot_glmnet(x = model_lasso_policy, label = TRUE, xvar = "lambda")
ggsave("Plots/Elisa/Model/lasso-path-coeff.png", width = 7,
height = 3)
coef_frame <- data.frame(coef_lasso_policy = names(coef_lasso_policy),
value = coef_lasso_policy)
# Forest plot: coefficients Lasso model
ggplot(coef_frame, aes(x=coef_lasso_policy, y=value)) +
coord_flip() +
geom_point(colour = "red", size = 3) +
theme(axis.title.y=element_blank(), axis.text =
element_text(size = 12)) + xlab("Coefficients")
ggsave("Plots/Elisa/Model/forest_plot_new.png", width = 7, height = 3)
```