-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path2_fit_models.R
executable file
·139 lines (113 loc) · 7.69 KB
/
2_fit_models.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
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
#!/usr/bin/env Rscript
library(tidyverse)
library(brms)
source("util.R")
set.seed(6789)
### CHILD MODELS
# Load data, turn responses into ordered factors, separate out child responses
d_child <- read_csv("data/combined_long_data.csv") %>%
filter(AgeGroup == "children") %>%
mutate(Age=scale(Age), response_has = factor(response_has, levels = c("No", "Partially", "Yes"), labels=c("No", "Partially", "Yes"), exclude=NA, ordered = TRUE),
response_similar = factor(response_similar, levels = c("Different", "Partially", "Same"), labels=c("Different", "Partially", "Same"), exclude=NA, ordered = TRUE))
# Model formulae
fixed_only_formula <- "outcome ~ type*Age + urb_rural_thoughts*Age + urb_rural_feelings*Age + gender_thoughts*Age + gender_feelings*Age"
mixed_formula <- paste(fixed_only_formula, "+ (1 + type*Age + urb_rural_thoughts*Age + urb_rural_feelings*Age + gender_thoughts*Age + gender_feelings*Age | Country_full) + (1 + type*Age + gender_thoughts + gender_feelings | Community_ID)")
null_formula <- "outcome ~ Age + Gender + (1 + Age + Gender | Country_full + Community_ID)"
# Model priors
int_prior <- prior(normal(0,5.0), class="b")
age_prior <- prior(normal(0,1.25), class="b", coef="Age")
age_int_prior<- prior(normal(0,0.625), class="b", coef="type:Age") +
prior(normal(0,0.625), class="b", coef="Age:gender_feelings") +
prior(normal(0,0.625), class="b", coef="Age:gender_thoughts") +
prior(normal(0,0.625), class="b", coef="Age:urb_rural_feelings") +
prior(normal(0,0.625), class="b", coef="Age:urb_rural_thoughts")
full_fixef_prior <- int_prior + age_prior + age_int_prior
ranef_prior <- prior(exponential(5), class="sd", group="Country_full") +
prior(exponential(10), class="sd", group="Community_ID")
full_mixed_prior <- full_fixef_prior + ranef_prior
child_null_prior <- int_prior + age_prior + ranef_prior
# Fit models
## Children has
child_has_full <- brm(as.formula(str_replace(mixed_formula, "outcome", "response_has")),
data=recode_for_models(d_child), family=cumulative(),
prior=full_mixed_prior, cores=4) %>% add_criterion("loo")
saveRDS(child_has_full, "m_children_has.rds")
child_has_fixed <- brm(as.formula(str_replace(fixed_only_formula, "outcome", "response_has")),
data=recode_for_models(d_child), family=cumulative(),
prior=full_fixef_prior, cores=4) %>% add_criterion("loo")
saveRDS(child_has_fixed, "m_children_has_fixed_only.rds")
child_has_null <- brm(as.formula(str_replace(null_formula, "outcome", "response_has")),
data=recode_for_models(d_child), family=cumulative(),
prior=child_null_prior, cores=4) %>% add_criterion("loo")
saveRDS(child_has_null, "m_children_has_null.rds")
## Children sim
child_sim_full <- brm(as.formula(str_replace(mixed_formula, "outcome", "response_similar")),
data=recode_for_models(d_child), family=cumulative(),
prior=full_mixed_prior, cores=4) %>% add_criterion("loo")
saveRDS(child_sim_full, "m_children_sim.rds")
child_sim_fixed <- brm(as.formula(str_replace(fixed_only_formula, "outcome", "response_similar")),
data=recode_for_models(d_child), family=cumulative(),
prior=full_fixef_prior, cores=4) %>% add_criterion("loo")
saveRDS(child_sim_fixed, "m_children_sim_fixed_only.rds")
child_sim_null <- brm(as.formula(str_replace(null_formula, "outcome", "response_similar")),
data=recode_for_models(d_child), family=cumulative(),
prior=child_null_prior, cores=4) %>% add_criterion("loo")
saveRDS(child_sim_null, "m_children_sim_null.rds")
### ADULT MODELS
# Load data, turn responses into ordered factors, separate out child responses
d_adult <- read_csv("data/combined_long_data.csv") %>%
filter(AgeGroup == "adults") %>%
mutate(response_has = factor(response_has, levels = c("No", "Partially", "Yes"), labels=c("No", "Partially", "Yes"), exclude=NA, ordered = TRUE),
response_similar = factor(response_similar, levels = c("Different", "Partially", "Same"), labels=c("Different", "Partially", "Same"), exclude=NA, ordered = TRUE))
# Model formulae
fixed_only_formula <- "outcome ~ type*Age + urb_rural_thoughts*Age + urb_rural_feelings*Age + gender_thoughts*Age + gender_feelings*Age"
mixed_formula <- paste(fixed_only_formula, "+ (1 + type*Age + urb_rural_thoughts*Age + urb_rural_feelings*Age + gender_thoughts*Age + gender_feelings*Age | Country_full) + (1 + type*Age + gender_thoughts + gender_feelings | Community_ID)")
null_formula <- "outcome ~ type + urb_rural_thoughts + urb_rural_feelings + gender_thoughts + gender_feelings + (1 + type + urb_rural_thoughts + urb_rural_feelings + gender_thoughts + gender_feelings | Country_full) + (1 + type + gender_thoughts + gender_feelings | Community_ID)"
super_null_formula <- "outcome ~ type + urb_rural_thoughts + urb_rural_feelings + gender_thoughts + gender_feelings"
# Model priors
int_prior <- prior(normal(0,5.0), class="b")
age_prior <- prior(normal(0,1.25), class="b", coef="Age")
age_int_prior<- prior(normal(0,0.625), class="b", coef="type:Age") +
prior(normal(0,0.625), class="b", coef="Age:gender_feelings") +
prior(normal(0,0.625), class="b", coef="Age:gender_thoughts") +
prior(normal(0,0.625), class="b", coef="Age:urb_rural_feelings") +
prior(normal(0,0.625), class="b", coef="Age:urb_rural_thoughts")
full_fixef_prior <- int_prior + age_prior + age_int_prior
ranef_prior <- prior(exponential(5), class="sd", group="Country_full") +
prior(exponential(10), class="sd", group="Community_ID")
full_mixed_prior <- full_fixef_prior + ranef_prior
adult_null_prior <- int_prior + ranef_prior
# Fit models
adult_has_full <- brm(as.formula(str_replace(mixed_formula, "outcome", "response_has")),
data=recode_for_models(d_adult), family=cumulative(),
prior=full_mixed_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_has_full, "m_adults_has.rds")
adult_has_fixed <- brm(as.formula(str_replace(fixed_only_formula, "outcome", "response_has")),
data=recode_for_models(d_adult), family=cumulative(),
prior=full_fixef_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_has_fixed, "m_adults_has_fixed_only.rds")
adult_has_null <- brm(as.formula(str_replace(null_formula, "outcome", "response_has")),
data=recode_for_models(d_adult), family=cumulative(),
prior=adult_null_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_has_null, "m_adults_has_null.rds")
adult_has_super_null <- brm(as.formula(str_replace(super_null_formula, "outcome", "response_has")),
data=recode_for_models(d_adult), family=cumulative(),
prior=int_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_has_super_null, "m_adults_has_super_null.rds")
## Adult sim
adult_sim_full <- brm(as.formula(str_replace(mixed_formula, "outcome", "response_similar")),
data=recode_for_models(d_adult), family=cumulative(),
prior=full_mixed_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_sim_full, "m_adults_sim.rds")
adult_sim_fixed <- brm(as.formula(str_replace(fixed_only_formula, "outcome", "response_similar")),
data=recode_for_models(d_adult), family=cumulative(),
prior=full_fixef_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_sim_fixed, "m_adults_sim_fixed_only.rds")
adult_sim_null <- brm(as.formula(str_replace(null_formula, "outcome", "response_similar")),
data=recode_for_models(d_adult), family=cumulative(),
prior=adult_null_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_sim_null, "m_adults_sim_null.rds")
adult_sim_super_null <- brm(as.formula(str_replace(super_null_formula, "outcome", "response_similar")),
data=recode_for_models(d_adult), family=cumulative(),
prior=int_prior, cores=4) %>% add_criterion("loo")
saveRDS(adult_sim_null, "m_adults_sim_null.rds")