-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathREADME.Rmd
341 lines (233 loc) · 11.2 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
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
---
title: "Hands on with {epicrop}"
author: "Adam H. Sparks"
date: "11/03/2024"
output: github_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
This will install {epicrop} from Codeberg if you do not already have it installed.
You can click "Knit" at the top of the window to see the final document or browse the file on GitHub, <https://github.com/openplantpathology/hands_on_with_epicrop>.
```{r, eval=FALSE}
if (!require("remotes"))
install.packages("remotes")
remotes::install_git("https://codeberg.org/adamhsparks/epicrop",
build_vignettes = TRUE
)
```
We have seen how the "EPIRICE" model works and discussed the theory behind the SEIR model that drives it.
Now we will work with the R package and examine some of these different factors to gain a better understanding of how the model works.
## The SEIR() function
The `SEIR()` function drives the EPIRICE model.
You can view the code on GitHub in the {epicrop} repository, <https://codeberg.org/adamhsparks/epicrop/src/branch/main/R/SEIR.R>.
If we look at the help file for `SEIR()` it will tell us how to use it.
```{r help_seir}
library("epicrop")
?SEIR()
```
### The parameters
* `wth`
The first parameter, `wth`, is a data.frame object of weather to drive the model.
{epicrop} has relatively simple needs for weather to run the model.
Daily temperature, relative humidity and rainfall are all that are required.
{epicrop} has a built-in function, `get_wth()` that will automatically fetch and format the weather data for use in `SEIR()` (or any of the other specialised functions for modelling disease) from the NASA POWER API, <https://power.larc.nasa.gov>.
* `emergence`
The date that the crop emerges, or if transplanted rice, is transplanted.
* `onset`
The expected number of days until the onset of disease after emergence date (day, integer).
Described in Table 1 of Savary _et al._ 2012.
* `duration`
Simulation duration _i.e._, growing season length (day, integer).
Described in Table 1 of Savary _et al._ 2012.
* `rhlim`
Relative humidity value threshold to decide whether leaves are wet or not (numeric).
Savary _et al._ 2012 used 90%.
* `rainlim`
Rainfall amount (mm) threshold to decide whether leaves are wet or not (numeric).
Savary _et al._ 2012 used 5mm.
* `H0`
Initial number of plant's healthy sites (integer).
Described in Table 1 of Savary _et al._ 2012.
* `I0`
Initial number of infective sites (integer).
Described in Table 1 of Savary _et al._ 2012.
* `RcA`
Modifier for _Rc_ (the basic infection rate corrected for removals) for crop age (numeric vector).
Described in Table 1 of Savary _et al._ 2012.
* `RcT`
Modifier for _Rc_ (the basic infection rate corrected for removals) for temperature (numeric vector).
Described in Table 1 of Savary _et al._ 2012.
* `RcOpt`
Potential basic infection rate corrected for removals (numeric).
Derived from Table 1 of Savary _et al._ 2012.
* `i`
Duration of infectious period (day, integer).
Described in Table 1 of Savary _et al._ 2012.
* `p`
Duration of latent period (day, integer).
Described in Table 1 of Savary _et al._ 2012.
* `Sx`
The maximum number of sites (integer).
Described in Table 1 of Savary _et al._ 2012.
* `a`
The aggregation coefficient, values are from 1 to >1 (numeric).
Described in Table 1 of Savary _et al._ 2012.
* `RRS`
Relative rate of physiological senescence (numeric).
Described in Table 1 of Savary _et al._ 2012.
* `RRG`
Relative rate of growth (numeric).
Described in Table 1 of Savary _et al._ 2012.
## The "predict" family of functions
Currently {epicrop} provides five pre-parameterised functions to predict rice diseases.
These functions are the EPICROP model as published by Savary _et al._ 2012.
* `predict_bacterial_blight()`,
* `predict_brown_spot()`,
* `predict_leaf_blast()`,
* `predict_sheath_blight()`, and
* `predict_tungro()`
You can view the help file for any of them using the same as above, `?predict_bacterial_blight`.
## Using the "predict" functions
I have attempted to make it as simple as possible to use these functions and generate model outputs.
Modelling and plotting the disease progress over time for bacterial leaf blight is as follows.
```{r plot-bb-no-eval, eval=TRUE}
# get weather for IRRI Zeigler Experiment Station for the year 2000
wth <- get_wth(
lonlat = c(121.25562, 14.6774),
dates = c("2000-01-01", "2000-12-31")
)
bb_wet <- predict_bacterial_blight(wth, emergence = "2000-07-01")
plot(x = bb_wet$dates, y = bb_wet$intensity, type = "l")
```
## Differences in the five functions
If you look at the code for the five functions, you will see that they each have their own values for the parameters found in `SEIR()`.
All of the functions can be found in the GitHub repository, <https://github.com/adamhsparks/epicrop/tree/main/R>.
We will look at `predict_bacterial_blight()` first.
On [line 92](https://github.com/adamhsparks/epicrop/blob/b3c85b997180c3f90e67cc02a5e58d62491d9edd/R/predict_bacterial_blight.R#L92) we see that the `age_coef_rc` is set up as:
```{r bb-age}
age_coef_rc <-
cbind(
0:12 * 10,
c(1, 1, 1, 0.9, 0.62, 0.43, 0.41, 0.42, 0.41, 0.41, 0.41, 0.41, 0.41)
)
```
and the `temp_coef_rc` is set up as:
```{r bb-temp}
temp_coef_rc <-
cbind(
16 + (0:8 * 3),
c(0, 0.29, 0.44, 0.90, 0.90, 1.0, 0.88, 0.01, 0)
)
```
These values are used in the `RcA` and `RcT` parameters of `SEIR()` to modify the disease's progress.
I find it easier to think about how this affects the model by visualising what this is creating.
```{r plot-bb-age}
plot(x = age_coef_rc[, 1], y = age_coef_rc[, 2], main = "BB age_coef_rc")
```
```{r plot-bb-temp}
plot(x = temp_coef_rc[, 1], y = temp_coef_rc[, 2], main = "BB temp_coef_rc")
```
### Exercise 1
Now, you compare these curves with the curves for leaf blast and answer the following questions.
1. How do they differ?
2. What does this imply about the pathogens and diseases that they cause?
Now compare one of the other disease models' curves with these two and see how they differ.
All of the code for the disease models are available [here](https://github.com/adamhsparks/epicrop/tree/main/R).
## RcOpt
We saw `RcOpt` in the presentation.
It is the reference potential value of the basic infection rate corrected for removals for each disease.
The default value for `predict_bacterial_blight()` is 0.87 while for `predict_leaf_blast()` it is 1.14.
### Exercise 2
Discuss how this value changing affects the disease progress.
What does a higher or lower value for `RcOpt` mean?
## H0, I0 and Sx
Consider the `H0` and `I0` parameters for the different diseases.
The model works for several types of disease, leaf area, whole leaf, tiller, whole plant all in a 1m^2 area.
### Exercise 3
1. What are `H0`, `I0` and `Sx`?
2. How are they related?
3. Consider how these values might differ between the functions/diseases.
1. Explain how/why these values differ.
2.
a. Are there other scales or plant organs that you would consider if you were adding new diseases?
b. What disease would these be (does not have to be rice)?
c. What might the default value be for the new disease?
## p
We did not see `p`, _per se_ in the presentation but we did cover the idea of latent and infectious sites.
`p` controls the latent period, that is how long a site remains in the `E` (exposed) state of SEIR before becoming infectious.
### Exercise 4
Compare and contrast the different values for `p` from each of the `predict` functions.
1. What does this tell you about the pathogens' life-cycle and reproductive strategies?
2. Which one has the shortest latent period?
3. Which one has the longest latent period?
## RRS and RRG
The last parameters in the model are `RRS` and `RRG`, or the relative rate of senescence and the relative rate of growth.
Both of these parameters relate to the crop itself.
Compare the values found in the different `predict` functions.
1. Which are different and which are the same?
2. How would changing these values affect the rate of disease progress?
# Working with {epicrop} outputs
`SEIR()` produces data tables that start on day 0, establishment or emergence and progress to the end of the season, `duration`.
By default in the `predict` functions, this is 120 days for a 120 day rice crop.
There are 16 possible columns, 14 will always be returned with `lat` and `lon` optionally returned if provided in the `wth` weather input object.
Any values calculated by the model are available at the end of the run.
## Description of fields and values from SEIR() output
1 **simday** Zero indexed day of simulation that was run
2 **dates** Date of simulation
3 **sites** Total number of sites present on day "x"
4 **latent** Number of latent sites present on day "x"
5 **infectious** Number of infectious sites present on day "x"
6 **removed** Number of removed sites present on day "x"
7 **senesced** Number of senesced sites present on day "x"
8 **ratinf** Rate of infection
9 **rtransfer** Rate of transfer from latent to infectious sites
10 **rgrowth** Rate of growth of healthy sites
11 **rsenesced** Rate of senescence of healthy sites
12 **rlex** Rate of lesion expansion
13 **diseased** Number of diseased (latent + infectious + removed) sites on
14 **intensity** Proportion of diseased (latent + infectious + removed) sites per total sites not including removed sites on day "x"
15 **lat** Latitude value if provided by `wth` object
16 **lon** Longitude value if provided by `wth` object
## Comparing seasons
Using the {epifitter} package (Alves and Del Ponte 2021) we can calculate the area under the disease progress curve (AUDPC) and compare disease progress in different seasons of the same year.
```{r install-epifitter}
if (!require("epifitter")) {
install.packages("epifitter")
}
library("epifitter")
```
The data are already loaded.
```{r recheck-wth}
wth
```
So we can proceed running the bacterial blight model for the dry season in the same location.
To do this, all that is necessary is to change the emergence date to January 1, 2000, which corresponds to a 120 day growing season that starts in the dry season.
```{r bb-dry-season}
bb_dry <- predict_bacterial_blight(wth, emergence = "2000-01-01")
plot(x = bb_dry$dates, y = bb_dry$intensity, type = "l")
```
### Calculating AUDPC
```{r audpc}
# wet season
AUDPC(
time = bb_wet$simday,
y = bb_wet$intensity,
y_proportion = FALSE,
type = "absolute"
)
# dry season
AUDPC(
time = bb_dry$simday,
y = bb_dry$intensity,
y_proportion = FALSE,
type = "absolute"
)
```
The AUDPC of the wet season is greater than that of the dry season.
So, this meets the expectations that the wet season AUDPC is higher than the dry season, which for bacterial blight makes sense as the pathogen needs leaf moisture and wind to move, infect and reproduce.
### Exercise 5
1. Explore different emergence dates to find a date with a higher or lower AUDPC than the current wet and dry season values.
There are two years of weather data in the `chirps` data that I've provided starting from 2000-01-01 to 2001-12-31, which allows for several establishment dates to be compared.
2. Explore the AUDPC for the other diseases in {epicrop}.
Which ones are worse in the wet season and which are worse in the dry season, do some not show as much difference between seasons?