-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathProjection_model_leastcost_territory_selection.R
439 lines (326 loc) · 21.3 KB
/
Projection_model_leastcost_territory_selection.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
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
#---
#"Projection model for baseline scenario, least cost path territory selection method"
#"Lisanne Petracca"
#"November 2023"
#---
library(tidyverse)
library(here)
library(nimble)
######------ PULLING IN FUNCTION, PACKID, OCCUPANCY, CONNECTIVITY, OTHER SPATIAL COMPONENTS ------######
#load categorical RSF territory selection information
load("data/Least_cost_territory_selection_information.RData") #LCP
#load all functions
source("functions/movement_function_singlelcp.R")
source("functions/removals_function.R")
source("functions/attraction_function.R")
#reading in spatial stuff from projection
load("data/Spatial_Information.RData")
#fixed values
proj <- 51 #100 #years of projection
nSims <- 100 #number of simulations per sample from the posterior
nSamples <- 500 #number of samples from the posterior
S <- 224 #territories
#setting up for new.guys array
newguys <- array(0, dim=c(nSamples,2,proj,S))
####### BASELINE SCENARIO ######
analysis <- "baseline"
#scenario 1 : baseline: removal rate at annual mean, immigration as estimated, no harvest, no translocation, no disease
####### THIS IS WHERE PROJECTION MODEL CODE BEGINS ######
#storage for sim loop
Nglobal_state.mean <- Nglobal_state_wmove.mean <- NAdult_state.mean <-
NAdult_EWash.mean <- NAdult_NCasc.mean <- NAdult_SCasc.mean <-
NSite_state.mean <- NSite_EWash.mean <- NSite_NCasc.mean <- NSite_SCasc.mean <- array(NA, dim = c(nSamples, proj, nSims))
Lambda.mean <- array(NA, dim = c(nSamples, proj-1, nSims))
BP_Presence <- Pack_Size <- Ntot.site <- Newguys.mean <- Two_Adult <- array(NA, dim = c(nSamples, proj, S, nSims))
lambda.mean <- lambda.upper <- lambda.lower <- p.quasiext <- p.recovery <- numeric(nSims)
for(sim in 1:nSims){
#reading in arrays needed for projection (incl. first year data)
load("data/Projection_Inputs.RData")
set.seed(37585+sim)
#the i loop is inherent here
#let's start where t==1
for (t in 1:(proj*2-1)) {
phiA.proj[,1,t] <- plogis(rnorm(nSamples, mean=int.surv1[,1], sd = sigma.period))
phiA.proj[,2,t] <- plogis(rnorm(nSamples, mean=int.surv1[,1], sd = sigma.period))
phiA.proj[,3,t] <- plogis(rnorm(nSamples, mean=int.surv1[,2], sd = sigma.period))
phiB.proj[,1,t] <- 0
phiB.proj[,2,t] <- plogis(rnorm(nSamples, mean=int.surv2[,1], sd = sigma.period))
phiB.proj[,3,t] <- plogis(rnorm(nSamples, mean=int.surv2[,2], sd = sigma.period))
}
#####---- STARTING MODEL WITH SECOND PERIOD OF YEAR 1 -----#####
#this is counting Ntot from the first time period [from other IPM run]
Ntot.proj[,1,] <-
N.proj[,1,1,1,] + N.proj[,2,1,1,] + N.proj[,3,1,1,]
#moving onto second period of first year now
for (s in 1:S){ #running through all 224 territories
#12-mo-olds
#we can actually have 12-mo movers now (edited June 2022 to add epsA)
N.stayers.proj[,1,2,1,s] <- rbinom(nSamples, N.proj[,1,1,1,s], phiA.proj[,1,1]*(1-epsA[,1]))
#new line June 2022
N.movers.proj[,1,2,1,s] <- rbinom(nSamples, N.proj[,1,1,1,s], phiA.proj[,1,1]*epsA[,1]*alpha)
#24-mo-olds
N.stayers.proj[,2,2,1,s] <- rbinom(nSamples, N.proj[,2,1,1,s], phiA.proj[,2,1]*(1-epsA[,1]))
#new and old movers - intermediate class - last period's 18-mo old residents survive and start moving but stay in state AND last period's 18-mo old movers continue moving but stay in state
#all are 1 here because from same year
N.movers.newmove.proj[,2,2,1,s] <- rbinom(nSamples, N.proj[,2,1,1,s], phiA.proj[,2,1]*epsA[,1]*alpha)
N.movers.oldmove.proj[,2,2,1,s] <- rbinom(nSamples, N.movers.proj[,2,1,1,s], phiB.proj[,2,1]*epsB[,1]*alpha)
N.movers.proj[,2,2,1,s] <- N.movers.newmove.proj[,2,2,1,s] + N.movers.oldmove.proj[,2,2,1,s]
#here are settlers - intermediate class - last periods 18-mo old new movers survive, stop moving, and settle at s
N.settlers.proj[,2,2,1,s] <- rbinom(nSamples, N.movers.proj[,2,1,1,s], phiB.proj[,2,1]*(1-epsB[,1]))
#36+-mo-olds
N.stayers.proj[,3,2,1,s] <- rbinom(nSamples, N.proj[,3,1,1,s], phiA.proj[,3,1]*(1-epsA[,2]))
#new movers - intermediate class - last period's 30-mo+ old residents survive and start moving but stay in state AND last period's 30-mo+ old movers continue moving but stay in state
#all 1 here bc same year
N.movers.newmove.proj[,3,2,1,s] <- rbinom(nSamples, N.proj[,3,1,1,s], phiA.proj[,3,1]*epsA[,2]*alpha)
N.movers.oldmove.proj[,3,2,1,s] <- rbinom(nSamples, N.movers.proj[,3,1,1,s], phiB.proj[,3,1]*epsB[,2]*alpha)
N.movers.proj[,3,2,1,s] <- N.movers.newmove.proj[,3,2,1,s] + N.movers.oldmove.proj[,3,2,1,s]
#settlers - intermediate class - last period's 30-mo+ old movers settle at s
N.settlers.proj[,3,2,1,s] <- rbinom(nSamples, N.movers.proj[,3,1,1,s], phiB.proj[,3,1]*(1-epsB[,2]))
#close s loop
}
#-----------------------------
##### FIRST MOVEMENT FUNCTION GOES HERE -----#####
n.settlers.for.fxn <- array(NA,dim = c(nSamples,2,224))
n.settlers.for.fxn <- N.settlers.proj[,c(2:3),2,1,] #just getting nSamples x 2 x site
n.res <- array(NA,dim = c(nSamples,224))
n.res <- N.stayers.proj[,1,2,1,] + N.stayers.proj[,2,2,1,] + N.stayers.proj[,3,2,1,]
new.guys <- get.move(n.settlers.for.fxn,n.res,site_check,mat_singleLCP)
for(i in 1:nSamples){ #nSamples samples
#which ids are occupied?
immig_id <- which(N.stayers.proj[i,1,2,1,] + N.stayers.proj[i,2,2,1,] + N.stayers.proj[i,3,2,1,]+
new.guys[[1]][i,] + new.guys[[2]][i,] >0)
#keeps total number of immigrants entering each year to the asymptote from the growth model (ie, limits immigration so it doesn't grow with increasing # of packs)
lambda.immig.t <- lambda.immig[i]*(assmp.immig/length(immig_id))
for(s in immig_id){
Tot.immig.proj[i,2,1,s] <- rpois(1, lambda.immig.t) #no .proj bc taken from data model
#there are no immigrant 6-11.99 mo olds
N.immig.proj[i,1,2,1,s] <- 0
#12-23.99 mo old class is the first class that can immigrate
N.immig.proj[i,2,2,1,s] <- rbinom(1, Tot.immig.proj[i,2,1,s], probImmig[1])
# group G: deterministic
N.immig.proj[i,3,2,1,s] <- Tot.immig.proj[i,2,1,s] - N.immig.proj[i,2,2,1,s]
}}
N.proj[,1,2,1,] <- N.stayers.proj[,1,2,1,]
N.proj[,2,2,1,] <- N.stayers.proj[,2,2,1,] + N.immig.proj[,2,2,1,] + new.guys[[1]] #these have rejected settlers and new guys
N.proj[,3,2,1,] <- N.stayers.proj[,3,2,1,] + N.immig.proj[,3,2,1,] + new.guys[[2]]
newguys[,2,1,] <- new.guys[[1]] + new.guys[[2]] #these are new guys only
##ATTRACTION FUNCTION HERE
n.wolves.solo.fxn <- array(NA,dim = c(nSamples,3,224))
n.wolves.solo.fxn <- N.proj[,,2,1,] #just getting nSamples x 3 x site
#solo function
group.neighbors <- get.solos(n.wolves.solo.fxn, neighbor_list)
N.proj[,,2,1,] <- group.neighbors
###### MOVING AHEAD TO T==2 ######
##### POPULATION PROJECTION MODEL #####
for (t in 1:(proj-1)){ #this is the big outer t loop
for (s in 1:S){ #running through all 224 territories
#Dec 18 mo
#last period's 12-mo olds survive and don't start moving
N.stayers.proj[,2,1,t+1,s] <- rbinom(nSamples, N.proj[,1,2,t,s], phiA.proj[,2,2*t]*(1-epsA[,1]))
#last periods 12-mo olds survive, initiate movement, but stay in state
#new lines June 2022
N.movers.newmove.proj[,2,1,t+1,s] <- rbinom(nSamples, N.proj[,1,2,t,s], phiA.proj[,2,2*t]*epsA[,1]*alpha )
N.movers.oldmove.proj[,2,1,t+1,s] <- rbinom(nSamples, N.movers.proj[,1,2,t,s], phiB.proj[,2,2*t]*epsB[,1]*alpha )
N.movers.proj[,2,1,t+1,s] <- N.movers.newmove.proj[,2,1,t+1,s] + N.movers.oldmove.proj[,2,1,t+1,s]
#can also have settlers at 18 mo
N.settlers.proj[,2,1,t+1,s] <- rbinom(nSamples, N.movers.proj[,1,2,t,s], phiB.proj[,2,2*t]*(1-epsB[,1]))
#-----------------------------
#Dec 30 mo+ (i.e., 30 mo, 42 mo, 54 mo, 66 mo, 78 mo, 90 mo...)
#key here is that this group is 30 PLUS months; has 30, 42, 54 mos
#stayers - intermediate class - last period's 24-mo old residents and 36-mo+ old residents survive and don't initiate movement
N.stayers.proj[,3,1,t+1,s] <- rbinom(nSamples, N.proj[,2,2,t,s] + N.proj[,3,2,t,s],
phiA.proj[,3,2*t]*(1-epsA[,2]))
#new and old movers - intermediate class - last period's 24-mo old and 36-mo+ old residents survive and start moving but stay in state AND last period's 24-mo old and 36-mo+ old movers continue moving but stay in state
N.movers.newmove.proj[,3,1,t+1,s] <- rbinom(nSamples, N.proj[,2,2,t,s] + N.proj[,3,2,t,s], phiA.proj[,3,2*t]*epsA[,2]*alpha )
N.movers.oldmove.proj[,3,1,t+1,s] <- rbinom(nSamples, N.movers.proj[,2,2,t,s]+ N.movers.proj[,3,2,t,s], phiB.proj[,3,2*t]*epsB[,2]*alpha )
N.movers.proj[,3,1,t+1,s] <- N.movers.newmove.proj[,3,1,t+1,s] + N.movers.oldmove.proj[,3,1,t+1,s]
#settlers - intermediate class - last periods 24-mo old movers and 36-mo+ old movers settle at s
N.settlers.proj[,3,1,t+1,s] <- rbinom(nSamples, N.movers.proj[,2,2,t,s]+N.movers.proj[,3,2,t,s], phiB.proj[,3,2*t]*(1-epsB[,2]))
} #close s loop
#Dec 6 mo olds
for(i in 1:nSamples){
for (s in 1:S){ #start s loop again
if (N.stayers.proj[i,3,1,t+1,s] >= 2){ #ok to be N.proj bc no breeding w solo indivs anyway
lambda.pups.proj[i,t+1,s] <- rcat(1,probs.pup[i,])-1 #need to subtract 1 to make it btw 0 and 6 pups
}
else
{lambda.pups.proj[i,t+1,s] <- 0}
}}
for (s in 1:S){ #assign 6-mo olds in Dec of that year
N.proj[,1,1,t+1,s] <- lambda.pups.proj[,t+1,s]
}
##### SECOND MOVEMENT FUNCTION GOES HERE -----#####
n.settlers.for.fxn <- array(NA,dim = c(nSamples,2,224))
n.settlers.for.fxn <- N.settlers.proj[,c(2:3),1,t+1,] #just getting nSamples x 2 x site
n.res <- array(NA,dim = c(nSamples,224))
n.res <- N.proj[,1,1,t+1,] + N.stayers.proj[,2,1,t+1,] + N.stayers.proj[,3,1,t+1,]
#call function
new.guys <- get.move(n.settlers.for.fxn,n.res,site_check,mat_singleLCP)
##### WE CAN ADD IMMIGRANTS HERE FOR DECEMBER
for(i in 1:nSamples){
#which ids are occupied?
immig_id <- which(N.proj[i,1,1,t+1,] + N.stayers.proj[i,2,1,t+1,] + N.stayers.proj[i,3,1,t+1,] +
new.guys[[1]][i,] + new.guys[[2]][i,] >0)
#keeps total number of immigrants entering each year to the asymptote from the growth model (ie, limits immigration so it doesn't grow with increasing # of packs)
lambda.immig.t <- lambda.immig[i]*(assmp.immig/length(immig_id))
for(s in immig_id){
Tot.immig.proj[i,1,t+1,s] <- rpois(1, lambda.immig.t) #no .proj bc taken from data model
#there are no immigrant 6-11.99 mo olds
N.immig.proj[i,1,1,t+1,s] <- 0
#12-23.99 mo old class is the first class that can immigrate
N.immig.proj[i,2,1,t+1,s] <- rbinom(1, Tot.immig.proj[i,1,t+1,s], probImmig[1])
# group G: deterministic
N.immig.proj[i,3,1,t+1,s] <- Tot.immig.proj[i,1,t+1,s] - N.immig.proj[i,2,1,t+1,s]
}}
N.proj[,2,1,t+1,] <- N.stayers.proj[,2,1,t+1,] + N.immig.proj[,2,1,t+1,] + new.guys[[1]] #these have rejected settlers and new guys
N.proj[,3,1,t+1,] <- N.stayers.proj[,3,1,t+1,] + N.immig.proj[,3,1,t+1,] + new.guys[[2]]
newguys[,1,t+1,] <- new.guys[[1]] + new.guys[[2]] #new guys only
##### REMOVAL FUNCTION GOES HERE -----#####
#THIS IS WHERE REMOVALS HAPPEN; HAPPEN ANNUALLY IN TIME PERIOD 1 (DECEMBER)
n.wolves.all.fxn <- array(NA,dim = c(nSamples,3,224))
n.wolves.EWash.fxn <- array(NA,dim = c(nSamples,3,length(EWash)))
#N.proj numbers should be going into the removal function because only sites with 2+ adults can get removed anyway
n.wolves.all.fxn <- N.proj[,,1,t+1,] #just getting nSamples x 3 x site
n.wolves.EWash.fxn <- N.proj[,,1,t+1,EWash] #just getting nSamples x 3 x site
#call function
n.postremove.EWash <- get.removals(n.wolves.all.fxn, n.wolves.EWash.fxn, removal_rate)
N.proj[,,1,t+1,EWash] <- n.postremove.EWash
##ATTRACTION FUNCTION HERE
n.wolves.solo.fxn <- array(NA,dim = c(nSamples,3,224))
n.wolves.solo.fxn <- N.proj[,,1,t+1,] #just getting nSamples x 3 x site
#solo function
group.neighbors <- get.solos(n.wolves.solo.fxn, neighbor_list)
N.proj[,,1,t+1,] <- group.neighbors
Ntot.proj[,t+1,] <-
N.proj[,1,1,t+1,] + N.proj[,2,1,t+1,] + N.proj[,3,1,t+1,]
for (s in 1:S){ #start s loop again
#-----------------------------
#Jun 12 mo
#last period's 6-mo olds survive (added epsA Jun 2022)
N.stayers.proj[,1,2,t+1,s] <- rbinom(nSamples, N.proj[,1,1,t+1,s], phiA.proj[,1,(2*t+1)]*(1-epsA[,1]))
#they can also move now (new Jun 2022)
N.movers.proj[,1,2,t+1,s] <- rbinom(nSamples, N.proj[,1,1,t+1,s], phiA.proj[,1,(2*t+1)]*epsA[,1]*alpha)
#-----------------------------
#Jun 24 mo
#stayers - intermediate class - last period's 18-mo old residents survive and don't initiate movement
N.stayers.proj[,2,2,t+1,s] <- rbinom(nSamples, N.proj[,2,1,t+1,s], phiA.proj[,2,(2*t+1)]*(1-epsA[,1]))
#new and old movers - intermediate class - last period's 18-mo old residents survive and start moving but stay in state AND last period's 18-mo old movers continue moving but stay in state
#all are t+1 here because from same year
N.movers.newmove.proj[,2,2,t+1,s] <- rbinom(nSamples, N.proj[,2,1,t+1,s], phiA.proj[,2,(2*t+1)]*epsA[,1]*alpha)
N.movers.oldmove.proj[,2,2,t+1,s] <- rbinom(nSamples, N.movers.proj[,2,1,t+1,s], phiB.proj[,2,(2*t+1)]*epsB[,1]*alpha)
N.movers.proj[,2,2,t+1,s] <- N.movers.newmove.proj[,2,2,t+1,s] + N.movers.oldmove.proj[,2,2,t+1,s]
#here are our FIRST settlers - intermediate class - last periods 18-mo old new movers survive, stop moving, and settle at s
N.settlers.proj[,2,2,t+1,s] <- rbinom(nSamples, N.movers.proj[,2,1,t+1,s], phiB.proj[,2,(2*t+1)]*(1-epsB[,1]))
#-----------------------------
#Jun 36 mo+ (i.e., 36 mo, 48 mo, 60 mo, 72 mo, 84 mo, 96 mo...)
#stayers - intermediate class - last period's 30-mo+ old residents survive and don't initiate movement
N.stayers.proj[,3,2,t+1,s] <- rbinom(nSamples, N.proj[,3,1,t+1,s], phiA.proj[,3,(2*t+1)]*(1-epsA[,2]))
#new movers - intermediate class - last period's 30-mo+ old residents survive and start moving but stay in state AND last period's 30-mo+ old movers continue moving but stay in state
#all t+1 here bc same year
N.movers.newmove.proj[,3,2,t+1,s] <- rbinom(nSamples, N.proj[,3,1,t+1,s], phiA.proj[,3,(2*t+1)]*epsA[,2]*alpha)
N.movers.oldmove.proj[,3,2,t+1,s] <- rbinom(nSamples, N.movers.proj[,3,1,t+1,s], phiB.proj[,3,(2*t+1)]*epsB[,2]*alpha)
N.movers.proj[,3,2,t+1,s] <- N.movers.newmove.proj[,3,2,t+1,s] + N.movers.oldmove.proj[,3,2,t+1,s]
#settlers - intermediate class - last period's 30-mo+ old movers settle at s
N.settlers.proj[,3,2,t+1,s] <- rbinom(nSamples, N.movers.proj[,3,1,t+1,s], phiB.proj[,3,(2*t+1)]*(1-epsB[,2]))
#-----------------------------
} #close s loop
##### THIRD MOVEMENT FUNCTION GOES HERE -----#####
dim(N.settlers.proj)
n.settlers.for.fxn <- array(NA,dim = c(nSamples,2,224))
n.settlers.for.fxn <- N.settlers.proj[,c(2:3),2,t+1,] #just getting nSamples x 2 x site
n.res <- array(NA,dim = c(nSamples,224))
n.res <- N.stayers.proj[,1,2,t+1,] + N.stayers.proj[,2,2,t+1,] + N.stayers.proj[,3,2,t+1,]
new.guys <- get.move(n.settlers.for.fxn,n.res,site_check,mat_singleLCP)
##### WE CAN ADD IMMIGRANTS HERE FOR JUNE
for(i in 1:nSamples){
#which ids are occupied?
immig_id <- which(N.stayers.proj[i,1,2,t+1,] + N.stayers.proj[i,2,2,t+1,] + N.stayers.proj[i,3,2,t+1,]+
new.guys[[1]][i,] + new.guys[[2]][i,] >0)
#keeps total number of immigrants entering each year to the asymptote from the growth model (ie, limits immigration so it doesn't grow with increasing # of packs)
lambda.immig.t <- lambda.immig[i]*(assmp.immig/length(immig_id))
for(s in immig_id){
Tot.immig.proj[i,2,t+1,s] <- rpois(1, lambda.immig.t) #no .proj bc taken from data model
#there are no immigrant 6-11.99 mo olds
N.immig.proj[i,1,2,t+1,s] <- 0
#12-23.99 mo old class is the first class that can immigrate
N.immig.proj[i,2,2,t+1,s] <- rbinom(1, Tot.immig.proj[i,2,t+1,s], probImmig[1])
# group G: deterministic
N.immig.proj[i,3,2,t+1,s] <- Tot.immig.proj[i,2,t+1,s] - N.immig.proj[i,2,2,t+1,s]
}}
N.proj[,1,2,t+1,] <- N.stayers.proj[,1,2,t+1,]
N.proj[,2,2,t+1,] <- N.stayers.proj[,2,2,t+1,] + N.immig.proj[,2,2,t+1,] + new.guys[[1]] #these have new guys and rejected settlers
N.proj[,3,2,t+1,] <- N.stayers.proj[,3,2,t+1,] + N.immig.proj[,3,2,t+1,] + new.guys[[2]]
newguys[,2,t+1,] <- new.guys[[1]] + new.guys[[2]] #new guys only
#we don't do Ntot here bc we only do that for first period
##ATTRACTION FUNCTION HERE
n.wolves.solo.fxn <- array(NA,dim = c(nSamples,3,224))
n.wolves.solo.fxn <- N.proj[,,2,t+1,] #just getting nSamples x 3 x site
#solo function
group.neighbors <- get.solos(n.wolves.solo.fxn, neighbor_list)
N.proj[,,2,t+1,] <- group.neighbors
} #close big t loop
#Ntot
for(i in 1:nSamples){
for(t in 1:proj){
#all of these are counting packs where n > 2 per site
Nglobal_state.proj[i,t] <- sum(Ntot.proj[i,t,])
Nglobal_state_wmove.proj[i,t] <- sum(Ntot.proj[i,t,]) + sum(N.movers.proj[i,,1,t,])
NAdult_state.proj[i,t] <- sum(N.proj[i,2,1,t,]) + sum(N.proj[i,3,1,t,])
NAdult_EWash.proj[i,t] <- sum(N.proj[i,2,1,t,EWash]) + sum(N.proj[i,3,1,t,EWash])
NAdult_NCasc.proj[i,t] <- sum(N.proj[i,2,1,t,NorthCasc]) + sum(N.proj[i,3,1,t,NorthCasc])
NAdult_SCasc.proj[i,t] <- sum(N.proj[i,2,1,t,SouthCasc]) + sum(N.proj[i,3,1,t,SouthCasc])
NSite_state.proj[i,t] <- length(which((N.proj[i,2,1,t,] + N.proj[i,3,1,t,]>=2) & N.proj[i,1,1,t,]>=2))
NSite_EWash.proj[i,t] <- length(which((N.proj[i,2,1,t,EWash] + N.proj[i,3,1,t,EWash]>=2) & N.proj[i,1,1,t,EWash]>=2))
NSite_NCasc.proj[i,t] <- length(which((N.proj[i,2,1,t,NorthCasc] + N.proj[i,3,1,t,NorthCasc]>=2) & N.proj[i,1,1,t,NorthCasc]>=2))
NSite_SCasc.proj[i,t] <- length(which((N.proj[i,2,1,t,SouthCasc] + N.proj[i,3,1,t,SouthCasc]>=2) & N.proj[i,1,1,t,SouthCasc]>=2))
for(s in 1:S){
N_newguys.proj[i,t,s] <- sum(newguys[i,,t,s])
BP_presence.proj[i,t,s] <- ifelse(N.proj[i,2,1,t,s] + N.proj[i,3,1,t,s]>=2 & N.proj[i,1,1,t,s]>=2,1,0)
Two_Adult.proj[i,t,s] <- ifelse(N.proj[i,2,1,t,s] + N.proj[i,3,1,t,s]>=2,1,0)
Pack_Size.proj[i,t,s] <- N.proj[i,1,1,t,s] + N.proj[i,2,1,t,s] + N.proj[i,3,1,t,s]
}}}
#growth rate over the whole study period
lambda.proj <- matrix(NA, nrow=nSamples, ncol=proj-1)
#need to replace Inf with NA in Nglobal_state.proj
for (t in 2:proj) {
# mean and quantiles per year across mcmc samples; leave na.rm for zero size pops in mcmc samples
lambda.proj[,t-1] <- Nglobal_state.proj[,t]/Nglobal_state.proj[,t-1]
# we need the Infs to become NAs
lambda.proj[,t-1][is.infinite(lambda.proj[,t-1])] <- NA
lambda.proj[,t-1][is.nan(lambda.proj[,t-1])] <- NA
} #closes t on lambda
#### derive and store values for each simulation
Lambda.mean[1:nSamples,,sim] <- as.matrix(lambda.proj)
Nglobal_state.mean[1:nSamples,,sim] <- as.matrix(Nglobal_state.proj)
Nglobal_state_wmove.mean[1:nSamples,,sim] <- as.matrix(Nglobal_state_wmove.proj)
NAdult_state.mean[1:nSamples,,sim] <- as.matrix(NAdult_state.proj)
NAdult_EWash.mean[1:nSamples,,sim] <- as.matrix(NAdult_EWash.proj)
NAdult_NCasc.mean[1:nSamples,,sim] <- as.matrix(NAdult_NCasc.proj)
NAdult_SCasc.mean[1:nSamples,,sim] <- as.matrix(NAdult_SCasc.proj)
NSite_state.mean[1:nSamples,,sim] <- as.matrix(NSite_state.proj)
NSite_EWash.mean[1:nSamples,,sim] <- as.matrix(NSite_EWash.proj)
NSite_NCasc.mean[1:nSamples,,sim] <- as.matrix(NSite_NCasc.proj)
NSite_SCasc.mean[1:nSamples,,sim] <- as.matrix(NSite_SCasc.proj)
Newguys.mean[,,,sim] <- as.array(N_newguys.proj)
Ntot.site[,,,sim] <- as.array(Ntot.proj)
BP_Presence[,,,sim] <- as.array(BP_presence.proj)
Two_Adult[,,,sim] <- as.array(Two_Adult.proj)
Pack_Size[,,,sim] <- as.array(Pack_Size.proj)
} #close sim
dim(BP_Presence)
#this will give max pack size across x samples
Pack_Size_max <- apply(Pack_Size,c(1),max)
#this will get probability of having BP by site and year
BP_Presence_summary <- apply(BP_Presence,c(2,3),mean)
Two_Adult_summary <- apply(Two_Adult,c(2,3),mean)
#this will get mean and median wolves by site and year, and mean new guys
Ntot.site_mean <- apply(Ntot.site,c(2,3),mean)
Ntot.site_median <- apply(Ntot.site,c(2,3),mean)
Newguys.mean <- apply(Newguys.mean,c(2,3),mean)
save(Lambda.mean,
Ntot.site_mean, Ntot.site_median, Newguys.mean,
BP_Presence_summary, BP_Presence, Pack_Size_max,Two_Adult_summary,Two_Adult,
Nglobal_state.mean, Nglobal_state_wmove.mean,
NAdult_state.mean,
NAdult_EWash.mean, NAdult_NCasc.mean, NAdult_SCasc.mean, Newguys.mean,
NSite_state.mean, NSite_EWash.mean, NSite_NCasc.mean, NSite_SCasc.mean, file="baseline_leastcost.RData")