-
Notifications
You must be signed in to change notification settings - Fork 1
/
how-to.Rmd
248 lines (177 loc) · 6.69 KB
/
how-to.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
---
title: "Notes on coding a Lotka-Volterra Model"
output: github_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Let's start with just figuring out how we we code up two populations if they had equal per-capita probabilities of birth and death. Hey! This is just a simplified kind of neutral theory!
```{r non-loop}
# initial population sizes
n <- c(10, 10)
# a matrix to hold populations sizes through time (rows are time points)
nmat <- matrix(0, nrow = 10, ncol = 2)
nmat[1, ] <- n # record the initial conditions
# birth or death?
# `event = 0` is death and `event = 1` is birth
event <- sample(0:1, 1)
if(event == 0) { # someone dies
# who dies?
thisOneDead <- sample(1:2, size = 1, prob = n)
# update populations
n[thisOneDead] <- n[thisOneDead] - 1
} else { # someone is born
thisOneBorn <- sample(1:2, size = 1, prob = n)
# update populations
n[thisOneBorn] <- n[thisOneBorn] + 1
}
# now `n` has changed
n
# update the matrix by hand
nmat[2, ] <- n
```
The issue here is that filling in that `nmat` matrix by hand would take a *long* time. Luckily we can ask the computer to do the work for us instead (and a lot faster). We use a for loop. Check out what a for loop does
```{r what-loop}
# initial population sizes
n <- c(10, 10)
# a matrix to hold populations sizes through time (rows are time points)
nmat <- matrix(0, nrow = 10, ncol = 2)
nmat[1, ] <- n # record the initial conditions
# here's the for loop, so called cause we use the word "for"
for(i in 2:10) {
nmat[i, 1] <- i
}
nmat
```
So the for loop iterates over the values 2 through 10 and *for* each value it assigns the value to `i` and then let's us do anything we want with `i`. Here we used `i` both to index the matrix, and to provide a new value.
We can use a for loop to run our neutral simulation.
```{r neutral}
# initial population sizes
n <- c(10, 10)
# a matrix to hold populations sizes through time (rows are time points)
nmat <- matrix(0, nrow = 10, ncol = 2)
nmat[1, ] <- n # record the initial conditions
# here's the for loop, so called cause we use the word "for"
for(i in 2:10) {
# birth or death?
# `event = 0` is death and `event = 1` is birth
event <- sample(0:1, 1)
if(event == 0) { # someone dies
# who dies?
thisOneDead <- sample(1:2, size = 1, prob = n)
# update populations
n[thisOneDead] <- n[thisOneDead] - 1
} else { # someone is born
thisOneBorn <- sample(1:2, size = 1, prob = n)
# update populations
n[thisOneBorn] <- n[thisOneBorn] + 1
}
# update the matrix
nmat[i, ] <- n
}
nmat
```
Now we can see `nmat` got totally filled in with a real simulation! We could visualize it like this
```{r first-viz}
plot(1:10, nmat[, 1], ylim = c(0, max(nmat)), col = 'blue', type = 'l',
xlab = 'Time steps', ylab = 'n')
lines(1:10, nmat[, 2], col = 'orange')
```
The nice thing about a for loop is we can easily run it a lot longer to see genuinely interesting dynamics
```{r neutral-long}
# initial population sizes
n <- c(10, 10)
# a matrix to hold populations sizes through time (rows are time points)
ntstep <- 500 # number of time steps
nmat <- matrix(0, nrow = ntstep, ncol = 2)
nmat[1, ] <- n # record the initial conditions
# here's the for loop, so called cause we use the word "for"
for(i in 2:ntstep) {
# everything could going extinct, so this is a safety precaution
# it will stop the loop if everything is dead
if(all(n <= 0)) break
# birth or death?
# `event = 0` is death and `event = 1` is birth
event <- sample(0:1, 1)
if(event == 0) { # someone dies
# who dies?
thisOneDead <- sample(1:2, size = 1, prob = n)
# update populations
n[thisOneDead] <- n[thisOneDead] - 1
} else { # someone is born
thisOneBorn <- sample(1:2, size = 1, prob = n)
# update populations
n[thisOneBorn] <- n[thisOneBorn] + 1
}
# update the matrix
nmat[i, ] <- n
}
plot(1:ntstep, nmat[, 1], ylim = c(0, max(nmat)), col = 'blue', type = 'l',
xlab = 'Time step', ylab = 'n')
lines(1:ntstep, nmat[, 2], col = 'orange')
```
But we're not trying for a neutral model, we're trying for Lotka-Voltera. What do we need to do? We need to have per capita death rates that are not all equal. Specifically for a species $i$ competing with a species $j$ we need per capita death rate $D_i$ to be
$$
D_i = \mu_{i0} + \mu_{ii}n_i + \mu_{ij}n_j
$$
We also have to have a birth rate, but this is inconsequential in terms of breaking neutrality. In fact we can just have both birth rates be the same
$$
B_i = B_j = \lambda
$$
Let's set values to these parameters
```{r full-param}
# birth
lambda <- 4
# intrinsic death rates
mu10 <- 0.2
mu20 <- 0.4
# intra-specific competition
mu11 <- 0.04
mu22 <- 0.06
# inter-specific competition
mu12 <- 0.002
mu21 <- 0.003
```
Now we can write a new for loop with Lotka-Volterra death
```{r lv-full}
# initial population sizes
n <- c(10, 10)
# a matrix to hold populations sizes through time (rows are time points)
ntstep <- 5000 # number of time steps
nmat <- matrix(0, nrow = ntstep, ncol = 2)
nmat[1, ] <- n # record the initial conditions
# here's the for loop, so called cause we use the word "for"
for(i in 2:ntstep) {
# everything could going extinct, so this is a safety precaution
# it will stop the loop if everything is dead
if(all(n <= 0)) break
# birth or death?
# death probabilities
dprob1 <- n[1] * (mu10 + mu11 * n[1] + mu12 * n[2])
dprob2 <- n[2] * (mu20 + mu22 * n[2] + mu21 * n[1])
# birth probabilities
bprob1 <- lambda * n[1]
bprob2 <- lambda * n[2]
# `event = 0` is death and `event = 1` is birth
event <- sample(0:1, size = 1, prob = c(dprob1 + dprob2, bprob1 + bprob2))
if(event == 0) { # someone dies
# who dies?
thisOneDead <- sample(1:2, size = 1, prob = c(dprob1, dprob2))
# update populations
n[thisOneDead] <- n[thisOneDead] - 1
} else { # someone is born
thisOneBorn <- sample(1:2, size = 1, prob = c(bprob1, bprob2))
# update populations
n[thisOneBorn] <- n[thisOneBorn] + 1
}
nmat[i, ] <- n
}
plot(1:ntstep, nmat[, 1], type = 'l', ylim = c(0, max(nmat)), col = 'blue',
xlab = 'Time steps', ylab = 'n')
lines(1:ntstep, nmat[, 2], col = 'orange')
```
Now talk about modern coexist insight: $\rho < f/f < 1/\rho$
and bio interpretations
Then talk about how multi-spp falls apart in terms of predicting coexist
Can students find an example of where pairwise coexist is predicted, but 3-way doesn't work?
simualte species rich community and look at SAD, look familiar?