-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path04b_model_testing.R
167 lines (149 loc) · 9.03 KB
/
04b_model_testing.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
155
156
157
158
159
160
161
162
163
164
165
166
167
library(dplyr)
library(tidyr)
library(brms)
library(rstan)
sim_func <- function(nsquare = 50, nplotpersquare = 5){
nplot <- nsquare*nplotpersquare
sim_data <- data.frame(SQUARE = gl(nplot/nplotpersquare,nplotpersquare)) %>%
mutate(REP_ID = paste0(SQUARE,"X",1:nplot),
Improved = rbinom(nplot,1,0.007*as.numeric(SQUARE))) %>%
mutate(pH = rnorm(nplot,mean = 4.5 + Improved + rep(rnorm(nplot/nplotpersquare,0,0.5),each = nplotpersquare), 1)) %>%
mutate(Ell = rnorm(nplot, mean = pH + rep(rnorm(nplot/nplotpersquare,0,1),each = nplotpersquare), 1)) %>%
mutate(Sdep_year12 = rep(rnorm(nplot/nplotpersquare,0,1),each = nplotpersquare)) %>%
mutate(Ndep_year12 = rep(rnorm(nplot/nplotpersquare,-Sdep_year12,1),each = nplotpersquare)) %>%
mutate(N_year12 = rnorm(nplot,
0.5*Ndep_year12 + rep(rnorm(nplot/nplotpersquare,0,0.2),each = nplotpersquare), 0.5)) %>%
mutate(rain_year12 = rep(rnorm(nplot/nplotpersquare,0,1),each = nplotpersquare),
PH_SE_year12 = 0.2,
ELL_SE_year12 = 0.15 + Improved*0.05,
PH_SE_year23 = 0.15,
ELL_SE_year23 = 0.2 + Improved*0.05,
PH_SE_year34 = 0.15,
ELL_SE_year34 = 0.15 + Improved*0.05) %>%
mutate(pH_diffYear12 = rstudent_t(nplot, 4,
0.5*Improved+(1-Improved)*0.5*Sdep_year12 +
0.3*rain_year12,PH_SE_year12))%>%
mutate(Ell_diffYear12 = rstudent_t(nplot, 4,
pH_diffYear12*(1-Improved) +
Improved*0.5*pH_diffYear12 +
(pH>5.5)*-0.5*N_year12,ELL_SE_year12)) %>%
mutate(pH_year2 = pH + pH_diffYear12,
Ell_year2 = Ell + Ell_diffYear12,
rain_year23 = rep(rnorm(nplot/nplotpersquare,0,1),each = nplotpersquare),
Sdep_year23 = rep(rnorm(nsquare, Sdep_year12[seq(1,nplot,nplotpersquare)], 0.1),each=5),
Ndep_year23 = rep(rnorm(nsquare, Ndep_year12[seq(1,nplot,nplotpersquare)], 0.1), each = 5)) %>%
mutate(N_year23 = N_year12 + rnorm(nplot, -0.1*N_year12 +
0.5*Ndep_year23,
0.1)) %>%
mutate(pH_diffYear23 = rstudent_t(nplot, 4,
0.5*Improved+(1-Improved)*0.5*Sdep_year23 +
0.3*rain_year23, PH_SE_year23)) %>%
mutate(Ell_diffYear23 = rstudent_t(nplot, 4,
pH_diffYear23*(1-Improved) +
Improved*0.5*pH_diffYear23 +
(pH_year2>5.5)*-0.5*N_year23,ELL_SE_year23)) %>%
mutate(pH_year3 = pH_year2 + pH_diffYear23,
Ell_year3 = Ell_year2 + Ell_diffYear23,
rain_year34 = rep(rnorm(nplot/nplotpersquare,-0.5,1),each = nplotpersquare),
Sdep_year34 = rep(rnorm(nsquare, Sdep_year23[seq(1,nplot,nplotpersquare)], 0.1),each=5),
Ndep_year34 = rep(rnorm(nsquare, Ndep_year23[seq(1,nplot,nplotpersquare)], 0.1), each = 5)) %>%
mutate(N_year34 = N_year23 + rnorm(nplot, -0.1*N_year23 +
0.5*Ndep_year34,
0.1)) %>%
mutate(pH_diffYear34 = rstudent_t(nplot, 4,
0.5*Improved+(1-Improved)*0.5*Sdep_year34 +
0.3*rain_year34, PH_SE_year34)) %>%
mutate(Ell_diffYear34 = rstudent_t(nplot, 4,
pH_diffYear34*(1-Improved) +
Improved*0.5*pH_diffYear34 +
(pH_year3>5.5)*-0.5*N_year34,ELL_SE_year34)) %>%
select(SQUARE, REP_ID, Improved,
Sdep_year12, Sdep_year23, Sdep_year34, Ndep_year12, Ndep_year23, Ndep_year34,
rain_year12, rain_year23, rain_year34,
N_year12, N_year23, N_year34,
PH_SE_year12, PH_SE_year23, PH_SE_year34,
ELL_SE_year12, ELL_SE_year23, ELL_SE_year34,
PH_year12 = pH_diffYear12, PH_year23 = pH_diffYear23, PH_year34 = pH_diffYear34,
Year1_pH_year12 = pH, Year1_pH_year23 = pH_year2, Year1_pH_year34 = pH_year3,
Ell_year12 = Ell_diffYear12, Ell_year23 = Ell_diffYear23, Ell_year34 = Ell_diffYear34) %>%
pivot_longer(contains("_year"), names_to = c("Variable","Time_period"),
names_sep = "_year") %>%
pivot_wider(names_from = Variable, values_from = value) %>%
mutate(YRnm = as.numeric(as.factor(Time_period)))
mod_pr <- c(prior(normal(0,0.5), class = "b", resp = "Ell"),
prior(normal(0,0.5), class = "b", resp = "PH"),
prior(normal(0,0.5), class = "b", resp = "N"),
prior(student_t(3, 0, 2.5), class = "sds", resp = "Ell"),
prior(normal(0,0.25), class = "Intercept", resp = "Ell"),
prior(normal(0,0.25), class = "Intercept", resp = "PH"),
prior(student_t(3,0,1), class = "Intercept", resp = "N"),
prior(gamma(4,1), class = "nu", resp = "Ell"),
prior(gamma(4,1), class = "nu", resp = "PH"),
prior(normal(0,0.1), class = "ar", resp = "Ell"),
prior(normal(0,0.1), class = "ar", resp = "PH"),
prior(normal(0,0.1), class = "ar", resp = "N"),
prior(student_t(3,0,0.5), class = "sd", resp = "Ell"),
prior(student_t(3,0,0.5), class = "sd", resp = "PH"),
prior(student_t(3,0,0.5), class = "sd", resp = "N"),
prior(student_t(3,0,0.5), class = "sigma", resp = "Ell"),
prior(student_t(3,0,0.5), class = "sigma", resp = "PH"),
prior(student_t(3,0,0.5), class = "sigma", resp = "N"))
sim_mod <- brm(bf(Ell | mi(ELL_SE) ~ Improved*mi(PH) + s(Year1_pH, N, Improved) +
(1|SQUARE) +
ar(time = YRnm, gr = REP_ID),
family = "student") +
bf(PH | mi(PH_SE) ~ Improved*Sdep + rain + (1|SQUARE) +
ar(time = YRnm, gr = REP_ID), family = "student") +
bf(N ~ Ndep + (1|SQUARE) +
ar(time = YRnm, gr = REP_ID)) +
set_rescor(FALSE), data = sim_data, prior = mod_pr,
save_pars = save_pars(all = TRUE, latent = TRUE),
cores = 4, iter = 4000)
sampler_params <- rstan::get_sampler_params(sim_mod$fit, inc_warmup = FALSE)
divergent <- c(sampler_params[[1]][,'divergent__'],
sampler_params[[2]][,'divergent__'],
sampler_params[[3]][,'divergent__'],
sampler_params[[4]][,'divergent__'])
# rerun if there are divergent transitions
# happens once, if setting adapt delta to 0.99 doesn't fix it then discard model results later
if(sum(divergent)>0 & max(rhat(sim_mod))<1.05){
sim_mod <- update(sim_mod, cores = 4, iter = 4000,
save_pars = save_pars(all = TRUE, latent = TRUE),
control = list(adapt_delta = 0.99))
divergent <- rstan::get_sampler_params(sim_mod$fit, inc_warmup=FALSE)[[1]][,'divergent__']
}
# rerun with more iterations if neff too low - only do this if Rhat acceptable
if(min(neff_ratio(sim_mod))<0.05 & min(neff_ratio(sim_mod))>0.02 & max(rhat(sim_mod))<1.05){
sim_mod <- update(sim_mod, cores = 4, iter = 2000,
save_pars = save_pars(all = TRUE, latent = TRUE))
}
# if divergent okay, rhat okay and neff okay then calculate if estimates
# within the 50 and 90% CI
# Parameters:
# [1] "Ell_Intercept" "PH_Intercept" "N_Intercept"
# [4] "Ell_Improved" "PH_Improved" "PH_Sdep"
# [7] "PH_rain" "PH_Improved:Sdep" "N_Ndep"
# [10] "Ell_sYear1_pHNImproved_1" "Ell_sYear1_pHNImproved_2" "Ell_sYear1_pHNImproved_3"
# [13] "Ell_sYear1_pHNImproved_4" "Ell_sYear1_pHNImproved_5" "Ell_sYear1_pHNImproved_6"
# [16] "Ell_sYear1_pHNImproved_7" "Ell_sYear1_pHNImproved_8" "Ell_sYear1_pHNImproved_9"
# [19] "Ell_miPH" "Ell_miPH:Improved"
sim_pars <- c(0,0,0,0.5,0.5,0.5,0.3,-0.5,0.5,0,0,0,0,0,0,0,0,0,1,-0.5)
if(sum(divergent) == 0 & max(rhat(sim_mod))<1.05 & min(neff_ratio(sim_mod))>0.05){
par_vals <- fixef(sim_mod, probs = c(0.05,0.25,0.75,0.95))
CI90 <- sim_pars > par_vals[,"Q5"] & sim_pars < par_vals[,"Q95"]
CI50 <- sim_pars > par_vals[,"Q25"] & sim_pars < par_vals[,"Q75"]
} else{
CI90 <- rep(NA, 21)
CI50 <- rep(NA, 21)
}
return(cbind(CI50,CI90))
}
set.seed(7)
brm_model_test <- replicate(10, sim_func(), simplify = FALSE)
test <- make_conditions(sim_data, vars = c("Improved","Year1_pH"))[1:6,]
test[,"Improved"] <- c(0,0,0,1,1,1)
test[,"cond__"] <- gsub("-0.2","0",test[,"cond__"])
test[,"cond__"] <- gsub("0.2","1",test[,"cond__"])
plot(conditional_effects(sim_mod, conditions = test, effects = "N",
select_points = 0.5),
points = TRUE, plot = FALSE)[[1]]