generated from EcologyR/templateRpackage
-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
193 lines (130 loc) · 10.7 KB
/
README.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
---
output: github_document
bibliography: references.bib
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# mappestRisk
<!-- badges: start -->
[](https://github.com/EcologyR/templateRpackage/actions/workflows/R-CMD-check.yaml) [](https://app.codecov.io/gh/EcologyR/templateRpackage?branch=master) `r badger::badge_lifecycle("experimental")` `r badger::badge_repostatus("WIP")` <!-- `r badger::badge_codefactor("ecologyr/templaterpackage")` -->
<!-- badges: end -->
The goal of `mappestRisk` package is to facilitate the transition from development data of pests obtained in lab-controlled conditions to understandable forecasts assessing risk of pest occurrence in a given region.
For that purpose, mappestRisk is built upon previous efforts such as `devRate` [@rebaudo2018], `rTPC` and `nls.multstart` packages [@padfield2021] and a methodology for predicting climatic suitability based on fundamental thermal niche as estimated by deductive, process-based approaches (i.e., those based on species performance variation across temperatures), as suggested in @taylor2019 . For now, `mappestRisk` is built for modelling developmental thermal performance curves, as this is the most commonly measured life-history trait in experimental approaches and it has been observed to have major contribution to fitness dependence on temperature [@pawar2024].
Therefore, mappestRisk has three different modules: *(1) model fitting & selection* using a set of the most commonly used equations describing developmental responses to temperature under the `nls.multstart`framework [@nls.multstart] using equation helpers from `rTPC`[@rTPC] and `devRate` [@devRate], with visualization of model fitting to help model selection by the user; (2) *calculation of suitability thermal limits,* which consist on a temperature interval delimiting the optimal performance zone or suitability; and (3) *climatic data extraction & visualization* with either exportable rasters or static or interactive map figures.
## Installation
```{r, warning = FALSE, message = FALSE}
# install.packages("devtools")
# devtools::install_github("EcologyR/mappestRisk")
#or alternatively
# remotes::install_github("EcologyR/mappestRisk")
#and load the package
#library(mappestRisk)
devtools::load_all() #for now, provisionally
```
If you want to clone or fork the repository or open and read some issues, you can find the code [here](https://github.com/EcologyR/mappestRisk).
## Example: `mappestRisk` workflow
### 1. Fit a thermal performance curve (TPC) to your data:
In this example, we'll show how to fit one to several thermal performance curves to a data set of development rate variation across temperatures[^1]. The following code provides an example as given in `fit_devmodels()` function documentation, with a data table showing the output of fitted models.
[^1]: At least 4 unique temperatures are required. Fore more details, see documentation and vignettes.
```{r, warning=FALSE}
data("aphid")
fitted_tpcs_aphid <- fit_devmodels(temp = aphid$temperature,
dev_rate = aphid$rate_value,
model_name = "all")
print(fitted_tpcs_aphid)
```
### 2. Plot the fitted TPCs and select the most appropriate:
To help select which model might be more appropriate, the function `plot_devmodels()` draws the predicted TPCs for each adequately-converged model. This step aims to improve model selection based not only on statistical criteria (i.e. AIC and number of parameters) but also on ecological realism, since curves can be graphically checked to select realistic shapes and thermal traits --vertical cuts with x-axis such as $CT_\min$, $CT_\max$ and $T_\text{opt}$ .
```{r, warning=FALSE}
plot_devmodels(temp = aphid$temperature,
dev_rate = aphid$rate_value,
fitted_parameters = fitted_tpcs_aphid,
species = "Brachycaudus schwartzi",
life_stage = "Nymphs")
```
### 3. Calculate predictions and bootstrap 100 TPCs for propagating parameter uncertainty
After careful examination of fitted TPCs in `plot_devmodels()`, the predictions of the selected models must be allocated in a `data.frame` in order to continue the package workflow towards climatic projection for pest risk analysis. Additionally, we recommend here to propagate uncertainty in parameter estimation of the fitted and selected TPC models using bootstrap procedures with residual resampling, following vignettes of `rTPC` package [@padfield2021]. This can be done with the function `predict_curves()` by setting the argument `propagate_uncertainty` to be `TRUE` and by providing a number of bootstrap sampling iterations in the argument `n_boots_samples`.
This function requires some time to compute bootstraps and predictions of the multiple TPCs it generates. Hence, it is important to use only one or few TPC models in the argument `model_name_2boot` based on visual examination from `plot_devmodels()` with reasonable selection. We discourage to bootstrap all TPC models for computational and time consuming reasons.
In our example, we choose three models that performed similarly and, unlike others with similar or even higher AICs, had a markedly higher left-skewness, e.g., *briere2, thomas* and *lactin2.* After a few minutes, we obtain the `data.frame`:
```{r, warning=FALSE}
preds_boots_aphid <-predict_curves(temp = aphid$temperature,
dev_rate = aphid$rate_value,
fitted_parameters = fitted_tpcs_aphid,
model_name_2boot = c("briere2", "thomas", "lactin2"),
propagate_uncertainty = TRUE,
n_boots_samples = 100)
```
### 4. Plot uncertainty TPCs:
The bootstrapped uncertainty curves can be plotted easily with the `plot_uncertainties()` function.
```{r, message=FALSE, warning=FALSE}
plot_uncertainties(bootstrap_uncertainties_tpcs = preds_boots_aphid,
temp = aphid$temperature,
dev_rate = aphid$rate_value,
species = "Brachycaudus schwartzi",
life_stage = "Nymphs")
```
### 5. Calculate thermal suitability bounds:
Once a model have been selected under both ecological realism and statistical criteria, the user can estimate the thermal boundaries defining the suitable range of temperatures for the studied population. The `thermal_suitability_bounds()` function calculate these values given the `tibble` output from `predict_curves()` function and the selected model name (note that only one model is allowed each time for this function). Additionally, a value of suitability defining the quantile-upper part of the curve can be provided by the user ($\text{Q}_{75}$ by default). If the user has propagated uncertainty in `predict_curves()` function, this function inherits the bootstrapped TPCs to calculate thermal suitability boundaries for each bootstrapped curved in addition to the estimated curve.
```{r, warning=FALSE}
boundaries_aphid <- therm_suit_bounds(preds_tbl = preds_boots_aphid,
model_name = "lactin2",
suitability_threshold = 80)
print(boundaries_aphid)
```
### 4. Climatic data extraction and projection
Using the thermal boundaries provided by the previous function and a set of raster maps of monthly temperatures for a given region (which can be provided by the user or downloaded by the function), two maps are plotted. The left panel (*risk map*) shows how many months a year (on average across simulations) thermal conditions are suitable for the development of the pest. The right panel (*uncertainty map*) shows the standard deviation of the risk index for each pixel.
Independently of whether `plot` is set to `TRUE` or `FALSE`, the function returns a raster that can be later exported by the user.
```{r}
# Downloading temperature data with geodata::wordlclim_global().
# May take time the first time you use the function on the same 'path'.
risk_rast <- map_risk(t_vals = boundaries_aphid,
path = "~/downloaded_maps", # directory to download data
region = "Morocco",
mask = TRUE,
plot = TRUE,
interactive = FALSE,
verbose = TRUE)
# We can also save the raster with:
# terra::writeRaster(risk_rast, filename = "~/output_maps/risk_rast.tif")
# Alternativele 1: if you already have a raster of monthly average temperatures for your region of interest, you can use that as input:
## load it (here Luxembourg data)
tavg_file <- system.file("extdata/tavg_lux.tif", package = "mappestRisk")
## import it with `terra`
tavg_rast <- terra::rast(tavg_file)
## and apply the function
risk_rast <- map_risk(t_vals = boundaries_aphid,
t_rast = tavg_rast,
mask = TRUE,
plot = TRUE,
interactive = FALSE,
verbose = TRUE)
# Alternative 2: if the region to project risk is given by the user as a spatial feature (sf) object.
andalucia_sf <- readRDS(system.file("extdata", "andalucia_sf.rds",
package = "mappestRisk"))
risk_rast <- map_risk(t_vals = boundaries_aphid,
path = "~/downloaded_maps", # directory to download data
region = andalucia_sf,
mask = TRUE,
plot = TRUE,
interactive = FALSE,
verbose = TRUE)
```
### 5. Interactive map
Additionally, the `map_risk()` can produce interactive maps based on `terra::plet()` by setting the argument `interactive` to `TRUE`. These maps are html objects that can be displayed in your browser, with one clickable tab for the risk map and another one for the uncertainty map.
You can try yourself this option in your RStudio Viewer.
### Citation
If using this package, please cite it:
```{r comment=NA}
citation("mappestRisk")
```
## Funding
The development of this software has been funded by Fondo Europeo de Desarrollo Regional (FEDER) and Consejería de Transformación Económica, Industria, Conocimiento y Universidades of Junta de Andalucía (proyecto US-1381388 led by Francisco Rodríguez Sánchez, Universidad de Sevilla).
## References: