abmR: An R Package for Agent-based Model Analysis of Large-scale Movements Across Taxa
Benjamin Gochanour, Javier Fernandez Lopez, Andrea Contina
2021-04-28
1 Getting Started
1.1 Installation
To use abmR, you must first install it from Github using devtools and load the library:
# devtools: install_github("bgoch5/abmR")
# If install gives errors, try running the following:
# Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="true")
library(abmR,quietly=TRUE,warn.conflicts=FALSE)While this package is still in development, it will be updated frequently, so please be sure to re-install frequently. Installing abmR will also automaticaly install its dependencies, if you don’t already have them installed. These include raster, sp, rgdal, table1, googledrive, swfscMisc, geosphere, kableExtra, and gtsummary, and ggplot.
1.2 Requirements for Environmental Data Rasters and Species objects
The modeling functions discussed in section 2 require two input objects to work, which we will discuss here.
1.2.1 Environmental raster stack
Example data that will be used in this Vignette and other package documentation can be downloaded using get_ex_data(). This function will download to your hard drive (to a location specified by your current working directory) the below NDVI raster files, totaling approximately 620 MB. The first time you use this function, you will be directed to your browser and required to sign in to your Google account to connect to the Tidyverse API. If you use the function a second time, you may simply follow the prompts and enter a number corresponding to the previous accounts listed.
# get_ex_data()Once you have the data downloaded, you must read it into R using commands like the following.
You may instead use your own environmental raster data - NDVI or otherwise - and read it into R in a similar way. Please refer to the file “Aggregating Environmental Rasters.R” on the abmR Github page for an example on how you might create a final rasterstack from a group of NOAA NDVI .nc files.
ndvi_raster_EU=stack("C:/Users/BGOCHANOUR/Documents/GitHub/move-model/Data/2013 NDVI/NDVI_2013_Europe")
ndvi_raster_EU_composite=raster("C:/Users/BGOCHANOUR/Documents/GitHub/move-model/Data/2013 NDVI/NDVI_2013_Europe_composite")We can now quickly examine our raster data to see if it read in correctly.
The first raster is a 27-layer stack (of which we will plot the first four layers), and the second is a composite raster formed by taking the mean of all layers.
plot(ndvi_raster_EU[[1:4]]) # First four layers of 27 layer RasterStack
plot(ndvi_raster_EU_composite)
1.2.2 Species object
This object will be constructed with the function as.species, and contains information on origin longitude (x) and latitude (y) as well as two morphological parameters, morphpar1 and morphpar2 for your agent. x and y are the only required parameters. However, if you choose to specify one of morphpar1 through morphpar2sign, you must also specify all other parameters. See ?as.species for more details.
EU.pop = as.species(x=10,
y=50,
morphpar1=24,
morphpar1mean=16,
morphpar1sd=5,
morphpar1sign="Pos",
morphpar2=12,
morphpar2mean=9,
morphpar2sd=3,
morphpar2sign="Neg")2 Running Agent-based models
abmR has two functions for modeling animal movement: moveSIM and energySIM. Both moveSIM and energySIM will output two things: (1) a dataframe of results ($results) and (2) a summary table of input parameters used ($run_params). For more details on how each function works, please see the documentation by using ?moveSIM or ?energySIM.
2.1 moveSIM
First, we will look at moveSIM, which runs basic Brownian/Ornstein-Uhlenbeck agent-based models for multiple replicates. Here, mortality occurs when an agent crosses fail_thresh on n_failures+ 1 consecutive timesteps. For more function details, please consult the documentation: ?moveSIM.
For each timestep, agents can have status “Alive”, “Stopped”, or “Died”. All agents start alive and may stop if, on a particular timestep, there are no non-NA raster values in the search region. This often occurs when agents are searching over an ocean or a large lake, for example. Once an agent stops, they remain stopped for the rest of the run. Similarly, once an agent dies, they retain this status for all subsequent timesteps. All timesteps with agent status “Stopped” or “Died” will have lat/lon=NA, so as to not affect subsequent analyses.
set.seed(1)
EU.move=moveSIM(replicates = 50,
days=27,
env_rast=ndvi_raster_EU,
search_radius=250,
sigma=0.1,
dest_x=999,
dest_y=999,
mot_x=0.7,
mot_y=0.7,
modeled_species=EU.pop,
optimum = 0.7,
direction="S",
write_results=F,
single_rast=F,
mortality = T,
n_failures = 3,
fail_thresh = 0.50)
#> [1] "Agent died"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Number of agents processed: 5"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 10"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 15"
#> [1] "Agent died"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 20"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Agent died"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 25"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 30"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Number of agents processed: 35"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 40"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 45"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 50"We can simply view our results in the console by typing NA.move or EU.move. However, we have created the function tidy_results() to provide cleaner output. Simply supply your moveSIM() or energySIM() output and specify type=“results” or type=“run_params”, depending on which portion of your output you’d like to display.
Let’s examine our moveSIM() output. Here, we view just the first 10 rows.
tidy_results(EU.move,
type="results",
nrows=10)| agent_id | day | lon | lat | curr_status | distance | |
|---|---|---|---|---|---|---|
| 1 | Agent_01 | 0 | 10.000 | 50.000 | Alive | NA |
| 2 | Agent_01 | 1 | 8.625 | 48.375 | Alive | 206.705 |
| 3 | Agent_01 | 2 | 8.375 | 45.875 | Alive | 278.942 |
| 4 | Agent_01 | 3 | 9.375 | 45.125 | Alive | 114.271 |
| 5 | Agent_01 | 4 | 8.275 | 44.275 | Alive | 128.562 |
| 6 | Agent_01 | 5 | 8.075 | 44.175 | Alive | 19.454 |
| 7 | Agent_01 | 6 | 4.875 | 43.775 | Alive | 260.174 |
| 8 | Agent_01 | 7 | 1.975 | 42.775 | Alive | 260.049 |
| 9 | Agent_01 | 8 | -0.025 | 42.675 | Alive | 163.929 |
| 10 | Agent_01 | 9 | 2.025 | 42.575 | Alive | 168.278 |
tidy_results(EU.move,
type="run_params")| Param value | |
|---|---|
| replicates | 50 |
| days | 27 |
| env_raster | env_raster |
| search_radius | 250 |
| sigma | 0.1 |
| dest_x | 999 |
| dest_y | 999 |
| mot_x | 0.7 |
| mot_y | 0.7 |
| modeled_species | EU.pop |
| optimum | 0.7 |
| n_failures | 3 |
| direction | S |
| mortality | TRUE |
| write_results | FALSE |
| single_rast | FALSE |
| missing_pct | 63.9285714285714 |
| mortality_pct | 18 |
2.2 energySIM
Next, we turn our attention to energySIM(), the more advanced function that simulates agent energy stores. Agents begin with energy=init_energy, and die when they reach energy=0. energy_adj mediates the energy gain or loss for each timestep. The 11 elements of this vector represent energy gain/loss for environmental values in the optimal range (1st element) and within 10, 20, …, 80, 90, and 90+ percent deviation (11th element) of the optimal range. We recommend using the default values for energy_adj unless you have a reason for changing them.
For each timestep, agents can have status “Alive”, “Stopped”, or “Died”. All agents start alive and may stop if, on a particular timestep, there are no non-NA raster values in the search region. This often occurs when agents are searching over an ocean or a large lake, for example. Once an agent stops, they remain stopped for the rest of the run. Similarly, once an agent dies, they retain this status for all subsequent timesteps. All timesteps with agent status “Stopped” or “Died” will have lat/lon=NA, so as to not affect subsequent analyses.
EU.energy = energySIM(replicates = 50,
days=27,
env_rast=ndvi_raster_EU,
search_radius=250,
sigma=0.1,
dest_x=999,
dest_y=999,
mot_x=0.7,
mot_y=0.7,
modeled_species=EU.pop,
optimum_lo = 0.65,
optimum_hi=0.75,
init_energy=100,
direction="S",
write_results=F,
single_rast=F,
mortality = T)
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 5"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 10"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 15"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 20"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Number of agents processed: 25"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 30"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 35"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Agent died"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 40"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 45"
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Can't find any non-NA cells. Searching over lake or ocean. Agent stopped."
#> [1] "Number of agents processed: 50"The tidy_results function will also work on energySIM() output. Here, we view just the first 10 rows of our results.
tidy_results(EU.energy,
type="results",
nrows=10)| agent_id | day | lon | lat | curr_status | energy | delta_energy | distance | |
|---|---|---|---|---|---|---|---|---|
| 1 | Agent_01 | 0 | 10.000 | 50.000 | Alive | 100 | NA | NA |
| 2 | Agent_01 | 1 | 8.875 | 48.125 | Alive | 85 | -15 | 224.268 |
| 3 | Agent_01 | 2 | 8.225 | 45.875 | Alive | 65 | -20 | 255.281 |
| 4 | Agent_01 | 3 | 8.125 | 44.625 | Alive | 90 | 25 | 139.370 |
| 5 | Agent_01 | 4 | 10.775 | 44.225 | Alive | 70 | -20 | 215.321 |
| 6 | Agent_01 | 5 | 11.425 | 44.075 | Alive | 90 | 20 | 54.537 |
| 7 | Agent_01 | 6 | 11.625 | 43.725 | Alive | 100 | 10 | 42.135 |
| 8 | Agent_01 | 7 | NA | NA | Stopped | NA | NA | NA |
| 9 | Agent_01 | 8 | NA | NA | Stopped | NA | NA | NA |
| 10 | Agent_01 | 9 | NA | NA | Stopped | NA | NA | NA |
tidy_results(EU.energy,
type="run_params")| Param value | |
|---|---|
| replicates | 50 |
| days | 27 |
| env_raster | env_raster |
| search_radius | 250 |
| sigma | 0.1 |
| dest_x | 999 |
| dest_y | 999 |
| mot_x | 0.7 |
| mot_y | 0.7 |
| modeled_species | EU.pop |
| optimum_lo | 0.65 |
| optimum_hi | 0.75 |
| init_energy | 100 |
| direction | S |
| mortality | TRUE |
| write_results | FALSE |
| single_rast | FALSE |
| missing_pct | 63.1428571428571 |
| mortality_pct | 4 |
3 Visualizing/Summarizing Results
The output from moveSIM and energySIM are simple dataframes, so one may easily visualize their results using custom codes. However, we have provided moveVIZ and energyVIZ, respectively, to provide pre-prepared solutions for visualizing and summarizing your movement simulation results.
3.0.1 Movement tracks
moveVIZ(EU.move,
type="plot",
title="moveSIM European results")
energyVIZ(EU.energy,
type="plot",
title="EnergySIM European results")
3.0.2 Summary table
moveVIZ(EU.move,
type="summary_table")| Characteristic | N = 1,4001 |
|---|---|
| lon | 8.4 (3.2, 10.0) |
| Unknown | 895 |
| lat | 44.63 (43.38, 45.98) |
| Unknown | 895 |
| curr_status | |
| Alive | 505 (36%) |
| Died | 132 (9.4%) |
| Stopped | 763 (55%) |
| distance | 204 (113, 241) |
| Unknown | 945 |
|
1
Median (IQR); n (%)
|
|
energyVIZ(EU.energy,
type="summary_table")| Characteristic | N = 1,4001 |
|---|---|
| lon | 8.3 (4.9, 10.0) |
| Unknown | 884 |
| lat | 44.60 (43.41, 46.03) |
| Unknown | 884 |
| curr_status | |
| Alive | 514 (37%) |
| Died | 31 (2.2%) |
| Stopped | 855 (61%) |
| energy | 85 (75, 100) |
| Unknown | 884 |
| delta_energy | 0 (-15, 10) |
| Unknown | 934 |
| distance | 155 (87, 223) |
| Unknown | 934 |
|
1
Median (IQR); n (%)
|
|
The next two types of output are available only for energySIM() output.
3.0.3 Stratified tables
energyVIZ(EU.energy,
type="strat_table")| -25 (N=5) |
-20 (N=71) |
-15 (N=61) |
-10 (N=41) |
-5 (N=24) |
0 (N=89) |
5 (N=46) |
10 (N=44) |
15 (N=38) |
20 (N=30) |
25 (N=17) |
Overall (N=1400) |
|
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| energy | ||||||||||||
| Mean (SD) | 63.0 (16.8) | 68.1 (13.2) | 71.1 (18.9) | 65.1 (27.1) | 77.3 (21.7) | 93.1 (16.1) | 83.5 (18.6) | 87.7 (15.3) | 84.9 (16.3) | 84.7 (12.7) | 90.0 (8.66) | 82.0 (19.9) |
| Median [Min, Max] | 75.0 [40.0, 75.0] | 65.0 [10.0, 80.0] | 80.0 [10.0, 85.0] | 75.0 [0, 90.0] | 82.5 [10.0, 95.0] | 100 [20.0, 100] | 85.0 [15.0, 100] | 90.0 [25.0, 100] | 85.0 [25.0, 100] | 85.0 [45.0, 100] | 90.0 [70.0, 100] | 85.0 [0, 100] |
| Missing | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 884 (63.1%) |
| day | ||||||||||||
| Mean (SD) | 13.0 (3.24) | 2.52 (2.44) | 4.75 (6.22) | 9.27 (6.38) | 10.5 (7.02) | 9.35 (5.07) | 9.59 (6.96) | 6.98 (5.35) | 7.21 (4.85) | 4.30 (1.49) | 4.47 (2.21) | 13.5 (8.08) |
| Median [Min, Max] | 13.0 [8.00, 17.0] | 2.00 [1.00, 14.0] | 2.00 [1.00, 27.0] | 8.00 [1.00, 23.0] | 9.00 [2.00, 27.0] | 8.00 [3.00, 25.0] | 6.50 [3.00, 26.0] | 5.00 [3.00, 26.0] | 5.00 [3.00, 23.0] | 4.00 [3.00, 8.00] | 3.00 [3.00, 10.0] | 13.5 [0, 27.0] |
| distance | ||||||||||||
| Mean (SD) | 171 (118) | 215 (47.8) | 193 (72.4) | 144 (86.2) | 147 (85.3) | 156 (85.2) | 137 (67.2) | 116 (65.1) | 113 (51.3) | 97.9 (56.0) | 93.4 (44.6) | 153 (79.0) |
| Median [Min, Max] | 166 [38.5, 355] | 227 [9.96, 258] | 218 [16.7, 357] | 148 [16.7, 311] | 163 [5.57, 266] | 152 [11.8, 346] | 146 [27.4, 265] | 91.6 [24.0, 275] | 118 [29.8, 267] | 86.9 [16.3, 238] | 93.6 [6.84, 166] | 155 [5.57, 357] |
| Missing | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 934 (66.7%) |
3.0.4 Gradient plots
energyVIZ(EU.energy,
type="gradient")
#> [inverse distance weighted interpolation]
4 Getting Help
Need help? Have suggestions to make abmR better? If so, please open an issue or pull request on Github (https://github.com/bgoch5/abmr), or drop me an email at ben.gochanour@ou.edu.