-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path1_Model construction.Rmd
351 lines (261 loc) · 18.2 KB
/
1_Model construction.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
---
title: |-
Document 1- R codes for the construction of the stage-structured matrix model for spotted hyenas infected with CDV
author: "S. Benhaiem, L. Marescot"
date: '`r Sys.time()`'
output:
pdf_document: default
html_document:
code_folding: hide
word_document: default
---
<style>
/* Style the linenumber div */
.linenumbers {
border: 1px solid #ccc;
border-radius: 4px;
background-color: #EBEBEB;
text-align: center;
padding: 0px 3px;
font-family: monospace;
float: left;
position: absolute;
transform:translate(-125%);
font-size: inherit !important;
}
.spoiler, .spoiler > * { transition: color 0.5s, opacity 0.5s }
.spoiler:not(:hover) { color: transparent }
.spoiler:not(:hover) > * { opacity: 0 }
/* fix weird transitions on Chrome: */
blockquote, blockquote > *:not(a) { color: black }
</style>
This code loads the input values for the 3 epidemic periods and builds the matrix model. It is structured as follows:
a) Loading input values (Multi-Event-Capture-Mark-Recapture (MECMR) parameter estimates)
b) Building the submatrices: Survival, State transition (demographic, social, infection) and Fecundity submatrices
c) Assembling the submatrices into the meta-matrix
d) Building the next generation matrix for the estimation of R0
Note that parameter and submatrix names may differ between main text and R codes.
--------------------------------------------------------------------------------------------------------
#Construction of the matrix model#
##a) Loading the input values (Multi-Event-Capture-Mark-Recapture (MECMR) parameter estimates)##
First, we load the R package popbio (Stubben & Milligan 2007) and the Maximum Likelihood Estimates (MLE) from the MECMR model fitted in E-surge, for each epidemic period. These estimates correspond to the values presented in Table 1 in the main text.
```{r, message=FALSE}
library(popbio)
source('Input_90-92(pre-epidem).R')
source('Input_93-94(epidem).R')
source('Input_95-99(post-epidem).R')
```
The Table 1 below shows the probability of surviving ($\phi$), becoming infected ($\beta$), a breeder ($\psi$), and of maintaining the current social state (r) in female spotted hyenas during each period, as estimated via the MECMR model. Demographic states were abbreviated as cubs (C), subadults (SA), breeders (B) or non-breeders (NB), social states as high status (H) or low status (L) and infection states as susceptible (S), infected (I) or recovered (R). This table is equivalent to Table 1 in the main text.
Parameter | pre-epidemic | epidemic | post-epidemic
------------- | -------------- | ------------------ | ----------
$\phi_{C.L.S}$ | 0.86 | 0.78 | 0.81
$\phi_{C.L.I}$ | 0.62 |0.48 | 0.54
$\phi_{C.H.S}$ | 0.90 | 0.84 | 0.72
$\phi_{C.H.I}$ | 0.79 | 0.67 | 0.86
$\phi_{SA.L.S}$ | 0.90 | 0.83 | 0.86
$\phi_{SA.LH.IR}$ | 0.69 |0.56 | 0.61
$\phi_{SA.H.S}$ | 0.97 | 0.95 | 0.96
$\phi_{B}$ | 0.95 | 0.91 | 0.93
$\phi_{NB}$ | 0.86 | 0.78 | 0.82
$r_{LL}$ | 0.97 | 0.97 | 0.97
$r_{HH}$ | 0.94 | 0.94 | 0.94
$\beta_{C.L}$ | 0.45 | 0.99 | 0.85
$\beta_{C.H}$ | 0.23 |0.96 | 0.67
$\beta_{SABNB.L}$ | 0.02 | 0.66 | 0.13
$\beta_{SABNB.H}$ | 0.45 | 0.90 | 0.42
$b_{SA.L}$ | 0.01 | 0.01 | 0.01
$b_{B.L}$ | 0.45 |0.45 | 0.45
$b_{NB}$ | 0.60 |0.60 | 0.60
$b_{SA.H}$ | 0.04 |0.04 | 0.04
$b_{B.H}$ | 0.60 |0.49 | 0.49
$b_{B.H}$ | 0.68 | 0.68 | 0.68
##b) Building the submatrices: Survival, State transition (demographic, social, infection) and Fecundity submatrices##
The definitions of the parameters are provided in the input files (see below section 2)a)).
# Survival
Below we construct 4 diagonal matrices representing the survival process for cubs, subadults, breeders and non breeders, respectively. For instance, elements in the di agonal of the matrix Surv.C represent the survival probabilites of cubs given their social and infection state.
$$Surv.C = \left[\begin{array}
{rrrrrr}
\phi_{C.L.S} & 0 & 0 & 0 & 0 & 0\\
0 & \phi_{C.L.I} & 0 & 0 & 0 & 0\\
0 & 0 & \phi_{C.L.R} & 0 & 0 & 0\\
0 & 0 & 0 & \phi_{C.H.S} & 0 & 0\\
0 & 0 & 0 & 0 & \phi_{C.H.I} & 0\\
0 & 0 & 0 & 0 & 0 & \phi_{C.H.R}\\
\end{array}\right]
$$
```{r}
Surv.C <- diag(c(phiCLS, phiCLI, phiCLR, phiCHS, phiCHI,phiCHR))
Surv.S <- diag(c(phiSLS, phiSLI, phiSLR, phiSHS, phiSHI, phiSHR))
Surv.B <- diag(c(phiBLS, phiBLI, phiBLR, phiBHS, phiBHI, phiBHR))
Surv.NB <- diag(c(phiNBLS, phiNBLI,phiNBLR, phiNBHS, phiNBHI, phiNBHR))
```
# State transitions
Each entry in the submatrices described below correspond to the probability of transition (or survival), or the fecundity, from a 'starting' combination of social and infection states to a 'following' combination of social and infection states. The 6 rows and 6 columns in these submatrices correspond to: Low social status - Susceptible, Low social status - Infected, Low social status - Recovered, High social status - Susceptible, High social status - Infected, High social status - Recovered.
*Demographic states*
We build six demography submatrices (i.e. recruitment and non-recruitment submatrices); three to account for the transition to the breeder state and three others to account for the transition to the non-breeder state. Both states are accessible from subadults, breeders and non-breeders.
```{r}
Rec.S <- diag(c(bSL, bSL, bSL, bSH, bSH,bSH))
Rec.B <- diag(c(bBL, bBL, bBL, bBH, bBH,bBH))
Rec.NB <- diag(c(bNBL, bNBL, bNBL, bNBH, bNBH,bNBH))
NRec.S <- diag(c((1-bSL), (1-bSL), (1-bSL),(1-bSH), (1-bSH),(1-bSH)))
NRec.B <- diag(c((1-bBL), (1-bBL), (1-bBL),(1-bBH), (1-bBH),(1-bBH)))
NRec.NB <- diag(c((1-bNBL), (1-bNBL), (1-bNBL),(1-bNBH), (1-bNBH),(1-bNBH)))
```
*Social states*
We build a social submatrix to account for the fact that subadults, breeders and non-breeders can either stay within their social state or change it. Cubs had the same social state as their mother.
```{r}
R <- matrix(c(
rLL.S, 0, 0, 1-rHH.S, 0, 0,
0, rLL.I, 0, 0, 1-rHH.I, 0,
0, 0, rLL.R, 0, 0, 1-rHH.R,
1-rLL.S, 0, 0, rHH.S, 0, 0,
0, 1-rLL.I, 0, 0, rHH.I, 0,
0, 0, 1-rLL.R, 0, 0, rHH.R)
,byrow=T,nrow=6)
```
*Infection states*
We build three infection submatrices: one for newborns, one for cubs, and one for subadults, breeders and non-breeders. In the matrix of infection for newborns (I.F), we considered that mothers only produced susceptible newborns, irrespective of their own infection status. Indeed, as adults do not excrete the virus actively as young individuals, we considered the transmission of CDV to be strictly horizontal.
```{r}
I.F <- matrix(c(
1, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 1, 1,
0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0)
,byrow=T,nrow=6)
```
In the matrix of infection for cubs (I.C), the only possible infection transitions were from a susceptible to a susceptible state or from a susceptible to an infected state. As recovered cubs were absent from our original CMR dataset, we set the columns for recovered cubs to 0.
```{r}
I.C <- matrix(c(
1-betaCL, 0, 0, 0, 0, 0,
betaCL, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0,
0, 0, 0, 1-betaCH, 0, 0,
0, 0, 0, betaCH, 0, 0,
0, 0, 0, 0, 1, 0)
,byrow=T,nrow=6)
```
The structure of the matrices of infection for subadults (I.SA), breeders (I.AD) and non-breeders (I.AD) is similar as the one for cubs (except that the transition from recovered to recovered occurs (3rd column and 3rd row for low-ranking and 6th column and 6th row for high-ranking).
For instance, the possible transitions in the infection matrix of subadults (I.SA) are: from suceptible (first (for low-ranking) and fourth (for high-ranking) columns) to susceptible (first (for low-ranking) and fourth (for high-ranking) rows) with probability $1-\beta_{SA.L.S}$ or $1-\beta_{SA.H.S}$, from susceptible to infected (second and fifth rows) with probability $\beta_{SA.L.S}$ or $\beta_{SA.H.S}$, and from infected (second and fifth columns)to recovered (third and sixth rows) with probabilities equal to 1 for the transitions I->R and R->R.
$$I.SA = \left[\begin{array}
{rrrrrr}
1-\beta_{SA.L} & 0 & 0 & 0 & 0 & 0\\
\beta_{SA.L} & 0 & 0 & 0 & 0 & 0\\
0 & 1 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 1-\beta_{SA.H} & 0 & 0\\
0 & 0 & 0 & \beta_{SA.H} & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 1\\
\end{array}\right]
$$
```{r}
I.SA <- matrix(c(
1-betaSL, 0, 0, 0, 0, 0,
betaSL, 0, 0, 0, 0, 0,
0, 1, 1, 0, 0, 0,
0, 0, 0, 1-betaSH, 0, 0,
0, 0, 0, betaSH, 0, 0,
0, 0, 0, 0, 1, 1)
,byrow=T,nrow=6)
I.AD <- matrix(c(
1-betaADL, 0, 0, 0, 0, 0,
betaADL, 0, 0, 0, 0, 0,
0, 1, 1, 0, 0, 0,
0, 0, 0, 1-betaADH, 0, 0,
0, 0, 0, betaADH, 0, 0,
0, 0, 0, 0, 1, 1)
,byrow=T,nrow=6)
```
# Fecundity
We construct the fecundity matrix for breeder females which is a diagonal matrix of 1, as all females which are in a breeder state a given year produce offspring with a probability of 1. The matrix is then multiplied by the litter size and the sex ratio (as the model is female-based), which are constant across social and infection states.
```{r}
fB <- ls* sr * diag(c(fLS, fLI, fLR, fHS, fHI,fHR))
```
##c) Assembling the submatrices into the meta-matrix##
All these processes described in 1.b) are then combined via multiplications of matrices (%*% in R) to generate successive events within a one year projection interval: 1) change in demographic state, 2) change in social state, 3) change in infection state and then 4) survival (the matrices below should be read from right to left). We assumed that all births and deaths occurred simultaneously at the end of the projection interval, that is, we modelled the population with a birth pulse reproduction and a pre-breeding census.
```{r}
NetF<-Surv.C%*%I.C%*%I.F%*% fB # Net Fecundity submatrix (called "Fertility" in the manuscript)
CS <- Surv.S%*%I.SA # the survival-transition submatrix to the subadult state from the cub state
SB <- Surv.B%*%I.AD%*%R%*%Rec.S # the survival-transition submatrix to the breeder state from the subadult state
BB <- Surv.B%*%I.AD%*%R%*%Rec.B # the survival-transition submatrix to the breeder state from the breeder state
NBB <- Surv.B%*%I.AD%*%R%*%Rec.NB # the survival-transition submatrix to the breeder state from the non-breeder state
SNB <- Surv.NB%*%I.AD%*%R%*%NRec.S # the survival-transition submatrix to the non-breeder state from the subadult state
BNB <- Surv.NB%*%I.AD%*%R%*%NRec.B # the survival-transition submatrix to the breeder state from the non-breeder state
NBNB<- Surv.NB%*%I.AD%*%R%*%NRec.NB # the survival-transition submatrix to the non-breeder state from the non-breeder state
```
Finally, we combine these survival-transition submatrices within a single meta-matrix M.final
$$M.final = \left[\begin{array}
{rrrr}
Z & Z & NetF & Z\\
CS & Z & Z & Z\\
Z & SB & BB & NB\\
Z & SB & BB & NBB\\
Z & SNB & BNB & NBNB\\
\end{array}\right]
$$
```{r}
Z <- matrix(0,6,6)
M <- rbind(cbind (Z, Z, NetF, Z),
cbind (CS, Z, Z, Z),
cbind (Z, SB, BB, NBB),
cbind (Z, SNB, BNB, NBNB))
rownames(M)<- c("CLS","CLI","CLR","CHS","CHI","CHR",
"SLS","SLI","SLR","SHS","SHI","SHR",
"BLS","BLI","BLR","BHS","BHI","BHR",
"NBLS","NBLI","NBLR","NBHS","NBHI","NBHR")
colnames(M)<- c("CLS","CLI","CLR","CHS","CHI","CHR",
"SLS","SLI","SLR","SHS","SHI","SHR",
"BLS","BLI","BLR","BHS","BHI","BHR",
"NBLS","NBLI","NBLR","NBHS","NBHI","NBHR")
M.final <- M[-c(3,6), -c(3,6)] # We remove the recovered cubs as they do not exist in the data
```
##d) Building the next generation matrix for the estimation of R0 ##
```{r}
### -- Building the transition matrix Tr.final:
Tr <- rbind(cbind (Z, Z, Z, Z),
cbind (CS, Z, Z, Z),
cbind (Z, SB, BB, NBB),
cbind (Z, SNB, BNB, NBNB))
Tr.final <- Tr[-c(3,6), -c(3,6)]
### -- Building the identity matrix I:
Id <- diag(22)
### --- Building N, the fundamental matrix as N = (Id-Tr.final)-^1:
N <- solve((lambda(M.final)*Id) -Tr.final) # Here we include the population's growth rate
## ---- Building F, the fertility matrix:
F <- matrix(c(
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaCL, 0, phiCHS*betaCL, 0, phiSLS*betaCL, 0, 0, phiSHS*betaCL, 0, 0, phiBLS*betaCL, 0, 0, phiBHS*betaCL, 0, 0, phiNBLS*betaCL, 0, 0, phiNBHS*betaCL,0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaCH, 0, phiCHS*betaCH, 0, phiSLS*betaCH, 0, 0, phiSHS*betaCH, 0, 0, phiBLS*betaCH, 0, 0, phiBHS*betaCH, 0, 0, phiNBLS*betaCH, 0, 0, phiNBHS*betaCH, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaSL, 0, phiCHS*betaSL, 0, phiSLS*betaSL, 0, 0, phiSHS*betaSL, 0, 0, phiBLS*betaSL, 0, 0, phiBHS*betaSL, 0, 0, phiNBLS*betaSL, 0, 0, phiNBHS*betaSL, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaSH, 0, phiCHS*betaSH, 0, phiSLS*betaSH, 0, 0, phiSHS*betaSH, 0, 0, phiBLS*betaSH, 0, 0, phiBHS*betaSH, 0, 0, phiNBLS*betaSH, 0, 0, phiNBHS*betaSH, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaADL, 0, phiCHS*betaADL, 0, phiSLS*betaADL, 0, 0, phiSHS*betaADL, 0, 0, phiBLS*betaADL, 0, 0, phiBHS*betaADL, 0, 0, phiNBLS*betaADL, 0, 0, phiNBHS*betaADL, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaADH, 0, phiCHS*betaADH, 0, phiSLS*betaADH, 0, 0, phiSHS*betaADH, 0, 0, phiBLS*betaADH, 0, 0, phiBHS*betaADH, 0, 0, phiNBLS*betaADH, 0, 0, phiNBHS*betaADH, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaADL, 0, phiCHS*betaADL, 0, phiSLS*betaADL, 0, 0, phiSHS*betaADL, 0, 0, phiBLS*betaADL, 0, 0, phiBHS*betaADL, 0, 0, phiNBLS*betaADL, 0, 0, phiNBHS*betaADL, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, phiCLS*betaADH, 0, phiCHS*betaADH, 0, phiSLS*betaADH, 0, 0, phiSHS*betaADH, 0, 0, phiBLS*betaADH, 0, 0, phiBHS*betaADH, 0, 0, phiNBLS*betaADH, 0, 0, phiNBHS*betaADH, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
),byrow=T,nrow=22)
### --- Next generation Matrix NGM:
NGM <- F %*% N
rownames(NGM)<- c("CLS","CLI","CHS","CHI",
"SLS","SLI","SLR","SHS","SHI","SHR",
"BLS","BLI","BLR","BHS","BHI","BHR",
"NBLS","NBLI","NBLR","NBHS","NBHI","NBHR")
colnames(NGM)<- c("CLS","CLI","CHS","CHI",
"SLS","SLI","SLR","SHS","SHI","SHR",
"BLS","BLI","BLR","BHS","BHI","BHR",
"NBLS","NBLI","NBLR","NBHS","NBHI","NBHR")
```
## Reference ##
Stubben, C. & Milligan, B. Estimating and analyzing demographic models using the popbio package in R. J. Stat. Softw 22, 1-23 (2007).