-
Notifications
You must be signed in to change notification settings - Fork 0
/
figures.R
348 lines (276 loc) · 18.1 KB
/
figures.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
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
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
#Fig. 1: time series####
#DM data
dmcont_dat=rodent_data%>%
mutate(abundance=rowSums(.[6:8]))%>%
select(newmoonnumber,treatment, abundance)%>%
replace_na(list(treatment='control'))%>%
filter(treatment=="control")
dmexcl_dat=rodent_data%>%
mutate(abundance=rowSums(.[6:8]))%>%
select(newmoonnumber,treatment, abundance)%>%
replace_na(list(treatment='exclosure'))%>%
filter(treatment=="exclosure")
dmcont_dat1=dmcont_dat%>%filter(!newmoonnumber<278, !newmoonnumber>526)%>%mutate(species="Dipodomys spp.")
dmexcl_dat1=dmexcl_dat%>%filter(!newmoonnumber<278, !newmoonnumber>526)%>%mutate(species="Dipodomys spp.")
dmcont_dat1$abundance=round_na.interp(dmcont_dat1$abundance)
dmexcl_dat1$abundance=round_na.interp(dmexcl_dat1$abundance)
dmdat=rbind(dmcont_dat1, dmexcl_dat1)
#PB data
pbcont_dat=rodent_data%>%
select(newmoonnumber,treatment, PB)%>%rename("abundance"="PB")%>%
replace_na(list(treatment='control'))%>%
filter(treatment=="control")
pbexcl_dat=rodent_data%>%
select(newmoonnumber,treatment, PB)%>%rename("abundance"="PB")%>%
replace_na(list(treatment='exclosure'))%>%
filter(treatment=="exclosure")
pbcont_dat1=pbcont_dat%>%mutate(species="C. baileyi")%>%
filter(!newmoonnumber<278, !newmoonnumber>526)
pbexcl_dat1=pbexcl_dat%>%mutate(species="C. baileyi")%>%
filter(!newmoonnumber<278, !newmoonnumber>526)
pbcont_dat1$abundance=round_na.interp(pbcont_dat1$abundance)
pbexcl_dat1$abundance=round_na.interp(pbexcl_dat1$abundance)
pbdat=rbind(pbcont_dat1, pbexcl_dat1)
#PP data
ppcont_dat=rodent_data%>%
select(newmoonnumber,treatment, PP)%>%rename("abundance"="PP")%>%
replace_na(list(treatment='control'))%>%
filter(treatment=="control")
ppexcl_dat=rodent_data%>%
select(newmoonnumber,treatment, PP)%>%rename("abundance"="PP")%>%
replace_na(list(treatment='exclosure'))%>%
filter(treatment=="exclosure")
ppcont_dat1=ppcont_dat%>%mutate(species="C. penicillatus")%>%
filter(!newmoonnumber<278, !newmoonnumber>526)
ppexcl_dat1=ppexcl_dat%>%mutate(species="C. penicillatus")%>%
filter(!newmoonnumber<278, !newmoonnumber>526)
ppcont_dat1$abundance=round_na.interp(ppcont_dat1$abundance)
ppexcl_dat1$abundance=round_na.interp(ppexcl_dat1$abundance)
ppdat=rbind(ppcont_dat1, ppexcl_dat1)
spdat=rbind(dmdat, pbdat, ppdat)
spdat1=spdat%>%mutate(trt=case_when(treatment=="control"~"control",
treatment=="exclosure"~"removal"))
#plot
pp11=ggplot(spdat1)+ geom_line(mapping=aes(x=newmoonnumber, y=abundance, col=species))+
theme(legend.text = element_text(face="italic"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+
facet_wrap(~trt, nrow=2)+ scale_color_viridis(discrete=TRUE)+
xlab("Sampling number (new moon number)")+ ylab("Abundance")+
geom_bracket(
xmin = 280 , xmax = 396, y.position = 60,
label = paste("italic(C.baileyi)"), type="expression"
)+
geom_bracket(
xmin = 411 , xmax = 526, y.position = 60,
label = paste("italic(C.penicillatus)"), type="expression"
)+
labs(color="Species")
#PP results####
##Fig. 2: coefficient overlap####
int_ov=ggdensity(coef_df_PP, x="intercept", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, ylab="Probability density",
add="median",xlab=F, size=1)+geom_vline(xintercept=0)+
scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio=1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
ar1_ov=ggdensity(coef_df_PP, x="beta1", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab=F, size=1, ylab=F
)+geom_vline(xintercept=0)+scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio=1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
ar12_ov=ggdensity(coef_df_PP, x="beta12", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab=F, size=1, ylab=F
)+geom_vline(xintercept=0)+scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio=1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
temp_ov=ggdensity(coef_df_PP, x="temp", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab=F, size=1, ylab=F
)+geom_vline(xintercept=0)+scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio=1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
wprec_ov=ggdensity(coef_df_PP, x="warm_precip", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab=F, size=1, ylab=F
)+geom_vline(xintercept=0)+scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
scale_y_continuous(n.breaks=3)+
theme(aspect.ratio=1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
cprec_ov=ggdensity(coef_df_PP, x="cool_precip", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab=F, size=1, ylab=F
)+geom_vline(xintercept=0)+scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio=1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
coef_ov=ggarrange(int_ov, ar1_ov, ar12_ov, temp_ov, wprec_ov, cprec_ov, common.legend = T, nrow=1, legend="bottom")+
theme(plot.margin = margin(0,0,0,0, "cm"))+annotate("text", label="C. penicillatus", fontface="italic", size=14)
coef_ov1=coef_ov%>%gridExtra::grid.arrange(get_legend(coef_ov), heights = unit(c(70, 5), "mm"))
annotate_figure(coef_ov1, left = text_grob("C. penicillatus", face = "italic", size = 16, rot=90))
##Fig. 3: forecasts####
preds_plot1=ggplot()+
geom_line(data=ppcont_covs_dat, aes(y=abundance, x=newmoonnumber), size=0.75)+
geom_line(data=PPpreds_cont_same, aes(y=preds_same, x=moon, group=m,color="same"), alpha=0.4)+
geom_line(data=PPpreds_cont_switch, aes(y=preds_switch, x=moon, group=m, color="switched"), alpha=0.4)+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.title = element_blank())+ ylab("Abundance")+
ggtitle("control")+xlab("Sampling number (new moon number)")+xlab("")+
scale_colour_manual(values=c(same="darkblue", switched="darkred"), labels=c("non-transferred", "transferred"))
preds_plot1
annotate_figure(preds_plot1,left = text_grob("C. penicillatus", face = "italic", size = 16, rot=90))
preds_plot2=ggplot()+geom_line(data=ppexcl_covs_dat, aes(y=abundance, x=newmoonnumber), size=0.75)+
geom_line(data=PPpreds_excl_same, aes(y=preds_same, x=moon, group=m,color="same"), alpha=0.4)+
geom_line(data=PPpreds_excl_switch, aes(y=preds_switch, x=moon, group=m, color="switched"), alpha=0.4)+
theme_classic()+ylab("Abundance")+xlab("Sampling number (new moon number)")+
ggtitle("removal")+
scale_colour_manual(values=c(same="darkblue", switched="darkred"), labels=c("non-transferred", "transferred"))
preds_plot2
annotate_figure(preds_plot2,left = text_grob("C. penicillatus", face = "italic", size = 16, rot=90))
pp_predsplot=ggarrange(preds_plot1, preds_plot2, common.legend = T, ncol=1, nrow=2, legend="bottom")
annotate_figure(pp_predsplot,left = text_grob("C. penicillatus", face = "italic", size = 16, rot=90))
##Fig. 4: RMSE####
pp_h1=ggdensity(ppevals1, x="score_diff", color="plot",fill="plot",
palette=c("#9900cc", "#FFcc00"), ylab=F,rug=T, add="median",xlab=F, size=1,
main="h=1")+geom_vline(xintercept=0)+xlim(-25, 25)+
annotate("text", x=-15, y=0.40, label="non-transferred better fit")+
annotate("text", x=13, y=0.40, label="transferred better fit")+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(axis.text.y = element_blank())
pp_h6=ggdensity(ppevals6, x="score_diff", color="plot",fill="plot",
palette=c("#9900cc", "#FFcc00"), rug=T,ylab=F, add="median", xlab=F,size=1,
main="h=6")+geom_vline(xintercept=0)+xlim(-25, 25)+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(axis.text.y = element_blank())
pp_h12=ggdensity(ppevals12, x="score_diff", color="plot",fill="plot", size=1,
palette=c("#9900cc", "#FFcc00"),ylab=F, rug=T, add="median", xlab="RMSE difference (non-transferred - transferred)",
main="h=12")+geom_vline(xintercept=0)+xlim(-25, 25)+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(axis.text.y = element_blank())
pph=ggarrange(pp_h1,pp_h6,pp_h12, common.legend = T, nrow=3, legend="bottom")
annotate_figure(pph, top = text_grob("C. penicillatus", face = "italic", size = 16),
left= text_grob("RMSE", size=16, rot=90))
#Brier scores###
ppbr1=ggdensity(pp_briers1, x="brier_diff", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"),ylab=F, rug=T, add="median",xlab=F, size=1,
main="h=1")+geom_vline(xintercept=0)+xlim(-0.20, .20)+
annotate("text", x=-.12, y=70, label="non-transferred better fit")+
annotate("text", x=0.10, y=70, label="transferred better fit")+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(axis.text.y = element_blank())
ppbr2=ggdensity(pp_briers6, x="brier_diff", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), ylab=F, rug=T, add="median",xlab=F, size=1,
main="h=6")+geom_vline(xintercept=0)+xlim(-0.20, .20)+
theme(axis.text.y = element_blank())+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))
ppbr3=ggdensity(pp_briers12, x="brier_diff", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"),ylab=F, rug=T, add="median", size=1,
main="h=12", xlab="Brier score difference (non-transferred - transferred)")+
geom_vline(xintercept=0)+xlim(-0.2, 0.2)+
theme(axis.text.y = element_blank())+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))
ppbh=ggarrange(ppbr1,ppbr2,ppbr3, common.legend = T, nrow=3, legend="bottom")
annotate_figure(ppbh, top = text_grob("C. penicillatus", face = "italic", size = 16),
left= text_grob("Brier score", size=16, rot=90))
#PB results####
##Fig. 2: coefficient overlap####
pbint_ov=ggdensity(coef_df_PB, x="intercept", color="treatment",fill="treatment", ylab="Probability density",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab="intercept", size=1
)+geom_vline(xintercept=0)+
scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio = 1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
pbar1_ov=ggdensity(coef_df_PB, x="beta1", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab="AR (1)", size=1, ylab=F
)+geom_vline(xintercept=0)+
scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio = 1, axis.text.x = element_text(angle=90), axis.text.y = element_blank())
pbar12_ov=ggdensity(coef_df_PB, x="beta12", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab=" AR (12)", size=1, ylab=F
)+geom_vline(xintercept=0)+
scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio = 1, axis.text.x = element_text(angle=90), axis.text.y = element_blank())
pbtemp_ov=ggdensity(coef_df_PB, x="temp", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab="mean temperature (lag=1)", size=1, ylab=F
)+geom_vline(xintercept=0)+
scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio = 1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
pbwprec_ov=ggdensity(coef_df_PB, x="warm_precip", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab="warm precipitation", size=1, ylab=F
)+geom_vline(xintercept=0)+
scale_y_continuous(n.breaks=3)+
scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio = 1,axis.text.x = element_text(angle=90), axis.text.y = element_blank())
pbcprec_ov=ggdensity(coef_df_PB, x="cool_precip", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab="cool precipitation", size=1, ylab=F
)+geom_vline(xintercept=0)+
scale_x_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(aspect.ratio = 1, axis.text.x = element_text(angle=90), axis.text.y = element_blank())
pbcoef_ov=ggarrange(pbint_ov, pbar1_ov, pbar12_ov, pbtemp_ov, pbwprec_ov, pbcprec_ov, common.legend = T, nrow=1, legend="bottom")+
theme(plot.margin = margin(0,0,0,0, "cm"))+annotate("text", label="C. baileyi", fontface="italic", size=14)
pbcoef_ov1=pbcoef_ov%>%
gridExtra::grid.arrange(get_legend(pbcoef_ov), heights = unit(c(70, 5), "mm"))
annotate_figure(pbcoef_ov1, left = text_grob("C. baileyi", face = "italic", size = 16, rot=90))
##Fig. 3: forecasts####
pbpreds_plot1=ggplot()+
geom_line(data=pbcont_covs_dat, aes(y=abundance, x=newmoonnumber), size=0.75)+
geom_line(data=PBpreds_cont_same, aes(y=preds_same, x=moon, group=m,color="same"), alpha=0.4)+
geom_line(data=PBpreds_cont_switch, aes(y=preds_switch, x=moon, group=m, color="switched"), alpha=0.4)+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.title=element_blank())+
ggtitle("control")+xlab("")+ylab("Abundance")+
scale_colour_manual(values=c(same="darkblue", switched="darkred"), labels=c("non-transferred", "transferred"))
pbpreds_plot1
annotate_figure(pbpreds_plot1,left = text_grob("C. baileyi", face = "italic", size = 16, rot=90))
pbpreds_plot2=ggplot()+geom_line(data=pbexcl_covs_dat, aes(y=abundance, x=newmoonnumber), size=0.75)+
geom_line(data=PBpreds_excl_same, aes(y=preds_same, x=moon, group=m,color="same"), alpha=0.4)+
geom_line(data=PBpreds_excl_switch, aes(y=preds_switch, x=moon, group=m, color="switched"), alpha=0.4)+
theme_classic()+xlab("Sampling number (new moon number)")+
ylab("Abundance")+
ggtitle("removal")+
scale_colour_manual(values=c(same="darkblue", switched="darkred"), labels=c("non-transferred", "transferred"))
pbpreds_plot2
annotate_figure(pbpreds_plot2,left = text_grob("C. baileyi", face = "italic", size = 16, rot=90))
pb_predsplot=ggarrange(pbpreds_plot1, pbpreds_plot2, common.legend = T, ncol=1, nrow=2, legend="bottom")
annotate_figure(pb_predsplot,left = text_grob("C. baileyi", face = "italic", size = 16, rot=90))
##Fig. 4: RMSE####
pb_h1=ggdensity(pbevals1, x="score_diff", fill="plot", color="plot",
palette=c("#9900cc", "#FFcc00"), rug=T, ylab=F,add="median",xlab=F, size=1,
main="h=1")+geom_vline(xintercept=0)+xlim(-45, 45)+
annotate("text", x=-28, y=0.12, label="non-transferred better fit")+
annotate("text", x=26, y=0.12, label="transferred better fit")+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(axis.text.y = element_blank())
pb_h6=ggdensity(pbevals6, x="score_diff", fill="plot", color="plot",
palette=c("#9900cc", "#FFcc00"), rug=T, ylab=F,add="median",xlab=F, size=1, main="h=6")+geom_vline(xintercept=0)+xlim(-45, 45)+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+
theme(axis.text.y = element_blank())
pb_h12=ggdensity(pbevals12, x="score_diff", color="plot",fill="plot", size=1,ylab=F,
palette=c("#9900cc", "#FFcc00"), rug=T, add="median", xlab="RMSE difference (non-transferred - transferred)",
main="h=12")+geom_vline(xintercept=0)+xlim(-45, 45)+
theme(axis.text.y = element_blank())+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))
pbhh=ggarrange(pb_h1,pb_h6,pb_h12, common.legend = T, nrow=3, legend="bottom")
annotate_figure(pbhh, top = text_grob("C. baileyi", face = "italic", size = 16),
left= text_grob("RMSE", size=16, rot=90))
#Brier scores###
pbbr1=ggdensity(pb_briers1, x="brier_diff", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T,ylab=F, add="median",xlab=F, size=1,
main="h=1")+geom_vline(xintercept=0)+xlim(-0.4, 0.4)+
annotate("text", x=-0.25, y=13, label="non-transferred better fit")+
annotate("text", x=.25, y=13, label="transferred better fit")+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+theme(axis.text.y = element_blank())
pbbr2=ggdensity(pb_briers6, x="brier_diff", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T, add="median",xlab=F,ylab=F, size=1,
main="h=6")+geom_vline(xintercept=0)+xlim(-0.4, 0.4)+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+theme(axis.text.y = element_blank())
pbbr3=ggdensity(pb_briers12, x="brier_diff", color="treatment",fill="treatment",
palette=c("#9900cc", "#FFcc00"), rug=T,ylab=F, add="median", size=1,
main="h=12", xlab="Brier score difference (non-transferred - transferred)")+
geom_vline(xintercept=0)+xlim(-0.4, 0.4)+
scale_y_continuous(n.breaks=3, labels = scales::label_number(accuracy = 0.01))+theme(axis.text.y = element_blank())
pbbh=ggarrange(pbbr1,pbbr2,pbbr3, common.legend = T, nrow=3, legend="bottom")
annotate_figure(pbbh, top = text_grob("C. baileyi", face = "italic", size = 16),
left= text_grob("Brier score", size=16, rot=90))