-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathS2 Simulate data.R
188 lines (143 loc) · 6.1 KB
/
S2 Simulate data.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
#############################################################################################################################
#####
##### THIS SCRIPT USED TO SIMULATE DATA FOR ANALYSIS IN THE
##### HIERARCHICAL INTEGRATED POP. MODEL PRESENTED IN "Nilsen, EB & Strand, O. (2017). Integrating data from several sources for
##### increased insight into demographic processes: Simulation studies and proof of concept for hierarchical change in ratio models. PlosOne"
SimPop <- function(m=1, phi1t="Nei", ft="Nei", p1t="Nei", p2t="Nei", T=31, PHI3 = 0.95, mean_F=0.85, f_sd=0.1, mean_phi1 = 0.85, phi1_sd=0.1,
p1=0.6, p1_sd=0.1, p2=0.4, p2_sd=0.1,
bias1 = 1, bias2 = 1, bias3=1, error.dist="pois", p.obs=NULL){
# m: Multiplier - for the initial poulation size
# T: Number of time steps to simulate
# PHI3: Adult annual survival probability
# mean_F: Mean number of calves pr. female,
# PHI1: Juvenile summer survival probability
# PHI1: JUvenile overwinter survival probability
# p1: Observation probability for calf surveys in spring
# p2: Observation probability for structur surveys post harvest
# bias1: Bias in observation data; Calf surveys in spring. bias1>1 implies higher detection probability
# for females/yearlings. bias1<1 implies higher detection probability
# for calves.
# bias2: Bias in observation data; Structural surveys post harvest. bias2>1 implies higher detection probability for
# females compared to calves and males. bias2<1 implies lower detection probability for calves than females and males
# bias3: Bias in observation data; Structural surveys post harvest. bias3>1 implies higher detection probability for
# males compared to calves and females
# error.dist: Error distribution for count data
# p.obs: Observation probaiblity for count data
## POPULATION VECTORS
N <- matrix(ncol=T, nrow=6) ## Pre harvest pop. vector. No monitoring
X <- matrix(ncol=T, nrow=6) ## Post harvest pop. vector. Monitored by str. surveys (calves, males, females) and potentially tot. surveys
N_juv <- matrix(ncol=T, nrow=2) ## Spring pop. vector. Monitored by recruitment surveys (calves vs yearlings and females)
N_tot <- X_tot <- N_obs <- J <- SU <- C0 <- Cf <- Cm <- matrix(ncol=T, nrow=1)
## OBSERVATION PARAMETERS; PRE-HARVEST SURVEYS
pa <- matrix(ncol=T, nrow=1)
if(p1t=="Ja"){
for(i in 1:T){
pa[i] <- betaval(p1, p1_sd)
}
}
if(p1t=="Nei"){
for(i in 1:T){
pa[i] <- p1
}
}
## OBSERVATION PARAMETERS; POST-HARVEST SURVEYS
pb <- matrix(ncol=T, nrow=1)
if(p2t=="Ja"){
for(i in 1:T){
pb[i] <- betaval(p1, p1_sd)
}
}
if(p2t=="Nei"){
for(i in 1:T){
pb[i] <- p2
}
}
## DEMOGRAPHIC PARAMETERS
# fecundity;
f <- matrix(ncol=T, nrow=1)
if(ft=="Ja"){
for(i in 1:T){
f[i] <- betaval(mean_F, f_sd)
}
}
if(ft=="Nei"){
for(i in 1:T){
f[i] <- mean_F
}
}
# Juvenile summer survival
phi1 <- matrix(ncol=T, nrow=1)
if(phi1t=="Ja"){
for(i in 1:T){
phi1[i] <- betaval(mean_phi1, phi1_sd)
}
}
if(phi1t=="Nei"){
for(i in 1:T){
phi1[i] <- mean_phi1
}
}
############################################################
## HARVEST OFF-TAKE
h <- matrix(ncol=T, nrow=6)
h[1,] <- sample(size=T, seq(0.2,0.3, by=0.01), replace=T)
h[2,] <- sample(size=T, seq(0.15,0.30, by=0.01), replace=T)
h[3,] <- sample(size=T, seq(0.1,0.15, by=0.01), replace=T)
h[4,] <- sample(size=T, seq(0.2,0.3, by=0.01), replace=T)
h[5,] <- sample(size=T, seq(0.05,0.1, by=0.01), replace=T)
h[6,] <- sample(size=T, seq(0.1,0.40, by=0.01), replace=T)
############################################################
## DERIVED HARVEST NUMBERS
H <- matrix(ncol=T, nrow=6)
#############################################################
### INITIAL POPULATION SIZE AND STRUCTURE - FOR X
N[1,1] <- 120*m
N[2,1] <- 120*m
N[3,1] <- 90*m
N[4,1] <- 80*m
N[5,1] <- 150*m
N[6,1] <- 180*m
for (t in 1:(T-1)){
###########################################################
# STATE PROCESS;
# PRE-HARVEST POPULATION VECTORS IN T+1
N[3,t+1] <- rbinom(1, round(N[1,t]*(1-h[1,t])), PHI3)
N[4,t+1] <- rbinom(1, round(N[2,t]*(1-h[2,t])), PHI3)
N[5,t+1] <- rbinom(1, round((N[3,t]*(1-h[3,t]))+(N[5,t]*(1-h[5,t]))), PHI3)
N[6,t+1] <- rbinom(1, round((N[4,t]*(1-h[4,t]))+(N[6,t]*(1-h[6,t]))), PHI3)
N[1,t+1] <- rbinom(1, N[5, t+1], phi1[t]*f[t]/2)
N[2,t+1] <- rbinom(1, N[5, t+1], phi1[t]*f[t]/2)
N_juv[1, t+1] <- round(N[5,t+1]*(f[t]/2)) # Female calves in spring
N_juv[2, t+1] <- round(N[5,t+1]*(f[t]/2)) # Male calves in spring
#############################################################
# POST-HARVEST POPULATION VECTORS IN T+1
X[,t] <- round(N[,t]*(1-h[,t]))
#############################################################
# DERIVED HARVEST NUMBERS
H[,t] <- round(N[,t]*h[,t])
X_tot[t] <- sum(X[,t])
N_tot[t] <- sum(N[,t])
if(error.dist=="pois"){
N_obs[t] <- rpois(1, N_tot[t])
}
if(error.dist=="bin"){
N_obs[t] <- rbinom(1, N_tot[t], p.obs)
}
if(error.dist=="count"){
N_obs[t] <- N_tot[t]
}
#############################################################
# DERIVED OBSERVATIONS
# Calf surveys in spring
J[t] <- rbinom(1, round((N[1,t]+N[2, t])/phi1[t]), pa[t])
SU[t] <- rbinom(1, (N[3, t]+N[4, t]+N[5, t]), pa[t]*bias1)
##############################################################
# Structural surveys - post harvest
C0[t] <- rbinom(1, (X[1,t]+X[2, t]), pb[t])
Cf[t] <- rbinom(1, (X[3,t]+X[5, t]), pb[t]*bias2)
Cm[t] <- rbinom(1, (X[4,t]+X[6, t]), pb[t]*bias3)
}
out <- list(X, N, N_juv, h, H, X_tot, N_tot, N_obs, J, SU, C0, Cf, Cm, pa, pb, f, phi1, PHI3)
names(out) <- c("X", "N", "N_juv", "h", "H", "X_tot", "N_tot","N_obs", "J", "SU", "C0", "Cf", "Cm", "p1", "p2", "f", "phi1", "phi3")
out
}