-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathRabbit_Multistate_DeadRecovery.R
343 lines (276 loc) · 11.3 KB
/
Rabbit_Multistate_DeadRecovery.R
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
#####################################################################################
## Multi-state, dead-recovery model for 18 years of data from Turretfield rabbit pop.
## Data collected by Biosecurity, South Australia, Department of Primary Industries and Regions
##
## Code accompanies paper in review.
##
## Barnett,LK, Prowse,TAA, Peacock, DE, Mutze, GJ, Sinclair, RG, Kovaliski,J,
## Cooke, BD, Bradshaw CJA. Previous exposure to myxomatosis reduces survival
## of European rabbits during outbreaks of rabbit haemorrhagic disease for
## biological control
##
##
## Louise Barnett
## November 2017
##
## Before running this script you will need to install
## Program MARK from
## http://www.phidot.org/software/mark/downloads/
##
## Mac and Linux users might find this post helpful:
## http://www.phidot.org/forum/viewtopic.php?f=21&t=3233&p=10967&hilit=install+RMark#p10967
##
##
##
##
## testing how previous disease exposure/immunity state and recurring outbreaks
## affect rabbit survival (S) and immunity state transitions (Psi)
##
## immunity state / previous exposure categories:
## N - Immunity to neither virus
## M - Immunity to myxoma virus only
## R - Immunity to rabbit haemorrhagic disease virus (RHDV) only
## B - Immunity to both viruses
##
##
## Age groups:
## Kittens <= 600 g (may have residual maternal immunity to RHD)
## Adults > 600 g (unlikely to have residual maternal immunity)
##
#####################################################################################
######################################################################
## Run Multistate Live-recapture dead-recovery Model on rabbit data ##
######################################################################
# clear the workspace
rm(list=ls())
# load required libraries
require(RMark)
# set the working directory to folder containing capture history
# and trip covariate files
# e.g.
# setwd("~/Dropbox/RabbitDisease")
######################################################################
#################### Import and Format Data #######################
######################################################################
# Import Trip Covariates
TripCovs<- read.csv(file="TripCovariates.csv")
# Import Capture History
CH<-read.csv(file="CaptureHist.csv")
# Subset to three columns
CH<- subset(CH, select=c("ch", "Rabbit", "InitAge"))
CH$ch<- as.character(CH$ch)
# set the number of trips
TripNo=107
######################################################################
################## Make design data for the model ####################
######################################################################
# Process data
mstrata.processed=process.data(CH, model="MSLiveDead",
strata.labels = c("N","M","R","B"),# labels for strata
groups=c("InitAge"),# set initial ages
age.var=1,
initial.ages =c(0,1),
time.intervals=TripCovs$Ints,# set intervals (prop. of mean)
nocc=TripNo) # number of capture occasions
# Make design data
ddl= make.design.data(mstrata.processed,
parameters=list(Psi=list(age.bins=c(0,1,(TripNo+3)), # bin ages
subtract.stratum=c("N","M","R","B")),
S=list(age.bins=c(0,1,(TripNo+2)))),
right=FALSE)
# Ignore warning
# Check age is working
head(ddl$S$age)
head(ddl$Psi$age)
# check the subtract stratum are correct
table(ddl$Psi[,c("stratum","tostratum")])
# should have 0s on the diagonal
# for transitions estimated by subtraction
################################################################
############# Tidy up age covariates for analysis ##############
################################################################
## Set the levels of age (in the default age column)
levels(ddl$Psi$age)=c("a1kitten", "c1adult")
levels(ddl$S$age)=c("a1kitten","c1adult")
head(ddl$S$age)
## Now make a new age categories
# S
ddl$S$AgeCat<- 0
ddl$S$AgeCat[ddl$S$age=="a1kitten"]<- "Kitten"
ddl$S$AgeCat[ddl$S$age=="c1adult"]<- "Adult"
# Psi
ddl$Psi$AgeCat<- 0
ddl$Psi$AgeCat[ddl$Psi$age=="a1kitten"]<- "Kitten"
ddl$Psi$AgeCat[ddl$Psi$age=="c1adult"]<- "Adult"
# Make numeric columns/dummy variables for Kitten and Adult
# S
ddl$S$Adult<- ifelse(ddl$S$AgeCat=="Adult",1,0)
ddl$S$Kitten<- ifelse(ddl$S$AgeCat=="Kitten",1,0)
# Psi
ddl$Psi$Kitten<- ifelse(ddl$Psi$AgeCat=="Kitten",1,0)
ddl$Psi$Adult<- ifelse(ddl$Psi$AgeCat=="Adult",1,0)
################################################################
##### Add Trip Covariates to design data for p, S and Psi ######
################################################################
## Add standardised Trap Effort to capture probability (p) design data
dfp=data.frame(time=unique(ddl$p$time),TrapEffort=TripCovs$SEffort[2:(TripNo)])
ddl$p=merge_design.covariates(ddl$p,dfp)
head(ddl$p)
## Add Myxo to S
dfp=data.frame(time=unique(ddl$S$time),MV=TripCovs$MV[1:(TripNo)]) # make
ddl$S=merge_design.covariates(ddl$S,dfp)
## Add Myxo to Psi
dfp=data.frame(time=unique(ddl$Psi$time),MV=TripCovs$MV[1:(TripNo)])
ddl$Psi=merge_design.covariates(ddl$Psi,dfp)
## Add RHD to S
dfp=data.frame(time=unique(ddl$S$time),RHDV=TripCovs$RHDV[1:(TripNo)]) # make
ddl$S=merge_design.covariates(ddl$S,dfp)
head(ddl$S) # check
## Add RHD to Psi
dfp=data.frame(time=unique(ddl$Psi$time),RHDV=TripCovs$RHDV[1:(TripNo)])
ddl$Psi=merge_design.covariates(ddl$Psi,dfp)
head(ddl$Psi)
################################################################
##### fix values for transitions that are not possible #####
################################################################
# Rabbits do not lose immunity
# so set transitions that can't occur to zero
ddl$Psi$fix[ddl$Psi$stratum=="M"&ddl$Psi$tostratum=="N"] <- 0
ddl$Psi$fix[ddl$Psi$stratum=="R"&ddl$Psi$tostratum=="N"] <- 0
ddl$Psi$fix[ddl$Psi$stratum=="B"&ddl$Psi$tostratum=="N"] <- 0
ddl$Psi$fix[ddl$Psi$stratum=="M"&ddl$Psi$tostratum=="R"] <- 0
ddl$Psi$fix[ddl$Psi$stratum=="R"&ddl$Psi$tostratum=="M"] <- 0
ddl$Psi$fix[ddl$Psi$stratum=="B"&ddl$Psi$tostratum=="M"] <- 0
ddl$Psi$fix[ddl$Psi$stratum=="B"&ddl$Psi$tostratum=="R"] <- 0
################################################################
### Make dummy variables to set equal survival effects ####
################################################################
# set equal survival effect for N & M during RHDV outbreaks
ddl$S$NMduringRHD<-0
ddl$S$NMduringRHD[ddl$S$stratum=="M"&ddl$S$RHDV==1]<-1
ddl$S$NMduringRHD[ddl$S$stratum=="N"&ddl$S$RHDV==1]<-1
# set equal survival effect for N & R during MV outbreaks
ddl$S$NRduringMV<-0
ddl$S$NRduringMV[ddl$S$stratum=="R"&ddl$S$MV==1]<-1
ddl$S$NRduringMV[ddl$S$stratum=="N"&ddl$S$MV==1]<-1
# Check this
head(ddl$S)
###############################################################
######################## Run Initial Model ####################
###############################################################
# Need to run a simple version of the model first
# to set intial estimates (helps model converge)
run.init=function(){
# Process data
mstrata.processed=mstrata.processed
# Create default design data
mstrata.ddl=ddl
# Create formula for capture probability
p.stratum=list(formula=~TrapEffort)
# additive effect of prev. exposure at all times
S.stratum.d=list(formula=~-1+
Kitten+
Kitten:MV+
Kitten:RHDV+
stratum:Adult+
B:Adult:MV+
R:Adult:MV+
NMduringRHD:Adult:RHDV+
B:Adult:MV+
M:Adult:MV+
NRduringMV:Adult:MV)
# Formula for Psi
Psi.stratum=list(formula=~-1+
Kitten+
Kitten:MV+
Kitten:RHDV+
stratum:tostratum:Adult +
stratum:tostratum:Adult:MV +
stratum:tostratum:Adult:RHDV,
link="logit")
model.list=create.model.list("MSLiveDead")
mstrata.results=mark.wrapper(model.list,data=mstrata.processed,
ddl=ddl)
return(mstrata.results) }
############################################################
init<- run.init()
################################################################
############ Set up a function to run the model ###############
###############################################################
########
run.mstrata=function(){
# Process data
mstrata.processed=mstrata.processed
# Create default design data
mstrata.ddl=ddl
# Create formula
p.stratum=list(formula=~TrapEffort)
# All synergistic
S.stratum.a=list(formula=~-1+
Kitten+
Kitten:MV+
Kitten:RHDV+
stratum:Adult+
stratum:Adult:MV+
stratum:Adult:RHDV)
# only synergistic for Myxo positive in RHDV outbreaks
# additive for during Myxo outbreaks
S.stratum.b=list(formula=~-1+
Kitten+
Kitten:MV+
Kitten:RHDV+
stratum:Adult+
stratum:Adult:RHDV+
B:Adult:MV+
M:Adult:MV+
NRduringMV:Adult:MV)
# # only synergistic for RHDV positive during myxo outbreaks
# # additive during RHD outbreaks
S.stratum.c=list(formula=~-1+
Kitten+
Kitten:MV+
Kitten:RHDV+
stratum:Adult+
stratum:Adult:MV+
B:Adult:RHDV+
R:Adult:RHDV+
NMduringRHD:Adult:RHDV)
# additive effect of prev. exposure at all times
S.stratum.d=list(formula=~-1+
Kitten+
Kitten:MV+
Kitten:RHDV+
stratum:Adult+
B:Adult:RHDV+
R:Adult:RHDV+
NMduringRHD:Adult:RHDV+
B:Adult:MV+
M:Adult:MV+
NRduringMV:Adult:MV)
# Formula for Psi
Psi.stratum=list(formula=~-1+
Kitten+
Kitten:MV+
Kitten:RHDV+
stratum:tostratum:Adult +
stratum:tostratum:Adult:MV +
stratum:tostratum:Adult:RHDV,
link="logit")
model.list=create.model.list("MSLiveDead")
mstrata.results=mark.wrapper(model.list,data=mstrata.processed,
ddl=ddl, initial=init[[1]])
return(mstrata.results) }
results<- run.mstrata()
# cleanup RMark Output files
cleanup(ask=FALSE)
################################################################
###################### View output #############################
################################################################
# View model table
results
# Select model to use for figures etc.
Best<- results[[2]]
#### Save output ####
save(file="RabbitDiseaseModels.RData", results, ddl, mstrata.processed)
## Now run Script 2: 'Rabbit_Multistate_Plotting_Output.R' using saved output
##################################################################