-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathutil.R
154 lines (145 loc) · 5.7 KB
/
util.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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
library(matrixStats)
ages <- read_csv("data/combined_long_data.csv") %>%
filter(AgeGroup == "children") %>%
pull(Age)
AGE_MEAN <- mean(ages)
AGE_SD <- sd(ages)
recode_for_models <- function(d) {
if(length(unique(d$Age)) > 1) {
d$Age <- (d$Age - AGE_MEAN) / AGE_SD
}
d %>%
mutate(
urb_rural_thoughts=case_when(type == "feelings" ~ 0,
Urb_rural == "urban" ~ -0.5,
Urb_rural == "rural" ~ +0.5,
Urb_rural == 0 ~ 0),
urb_rural_feelings=case_when(type == "thoughts" ~ 0,
Urb_rural == "urban" ~ -0.5,
Urb_rural == "rural" ~ +0.5,
Urb_rural == 0 ~ 0),
gender_thoughts=case_when(type == "feelings" ~ 0,
Gender == "diverse" ~ 0,
Gender == "male" ~ -0.5,
Gender == "female" ~ +0.5),
gender_feelings=case_when(type == "thoughts" ~ 0,
Gender == "diverse" ~ 0,
Gender == "male" ~ -0.5,
Gender == "female" ~ +0.5),
type=if_else(type == "thoughts", -0.5, 0.5),
Gender=case_when(Gender == "diverse" ~ 0,
Gender == "male" ~ -0.5,
Gender == "female" ~ +0.5)
) %>%
return()
}
get_preds <- function(model, nd, re_formula=NULL) {
preds <- posterior_linpred(model, newdata=nd, re_formula=re_formula)
draws <- as_draws_df(model)
ints1 <- draws$`b_Intercept[1]`
ints2 <- draws$`b_Intercept[2]`
positive_response <- pnorm(preds, ints2)
negative_response <- 1 - pnorm(preds, ints1)
partial_response <- pnorm(preds, ints1) - pnorm(preds, ints2)
stopifnot(all(abs((positive_response + negative_response + partial_response) - 1) < 1e5))
nd$pos_mean <- colMeans(positive_response)
nd$pos_lower <- colQuantiles(positive_response, probs=0.025)
nd$pos_upper <- colQuantiles(positive_response, probs=0.975)
nd$neg_mean <- colMeans(negative_response)
nd$neg_lower <- colQuantiles(negative_response, probs=0.025)
nd$neg_upper <- colQuantiles(negative_response, probs=0.975)
nd$part_mean <- colMeans(partial_response)
nd$part_lower <- colQuantiles(partial_response, probs=0.025)
nd$part_upper <- colQuantiles(partial_response, probs=0.975)
nd$pos_most <- colMeans((positive_response > partial_response) & (positive_response > negative_response))
nd$neg_most <- colMeans((negative_response > partial_response) & (negative_response > positive_response))
return(nd)
}
get_country_mean_preds <- function(data, model, type, ages, split_urb_rural=FALSE) {
result_t <- tibble(Country_full = character(), Age = numeric(),
pos_mean = numeric(), pos_lower = numeric(), pos_upper = numeric(),
neg_mean = numeric(), neg_lower = numeric(), neg_upper = numeric(),
part_mean = numeric(), part_lower = numeric(), part_upper = numeric()
)
if(split_urb_rural) {
result_t$Urb_rural <- " "
countcomm_key <- data %>%
filter(type == type) %>%
select(Country_full, Community_ID, Urb_rural) %>%
unique()
} else {
countcomm_key <- data %>%
filter(type == type) %>%
select(Country_full, Community_ID) %>%
unique()
}
countries <- unique(countcomm_key$Country_full)
for(country in countries) {
if(is.na(country)) { next }
if(split_urb_rural) {
for(urb_rural in c("urban", "rural")) {
communities <- countcomm_key %>%
filter(Country_full == country, Urb_rural == urb_rural) %>%
pull(Community_ID) %>%
unique
if(length(communities) == 0) { next }
nd <- expand_grid(Community_ID=communities, type=type,
Age=ages, gender_thoughts=0, gender_feelings=0)
nd <- nd %>%
left_join(countcomm_key, by="Community_ID")
print(c(country, urb_rural))
print(nd)
print(recode_for_models(nd))
preds <- posterior_epred(model, newdata=recode_for_models(nd))
for(age in ages) {
x <- which(nd$Age == age)
#if(length(x) == 1) {
meaning_func <- mean
result_t <- add_row(result_t, Country_full = country, Urb_rural=urb_rural, Age = age,
pos_mean = mean(preds[,x,3]),
pos_lower = quantile(preds[,x,3], p=0.025),
pos_upper = quantile(preds[,x,3], p=0.975),
neg_mean = mean(preds[,x,1]),
neg_lower = quantile(preds[,x,1], p=0.025),
neg_upper = quantile(preds[,x,1], p=0.975),
part_mean = mean(preds[,x,2]),
part_lower = quantile(preds[,x,2], p=0.025),
part_upper = quantile(preds[,x,2], p=0.975)
)
}
print(result_t)
}
} else {
communities <- countcomm_key %>%
filter(Country_full == country) %>%
pull(Community_ID) %>%
unique
if(length(communities) == 0) { next }
nd <- expand_grid(Community_ID=communities, type=type, Urb_rural=0,
Age=ages, gender_thoughts=0, gender_feelings=0)
nd <- nd %>%
left_join(countcomm_key, by="Community_ID")
preds <- posterior_epred(model, newdata=recode_for_models(nd))
for(age in ages) {
x <- which(nd$Age == age)
#if(length(x) == 1) {
meaning_func <- mean
result_t <- add_row(result_t, Country_full = country, Age = age,
pos_mean = mean(preds[,x,3]),
pos_lower = quantile(preds[,x,3], p=0.025),
pos_upper = quantile(preds[,x,3], p=0.975),
neg_mean = mean(preds[,x,1]),
neg_lower = quantile(preds[,x,1], p=0.025),
neg_upper = quantile(preds[,x,1], p=0.975),
part_mean = mean(preds[,x,2]),
part_lower = quantile(preds[,x,2], p=0.025),
part_upper = quantile(preds[,x,2], p=0.975)
)
}
}
}
result_t <- result_t %>%
mutate(result_t, Country_full = if_else(Country_full == "Galapagos_I", "Ecuador", Country_full)) %>%
rename(Country = Country_full)
return(result_t)
}