-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGLBsynchrony.Rmd
1257 lines (1173 loc) · 41.8 KB
/
GLBsynchrony.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
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
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: "Code and data for: 'Mapping shifts in spatial synchrony in grassland birds to inform conservation planning'"
author: "Mike Allen and Julie Lockwood"
date: "December 1, 2020"
output: html_document
---
# Section 0. README
This code (i.e., the GLBsynchrony.Rmd file along with associated scripts "gridcollapse.func.R", "run.gamyet.func.R", and "mapsync.func.R") can be used to reproduce the analyses and figures from the associated manuscript (Allen & Lockwood 2020). It processes route-level North American Breeding Bird Survey count data summarized into 2x2 degree grid-cells to analyze spatial synchrony for 19 species of grassland birds. It relies heavily on data and code provided within the bbsBayes package (Edwards & Smith 2020; version 2.3.3.2020) and a Bayesian hierarchical model implemented within that package and described in Smith and Edwards (2020).
IMPORTANT NOTE: this code uses North American Breeding Bird Survey data (Pardieck et al. 2020) that must be accessed through the BBSbayes package (Edwards & Smith 2020). The code was most recently tested on 2020-10-27 with R version 4.0.2 (2020-06-22), RStudio version 1.3.1056, and bbsBayes version 2.3.3.2020. Version information for other packages are provided in the comments of Section 1 and full session info is provided in Section 8 below.
ANOTHER IMPORTANT NOTE: Many of the computations performed by this code can take a long time to complete depending on your computer (e.g., up to 15 hours per species - time period combination). Therefore, we provide the option to skip Sections 1-3 below by pre-loading synchrony measurements calculated from each of 3000 data sets resampled from the BBS index posterior values (i.e., the values used in the analysis presented in Allen & Lockwood 2020). Alternatively, you can skip all four 'processing' sections (1-4) and proceed directly to the data visualization section (Sections 5 & 6). This is possible because we also pre-load the final summarized spatial synchrony data in Section 1 as an Rdata file (file name: sync.summaries.Rdata).
## 0.1. Minimum list of files required for all code to work properly
File names with relative file paths. Note this does not include files that are included as a means of skipping sections as described above.
"GLBsynchrony.Rmd" - this document, which acts as a guide for all analyses.
"scripts/gridcollapse.func.R" - function to collapse data into 2x2 grid cells.
"scripts/run.gamyet.func.R" - function to run JAGS models to generate BBS indices.
"scripts/mapsync.func.R" - a function to create many resampled data sets from BBS index posteriors and compute synchrony metrics for each.
"scripts/gamyet.bug" - the heavy-tailed GAMYE JAGS model script from bbsBayes package (Edwards & Smith 2020, Smith & Edwards 2020).
"data/jags/db.csv" - a list of 1x1 latitude/longitude grid cells with area information for JAGS model.
"data/shapefiles/BCR_Terrestrial_master_International_clip.shp" - a shapefile for mapping Bird Conservation Regions (also need associated .shx, .dbf, .prj files).
"data/shapefiles/province.shp" - a shapefile for mapping Canadian provinces (also need associated .shx, .dbf, .prj files).
## 0.2. Table of contents
Sections
1. Load libraries, functions, data
2. Setup and run JAGS model to generate BBS indices
3. Computing mean spatial synchrony by species and grid cell
4. Summarize spatial synchrony data for mapping and analysis
5. Create figures in the manuscript body
6. Create figures and recreate analysis from Supporting Information
7. References
8. Package versions and session info
# Section 1. Load libraries, functions, data
## 1.1 Load packages, functions, & species data
We started with 30 grassland-specialist species listed in Sauer et al. (1995). Eleven species did not have enough data for synchrony analysis after processing (see Supporting Information in Allen & Lockwood 2020).
```{r, warning = F}
library(tidyverse) # tidyverse_1.3.0
library(bbsBayes) # bbsBayes_2.3.3.2020
library(rgdal) # version 1.5-16
library(maps) # version 3.3.0
library(ggthemes) # version 4.2.0
library(mapproj) # version 1.2.7
library(forcats) # version 0.5.0
source("scripts/gridcollapse.func.R") # function to prepare raw BBS data within 2 x 2 degree grid cells formatted for the Bayesian (JAGS) model (modified bbsBayes code)
source("scripts/run.gamyet.func.R") # function to run the hierarchical JAGS model to generate BBS indices for 2 x 2 degree grid cells (modified from bbsBayes code)
source("scripts/mapsync.func.R") # function to perform synchrony calculations
load("data/sync.summaries.Rdata") # load synchrony summaries to skip section 4
# load functions for "does not contain" and for identifying odd numbers
`%notin%` <-
Negate(`%in%`) # function for "does not contain"
is.odd <-
function(x)
x %% 2 != 0 # functions for identifying odd numbers
# Create species and time period list
species = data.frame(
species = c(
"Northern Harrier",
"Ring-necked Pheasant",
"Upland Sandpiper",
"Long-billed Curlew",
"Horned Lark",
"Sedge Wren",
"Sprague's Pipit",
"Dickcissel",
"Cassin's Sparrow",
"Vesper Sparrow",
"Lark Bunting",
"Savannah Sparrow",
"Baird's Sparrow",
"Grasshopper Sparrow",
"LeConte's Sparrow",
"Chestnut-collared Longspur",
"Bobolink",
"Eastern Meadowlark",
"Western Meadowlark",
"Northern Harrier",
"Ring-necked Pheasant",
"Upland Sandpiper",
"Long-billed Curlew",
"Horned Lark",
"Sedge Wren",
"Sprague's Pipit",
"Dickcissel",
"Cassin's Sparrow",
"Vesper Sparrow",
"Lark Bunting",
"Savannah Sparrow",
"Baird's Sparrow",
"Grasshopper Sparrow",
"LeConte's Sparrow",
"Chestnut-collared Longspur",
"Bobolink",
"Eastern Meadowlark",
"Western Meadowlark"
),
spcode = c(
"noha",
"rnep",
"upsa",
"lbcu",
"hola",
"sewr",
"sppi",
"dick",
"casp",
"vesp",
"larb",
"savs",
"bais",
"grsp",
"lcsp",
"cclo",
"bobo",
"eame",
"weme",
"noha",
"rnep",
"upsa",
"lbcu",
"hola",
"sewr",
"sppi",
"dick",
"casp",
"vesp",
"larb",
"savs",
"bais",
"grsp",
"lcsp",
"cclo",
"bobo",
"eame",
"weme"
),
startyr = c(rep(1968, 19), rep(1994, 19)),
endyr = c(rep(1993, 19), rep(2019, 19))
)
# note: this list excludes 11 species which had 0-4 final qualifying cells in the final synchrony change map (see Supporting Information in Allen & Lockwood 2020).
# Load spatial data for mapping
# Canada
can = rgdal::readOGR("data/shapefiles/province.shp")
canada.fort = ggplot2::fortify(can) %>%
filter(lat < 60.0001)
rm(can)
# load Bird Conservation Region shapefile
bcr <-
rgdal::readOGR("data/shapefiles/BCR_Terrestrial_master_International_clip.shp")
bcr.fort = ggplot2::fortify(bcr) %>%
filter(lat < 60.0001, long > -150)
rm(bcr)
# make a base map of North America
na <- ggplot() +
borders(
"world",
xlim = c(-150,-60),
ylim = c(35, 40),
colour = "gray85",
fill = "gray80"
) +
borders("state",
colour = "gray85",
fill = "gray80",
size = 0.5) +
geom_polygon(
data = canada.fort,
aes(x = long, y = lat, group = group),
colour = "gray85",
fill = "gray80",
size = 0.5
) +
ggthemes::theme_map() +
coord_map('albers', lat0 = 30, lat1 = 60)
na.blank <- ggplot() +
borders(
"world",
xlim = c(-150,-60),
ylim = c(25, 40),
colour = "gray80",
fill = "gray80"
) +
# borders("state", colour = "gray85", fill = "gray80") +
geom_polygon(
data = canada.fort,
aes(x = long, y = lat, group = group),
colour = "gray80",
fill = "gray80"
) +
ggthemes::theme_map() +
coord_map('albers', lat0 = 30, lat1 = 60)
rm(canada.fort)
```
# Section 2. Setup and run JAGS model to generate BBS indices
## 2.1. Access BBS data via bbsBayes package
You only need to run this code once. You'll need to accept the terms of use.
```{r}
#bbs_data <- fetch_bbs_data() # get BBS data (only need to run this once ever)
```
## 2.2. Format BBS data & export csv files that will be used in model
This code generates a csv file for each species and time period (saved to the data/jags folder) that will be used to run hierarchical models to generate BBS indices for 2x2 degree grid cells. The code, functions, and model (GAMYE-t) are from the bbsBayes package (Edwards & Smith 2020; see also Smith & Edwards 2020). It was necessary to generate a large number of CSV's, rather than running it all in one script, so that the hierarchical models could be run in the JagsUI package on a supercomputer cluster (on which we could not install bbsBayes). This code also customizes the bbsBayes code to stratify by 2x2 degree instead of 1x1 grid cells. This step may take half an hour or longer to run for all species.
```{r}
#####################################
#### a loop to create all species JAGS data files (as csv files)
#####################################
for (i in 1:nrow(species)) {
gridcollapse(
spcode = species$spcode[i],
species = species$species[i],
startyr = species$startyr[i],
endyr = species$endyr[i],
save.loc = "data/jags/" # where to write .csv files
)
}
#####################################
#### OR create JAGS data for one species at a time (as csv files)
#####################################
#gridcollapse(
# spcode = "bais",
# species = "Baird's Sparrow",
# startyr = 1994,
# endyr = 2019,
# save.loc = "data/jags/" # where to write .csv file
#)
```
## 2.3. Run Jags models to generate BBS indices
This code runs the Bayesian hierarchical model to generate BBS indices for each species/grid/year combination. The model is the General Additive Model with Year (heavy-tailed option) from the bbsBayes package (Edwards & Smith 2020, Smith & Edwards 2020). These can take ~1-15+ hrs for each species/period (running on 3 cores in parallel) depending on how widespread the species is (e.g., on the author's laptop, Baird's Sparrow took ~45 min, while Savannah Sparrow took 10-15 hours). We ran the models in parallel on the Rutgers Amarel supercomputer cluster (https://oarc.rutgers.edu). MCMC parameters can be customized by using the chains, iter, burn, and thin arguments (currently they are set at the bbsBayes defaults).
```{r}
#####################################
#### to run hierarchical model for a single species and period
#####################################
run.gamyet(
spcode = "bais",
species = "Baird's Sparrow",
startyr = 1968,
endyr = 1993,
chains = 3, # number of chains
burn = 20000, # number of burn-in iterations per chain
iter = 30000, # total number of iterations per chain including the burn-in
thin = 10, # keep every X iterations
par = T, # run chains in parallel?
gamyet.loc = "scripts/gamyet.bug", # location of the JAGS code file
data.loc = "data/jags/", # location of the JAGS data prepared in the previous step
save.posteriors = "data/index_posteriors/" # location to save the posterior samples to
)
# csv files are exported to the locations specified in "save.posteriors"
#####################################
#### OR use this code to run for all species and periods (not recommended as this would take a _very_ long time)
#####################################
#for(i in 1:nrow(species)){
# run.gamyet(spcode = species$spcode[i],
# species = species$species[i],
# startyr = species$startyr[i],
# endyr = species$endyr[i]
# )
#}
```
# Section 3. Computing mean spatial synchrony by species and grid cell
## 3.1 Compute mean spatial synchrony (r)
Setting other.metrics to "F" computes synchrony metrics (mean r) only; setting this to "T" also computes mean abundance and detrended variance for cells. The former is quicker. Setting uncertainty to "T" computes the metrics numerous times (set using "n.iter") on data sets created by resampling from the posterior distributions; setting it to "F" only computes them once using the median values of the posteriors. The latter is quicker. For reference, this function took ~15 minutes for a widespread species (Savannah Sparrow) on the author's laptop with all.methods = F and uncertainty = T. The resulting csv files of mean correlations saves to the location indicated by "save.to".
```{r}
#####################################
#### code to compute mean spatial synchrony for a single species and time period
#####################################
mapsync(
spcode = "bais",
startyr = 1968,
endyr = 1993,
n.iter = 3000, # number of data sets to create by resampling from the posterior distributions
n.post = 3000, # number of posterior samples that exist for each population estimate
uncertainty = T, # if TRUE, creates n.iter data sets resampled from the posteriors; if FALSE, creates one data set from posterior medians
maxzeros = 5, # sets the maximum number of zeros (in raw count data) to allow within each grid cell time series
other.metrics = T, # logical indicating whether to calculate mean abundance and mean variance in addition to mean synchrony
posterior.data.loc = "data/index_posteriors/",
jags.data.loc = "data/jags/",
save.to = "data/mapsync_iterations/" # where to save resulting .csv files
)
#####################################
#### OR compute mean spatial synchrony for all species and time periods
#### note: doing all species at once can take multiple hours to complete
#### note2: we ran widespread species on the cluster due to memory constraints
#####################################
#for(i in 1:nrow(species)){
# mapsync(
# spcode = species$spcode[i],
# startyr = species$startyr[i],
# endyr = species$endyr[i],
# n.iter = 3000,
# n.post = 3000,
# uncertainty = T,
# maxzeros = 5,
# other.metrics = T,
# posterior.data.loc = "data/index_posteriors/",
# jags.data.loc = "data/jags/",
# save.to = "data/mapsync_iterations/"
# )
#}
```
# Section 4. Summarize spatial synchrony data for mapping and analysis
## 4.1. Read in files with mean synchrony data (mean synchrony calculated from many resampled data sets from the BBS index posterior distributions)
Synchrony iterations were generated based on 3000 data sets from random draws of the BBS index posterior distributions (3000 for each species/grid cell/year). Computed mean r for each species/cell for each data set, resulting in a distribution of synchrony and synchrony difference estimates for each species/cell/year.
```{r}
# set folder location of "mapsync" files
mapsync.loc = "data/mapsync_iterations/"
# make a list of files containing the mean synchrony data sets
mapsync.files = sprintf(
paste0(mapsync.loc,"%s.mapsync.iterations.csv"),
paste0(
species$spcode,
".",
substr(species$startyr, 3, 4),
".",
substr(species$endyr, 3, 4)
)
)
# read in all mapsync data into list form
mapsync.iterations = lapply(mapsync.files, function(x)
read.csv(x))
names(mapsync.iterations) <- paste0(species$spcode,
".",
substr(species$startyr, 3, 4),
".",
substr(species$endyr, 3, 4))
# get a list of the grid cells in common between periods for each species (for filtering)
gridcom = function(x) {
Reduce(intersect, list(
unique(mapsync.iterations[[x]]$grid),
unique(mapsync.iterations[[x + length(unique(species$spcode))]]$grid)
))
}
numsp = length(unique(species$spcode))
grid.common = lapply(1:numsp, gridcom)
names(grid.common) = unique(species$spcode)
grid.common2 = c(grid.common, grid.common)
# subset only the grid cells in common between the two time periods
mapsync.com = list()
for (i in 1:nrow(species)) {
mapsync.com[[i]] = filter(mapsync.iterations[[i]], grid %in% grid.common2[[i]])
}
names(mapsync.com) <- species$spcode
# separate out period one, period two, and subtract the two to get change in synchrony
mapsync.com1 = mapsync.com[1:numsp]
mapsync.com2 = mapsync.com[(numsp+1):nrow(species)]
mapsync.dif1 = Map(function(listper2, listper1)
listper2[, 3:5] - listper1[, 3:5],
mapsync.com2,
mapsync.com1)
# add in iteration and grid variables back to synchrony change
mapsync.dif = Map(function(mainlist, listnames)
cbind(listnames[, 1:2], mainlist[, 1:3]),
mapsync.dif1,
mapsync.com1)
# add in species names
mapsync.dif.names = list()
for (i in 1:numsp) {
mapsync.dif.names[[i]] = cbind.data.frame(species = species$spcode[i], mapsync.dif[[i]])
}
rm(mapsync.files, mapsync.loc, grid.common, grid.common2, mapsync.iterations)
```
## 4.2. Average grid cells across species
```{r}
# create giant data frame with all iterations of dif, per1, and per2 for all spp
syncdif.all = cbind(
do.call("rbind", mapsync.dif.names)[, c(1:6)],
# difference
do.call("rbind", mapsync.com1)[, c(3:5)],
# period 1
do.call("rbind", mapsync.com2)[, c(3:5)] # period 2
)
colnames(syncdif.all) = c(
"spcode",
"iteration",
"grid",
"sync.dif",
"mean.dif",
"var.dif",
"sync.one",
"mean.one",
"var.one",
"sync.two",
"mean.two",
"var.two")
# average by grid cell across species for each iteration
syncdif = syncdif.all %>%
group_by(iteration, grid) %>%
summarise(
sync.dif = mean(sync.dif),
mean.dif = mean(mean.dif),
var.dif = mean(var.dif),
sync.one = mean(sync.one),
mean.one = mean(mean.one),
var.one = mean(var.one),
sync.two = mean(sync.two),
mean.two = mean(mean.two),
var.two = mean(var.two),
spp.per.grid = length(spcode)
) %>%
ungroup() %>%
filter(spp.per.grid >= 3) %>%
dplyr::select(-spp.per.grid)
# summarize synchrony change iterations by quantiles (all species)
syncdif.split = split(syncdif, f = syncdif$grid)
syncdif.sum.list = lapply(syncdif.split, function(x)
apply(x[, 3:11], 2L, FUN = quantile, probs = c(0.025, 0.5, 0.975)))
syncdif.sum = data.frame(do.call("rbind", syncdif.sum.list)) %>%
mutate(quant = rep(c("p2.5", "med", "p97.5"), length(syncdif.sum.list)),
grid = c(mapply(
unique(syncdif$grid),
FUN = function(x)
rep(x, 3)
))) %>%
tidyr::pivot_wider(
names_from = quant,
values_from = c(sync.dif:var.two)
) %>%
mutate(
lat = as.numeric(substr(grid, 1, 2)) - 0.5,
lon = as.numeric(substr(grid, 3, 10)) + 0.5
)
rm(mapsync.com, mapsync.com1, mapsync.com2, mapsync.dif1, mapsync.dif, mapsync.dif.names, syncdif.split, syncdif.sum.list)
```
## 4.3. Summarize data for individual species
```{r}
# join averaged species syncdif data to individual species data
syncdif.spp.prep = syncdif.all %>%
select(
spcode,
sync.dif,
mean.dif,
var.dif,
sync.one,
mean.one,
var.one,
sync.two,
mean.two,
var.two,
iteration
) %>%
bind_rows(
select(
syncdif,
sync.dif,
mean.dif,
var.dif,
sync.one,
mean.one,
var.one,
sync.two,
mean.two,
var.two,
iteration
)
) %>%
mutate(spcode = as.factor(case_when(is.na(spcode) ~ "All species",
TRUE ~ spcode)))
var.list = c(
"sync.dif",
"mean.dif",
"var.dif",
"sync.one",
"mean.one",
"var.one",
"sync.two",
"mean.two",
"var.two"
)
# average synchrony change of grid cells within species and iterations
# also correlations between synchrony and mean/variance/trends too
syncdif.spp = cbind(
aggregate(
syncdif.spp.prep[, var.list],
by = list(syncdif.spp.prep$iteration, syncdif.spp.prep$spcode),
mean
),
aggregate(
syncdif.spp.prep[, c("sync.dif")],
by = list(syncdif.spp.prep$iteration, syncdif.spp.prep$spcode),
FUN = function(x)
length(x)
)
)[, c(1:11, 14)] %>%
rename(iteration = Group.1,
spcode = Group.2,
n = x)
# get the median and quantiles of the mean synchrony change iterations
# for plotting, individual species and all species averaged
syncdif.spp.sum = cbind(
aggregate(
syncdif.spp[, var.list],
by = list(syncdif.spp$spcode),
FUN = function(x)
quantile(x, 0.5)
),
aggregate(
syncdif.spp[, var.list],
by = list(syncdif.spp$spcode),
FUN = function(x)
quantile(x, 0.025)
),
aggregate(
syncdif.spp[, var.list],
by = list(syncdif.spp$spcode),
FUN = function(x)
quantile(x, 0.975)
),
aggregate(syncdif.spp[, c("n")],
by = list(syncdif.spp$spcode),
mean)
)[, c(1:10, 12:20, 22:30, 32)]
colnames(syncdif.spp.sum) = c(
"spcode",
paste0(var.list, "_med"),
paste0(var.list, "_p2.5"),
paste0(var.list, "_p97.5"),
"n"
)
syncdif.spp.sum = syncdif.spp.sum %>%
mutate(
spcode = as.character(spcode),
fullsp = case_when(
spcode == "eame" ~ "E. Meadowlark",
spcode == "sppi" ~ "Sprague's Pipit",
spcode == "noha" ~ "N. Harrier",
spcode == "cclo" ~ "Chestnut-col. Longspur",
spcode == "grsp" ~ "Grasshopper Sparrow",
spcode == "upsa" ~ "Upland Sandpiper",
spcode == "rnep" ~ "Ring-necked Pheasant",
spcode == "vesp" ~ "Vesper Sparrow",
spcode == "bobo" ~ "Bobolink",
spcode == "lcsp" ~ "LeConte's Sparrow",
spcode == "hola" ~ "Horned Lark",
spcode == "savs" ~ "Savannah Sparrow",
spcode == "weme" ~ "W. Meadowlark",
spcode == "larb" ~ "Lark Bunting",
spcode == "dick" ~ "Dickcissel",
spcode == "bais" ~ "Baird's Sparrow",
spcode == "sewr" ~ "Sedge Wren",
spcode == "casp" ~ "Cassin's Sparrow",
spcode == "lbcu" ~ "Long-billed Curlew",
spcode == "All species" ~ "All species"
)
) %>%
mutate(
order = case_when(spcode == "All species" ~ -0.5,
TRUE ~ sync.dif_med),
fullsp = forcats::fct_reorder(fullsp, order, .desc = T)
)
# synchrony change of each grid cell (median of 3000 iterations) for each species
# for mapping
syncdif.spp.map = cbind(
aggregate(
syncdif.all[, var.list],
by = list(syncdif.all$spcode, syncdif.all$grid),
median
),
aggregate(
syncdif.all[, var.list],
by = list(syncdif.all$spcode, syncdif.all$grid),
FUN = function(x)
quantile(x, 0.025)
),
aggregate(
syncdif.all[, var.list],
by = list(syncdif.all$spcode, syncdif.all$grid),
FUN = function(x)
quantile(x, 0.975)
)
)[, c(1:11, 14:22, 25:33)]
colnames(syncdif.spp.map) = c(
"spcode",
"grid",
paste0(var.list, "_med"),
paste0(var.list, "_p2.5"),
paste0(var.list, "_p97.5")
)
syncdif.spp.map$lat = as.numeric(substr(syncdif.spp.map$grid, 1, 2)) - 0.5
syncdif.spp.map$lon = as.numeric(substr(syncdif.spp.map$grid, 3, 10)) +
0.5
syncdif.spp.map = syncdif.spp.map %>%
left_join(syncdif.spp.sum[, c("spcode", "fullsp")], by = "spcode")
syncdif.spp.map = syncdif.spp.map %>%
mutate(
fullsp2 = case_when(
as.character(fullsp) == "Chestnut-col. Longspur" ~ "C. Longspur",
as.character(fullsp) == "Ring-necked Pheasant" ~ "R. Pheasant",
TRUE ~ as.character(fullsp)
),
fullsp2 = forcats::fct_reorder(fullsp2, as.numeric(fullsp), .desc = F)
)
rm(syncdif.spp.prep, syncdif.spp, syncdif, syncdif.all, var.list)
```
# Section 5. Create figures in the manuscript body
# 5.1. Figure 2. Forest plot of synchrony change by species
```{r}
dotplot_n = syncdif.spp.sum %>%
dplyr::select(spcode, sync.dif_med, n, order) %>%
arrange(desc(order))
#x11(10,9)
syncdif.spp.sum %>%
mutate(fullsp = fct_reorder(fullsp, order, .desc = F)) %>%
ggplot() +
geom_vline(
aes(xintercept = 0),
color = "red",
linetype = 2,
size = 1.5
) +
geom_errorbarh(aes(
xmin = sync.dif_p2.5,
xmax = sync.dif_p97.5,
y = fullsp,
height = .2
)) +
geom_point(
aes(x = sync.dif_med, y = fullsp, fill = sync.dif_med),
size = 4.5,
pch = 21,
color = "black"
) +
scale_fill_distiller(name = "Mean\nchange",
palette = "RdYlBu",
limits = c(-0.105, 0.21)) +
labs(x = 'Change in spatial synchrony', y = "") +
annotate(
"text",
x = 0.48,
y = c(20.4, nrow(dotplot_n):1),
label = c("", paste0(dotplot_n$n[1:nrow(dotplot_n)])),
hjust = 1,
color = "gray25",
size = 4
) +
annotate(
"text",
x = 0.43,
y = nrow(dotplot_n),
# label="italic(n)", parse=TRUE,
label = "n = ",
color = "gray25",
fontface = "italic"
) +
annotate(
"text",
x = -0.3,
y = 1:20,
label = rev(
c(
"LC",
"LC",
"LC",
"VU",
"VU",
"LC",
"LC",
"LC",
"LC",
"LC",
"LC",
"LC",
"LC",
"LC",
"LC",
"NT",
"LC",
"LC",
"LC",
""
)
),
color = "gray25",
size = 3
) +
xlim(c(-0.3, 0.48)) +
theme_bw() +
theme(
axis.text = element_text(size = 12, color = "black"),
axis.title = element_text(size = 12, color = "black")
)
# ggsave("Fig_2_dotplot_color.tiff", dpi = 600)
```
## 5.2. Figure 3. Map change in synchrony for all species individually
```{r}
na +
geom_tile(aes(x = lon, y = lat, fill = sync.dif_med),
data = filter(syncdif.spp.map, sync.dif_med < .9)) +
# note: filter affects outlier high value in one cell for horned lark
scale_fill_distiller(
palette = "RdYlBu",
limits = c(-0.69, 0.58),
breaks = c(-0.5, 0, 0.5),
labels = c(-0.5, 0, 0.5)
) +
theme_classic() +
facet_wrap(~ fullsp2, ncol = 4) +
labs(fill = "Mean\nchange",
x = "",
y = "") +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
strip.background = element_blank(),
strip.text = element_text(size = 10),
legend.title.align = 0.5,
legend.position = "bottom",
legend.title = element_text(size = 8),
legend.text = element_text(size = 8),
axis.line = element_line(color = "transparent")
) +
geom_point(
aes(y = lat, x = lon),
alpha = 1,
color = "black",
shape = "+",
size = 1.5,
stroke = 2,
data = filter(syncdif.spp.map, sync.dif_p2.5 > 0)
) +
geom_point(
aes(y = lat, x = lon),
alpha = 1,
color = "black",
shape = "-",
size = 2.5,
stroke = 2,
data = filter(syncdif.spp.map, sync.dif_p97.5 < 0)
) +
geom_point(
aes(y = lat, x = lon),
alpha = 1,
color = "white",
shape = "-",
size = 1.5,
stroke = 2,
data = filter(syncdif.spp.map, sync.dif_p97.5 < 0)
) +
theme(
panel.spacing.y = unit(-0.5, "lines"),
panel.spacing.x = unit(-4, "lines"),
panel.background = element_rect(fill = "transparent", colour = NA),
plot.background = element_rect(fill = "transparent", colour = NA),
strip.text = element_text(vjust = -1.4)
)
ggsave("Fig_3.tiff", width = 7, height = 10, dpi = 600)
```
## 5.3. Figure 4. Map mean change in synchrony for all species combined
```{r}
#x11(13,9)
na.blank +
geom_tile(aes(x = lon, y = lat, fill = sync.dif_med),
data = syncdif.sum,
alpha = 1) +
scale_fill_distiller(name = "Mean\nchange", palette = "RdYlBu") +
theme_classic() +
theme(text = element_text(size = 15)) +
labs(x = "", y = "") +
theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank()
) +
theme(legend.title.align = 0.5,
legend.position = c(.85, .5)) +
theme(axis.line = element_line(color = "transparent")) +
geom_polygon(
data = bcr.fort,
aes(x = long, y = lat, group = group),
color = "gray60",
fill = "transparent",
size = .5
) +
geom_point(
aes(y = lat, x = lon),
alpha = 1,
color = "black",
shape = "+",
size = 4,
stroke = 1.25,
data = filter(syncdif.sum, sync.dif_p2.5 > 0)
) +
geom_point(
aes(y = lat, x = lon),
alpha = 1,
color = "white",
shape = "-",
size = 4,
stroke = 1.25,
data = filter(syncdif.sum, sync.dif_p97.5 < 0)
)
#ggsave("Figure_4.tiff", dpi = 600)
```
## 5.4. Raster map of mean synchrony change to make a shapefile for Geoda
This step is intermediate to creating the 'hotspot' polygon shapes in Figure 4b. The analysis was performed in the Geoda program (see manuscript) and the polygon map (pannel b) was made in QGIS.
```{r}
# Make a grid of spatial polygons to match up with meansyncdiff data
rr <- raster::raster(ext = raster::extent(-177, -52, 24, 69), res=c(2,2))
raster::values(rr) <- 1:raster::ncell(rr)
p <- raster::rasterToPolygons(rr)
attributes = tibble(lon = coordinates(p)[,1], lat = coordinates(p)[,2]) %>%
mutate(join = 1, grid = paste(trunc(lat),trunc(lon),sep="")) %>%
left_join(syncdif.sum, by="grid")
p$grid = attributes$grid
p$sync.dif_med = attributes$sync.dif_med
syncdif.poly = p[is.na(p$sync.dif_med)!=T,]
#writeOGR(syncdif.poly, "C:/Users/Mike/mike_files/Research/Chapter 3 - Synchrony/GLBsynchrony/data", "syncdif.poly", driver="ESRI Shapefile")
```
# Section 6. Create figures and recreate analysis from Supporting Information
## 6.1. Figure S2A. Map mean synchrony for all species individually - period 1
```{r}
na +
geom_tile(
aes(x = lon, y = lat, fill = sync.one_med),
#alpha = .7,
data = syncdif.spp.map
) +
scale_fill_distiller(palette = "RdYlBu", limits = c(-0.3, 0.75), breaks = c(-0.5,0,0.5), labels = c(-0.5,0,0.5)) +
theme_classic() +
facet_wrap( ~ fullsp2, ncol = 4) +
labs(
fill = "Mean\nsynchrony",
x = "",
y = ""#,
#title = "1968-1993 vs. 1994-2019"
) +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
strip.background = element_blank(),
strip.text = element_text(size = 10),
legend.title.align = 0.5,
legend.position = "bottom",
legend.title = element_text(size = 8),
legend.text = element_text(size = 8),
axis.line = element_line(color = "transparent")) +
geom_point(
aes(y = lat, x = lon),
alpha = .7,
color = "black",
shape = "+",
size = 1.5,
data = filter(syncdif.spp.map, sync.one_p2.5 > 0)
) +
geom_point(
aes(y = lat, x = lon),
alpha = .7,
color = "white",
shape = "-",
size = 1.5,
# stroke = 1.5,
data = filter(syncdif.spp.map, sync.one_p97.5 < 0)
) +
theme(panel.spacing.y = unit(-0.5, "lines"),
panel.spacing.x = unit(-4, "lines"),
panel.background = element_rect(fill = "transparent", colour = NA),
plot.background = element_rect(fill = "transparent", colour = NA),
strip.text=element_text(vjust=-1.4))
#ggsave("Figure_S2A.tiff", width = 7, height = 10, dpi = 600)
```
## 6.2. Figure S2B. Map mean synchrony for all species individually - period 2
```{r}
na +
geom_tile(
aes(x = lon, y = lat, fill = sync.two_med),
data = syncdif.spp.map
) +
scale_fill_distiller(palette = "RdYlBu", limits = c(-0.3, 0.75), breaks = c(-0.5,0,0.5), labels = c(-0.5,0,0.5)) +
theme_classic() +
facet_wrap( ~ fullsp2, ncol = 4) +
labs(
fill = "Mean\nsynchrony",
x = "",
y = ""
) +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
strip.background = element_blank(),
strip.text = element_text(size = 10),
legend.title.align = 0.5,
legend.position = "bottom",
legend.title = element_text(size = 8),
legend.text = element_text(size = 8),
axis.line = element_line(color = "transparent")) +
geom_point(
aes(y = lat, x = lon),
alpha = .7,
color = "black",
shape = "+",
size = 1.5,
# stroke = 1.5,
data = filter(syncdif.spp.map, sync.two_p2.5 > 0)
) +
geom_point(
aes(y = lat, x = lon),
alpha = .7,
color = "white",
shape = "-",
size = 1.5,
# stroke = 1.5,
data = filter(syncdif.spp.map, sync.two_p97.5 < 0)
) +
theme(panel.spacing.y = unit(-0.5, "lines"),
panel.spacing.x = unit(-4, "lines"),
panel.background = element_rect(fill = "transparent", colour = NA),
plot.background = element_rect(fill = "transparent", colour = NA),
strip.text=element_text(vjust=-1.4))
#ggsave("Figure_S2B.tiff", width = 7, height = 10, dpi = 600)
```
## 6.3. Figure S3A. Map change in mean abundance for all species individually
```{r}
na +
geom_tile(
aes(x = lon, y = lat, fill = 100*mean.dif_med/mean.one_med),
data = syncdif.spp.map
) +
scale_fill_distiller(palette = "RdYlBu", limits = c(-100, 100)) +
theme_classic() +
facet_wrap( ~ fullsp2, ncol = 4) +
labs(
fill = "% change\nabund.",
x = "",
y = ""#,
) +
theme(
text = element_text(size = 15),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
strip.background = element_blank(),
strip.text = element_text(size = 10),
legend.title.align = 0.5,
legend.position = "bottom",
legend.title = element_text(size = 8),
legend.text = element_text(size = 8),
axis.line = element_line(color = "transparent")) +