Description Usage Arguments Examples
View source: R/GenerateCRData.R
Generates a capture-recapture data set that is consistent with the recruitment model. The function returns a list that includes summarises of the data needed when fitting the model. This function is useful for testing the fitting approach.
1 | d <- GenerateCRData(T, E, N, t.bar, sigma, alpha.0, alpha.1 = 0, beta.0, beta.1)
|
T |
Vector of length |
E |
Vector of length |
N |
The total number of recruits that are availbale for capture during the survey. |
t.bar |
The day that recruitment peaks. |
sigma |
Duration of the recruitment period: approximately 95% of recruitment occurs over 4*sigma days. |
alpha.0 |
Baseline catch rate. |
alpha.1 |
Age-dependent change in catch rate. |
beta.0 |
Baseline mortality rate. |
beta.1 |
Age-dependent change in mortality rate. |
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 | # Summary of survey (days of sampling and effort during each survey)
T <- c(1, 4, 7, 10, 13, 16, 19, 22, 25, 28) # sampling days
E <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1) # sampling hours per event
# True parameter values
N <- 250 # total number of recruiting animals
t.bar <- 3.0 # peak day of recruitment
sigma <- 3.0 # spread in recruitment
alpha.0 <- 0.15 # baseline catchability
alpha.1 <- 0.0 # age-depependent catchability
beta.0 <- 0.03 # baseline mortality rate
beta.1 <- 0.10 # age-dependent mortality rate
# Generate some capture-recapture data
Simulated <- GenerateCRData(T, E, N, t.bar, sigma, alpha.0, alpha.1,
beta.0, beta.1)
# Plot the true number of catchable animals for each survey event
plot(x = T, y = Simulated$N.catchable.obs,
xlab = "Day", ylab = "Catchable animals")
# set up necessary global variables for llRecruit
y <- Simulated$y
f <- Simulated$f
l <- Simulated$l
T.F <- Simulated$T.F
T.L <- Simulated$T.L
# Can the true parameters be recovered?
guess <- c(t.bar, sigma, alpha.0, alpha.1, beta.0, beta.1)
fit <- optim(par = guess, fn = llRecruit,
lower = c(0, 1.0, 0.10, -0.1, 0.0025, 0.001),
upper = c(6, 5.0, 0.35, 0.1, 0.1000, 0.190),
method = "L-BFGS-B",
control = list(fnscale = -1, trace = 1, maxit = 100, REPORT = 20))
fit$par # display maximum-likelihood parameter estimates
I <- dim(y)[1] # number of unique animals caught
I/(1-PrNotDetect(fit$par)) # estimate of population size (N)
# compare true and estimated values
plot(x = guess, y = fit$par, xlab = "True value", ylab = "Fitted value")
abline(a = 0, b = 1, col = "grey")
## The function is currently defined as
function (T, E, N, t.bar, sigma, alpha.0, alpha.1 = 0, beta.0,
beta.1)
{
last.sample.day <- max(T)
J <- length(E)
T.F <- min(T) - 14
T.L <- max(T)
Y.days <- last.sample.day - T.F + 1
Y.t <- T.F:last.sample.day
recruit.day <- seq(from = T.F, to = T.L, by = 1)
u.True <- exp(-0.5 * ((recruit.day - t.bar)/sigma)^2)
u.True <- u.True/sum(u.True)
U <- cumsum(u.True)
age <- 0:Y.days
if (beta.1 != 0) {
pr.alive <- exp(-(beta.0/beta.1) * (exp(beta.1 * age) -
1))
}
else {
pr.alive <- exp(-beta.0 * age)
}
pr.dead <- 1 - pr.alive
r.1 <- runif(n = N, min = 0, max = 1)
r.2 <- runif(n = N, min = 0, max = 1)
day.recruited <- rep(0, N)
day.died <- rep(0, N)
longevity <- rep(0, N)
for (i in 1:N) {
day.recruited[i] <- sum(r.1[i] > U) + T.F
longevity[i] <- sum(r.2[i] >= pr.dead) - 1
day.died[i] <- day.recruited[i] + longevity[i]
}
Y <- matrix(data = 0, nrow = N, ncol = Y.days)
for (j in 1:Y.days) {
t <- T.F - 1 + j
for (i in 1:N) {
if ((t >= day.recruited[i]) & (t <= day.died[i])) {
Y[i, j] <- 1
}
}
}
N.catchable <- rep(0, Y.days)
for (j in 1:Y.days) {
N.catchable[j] <- sum(Y[, j])
}
N.catchable.obs <- rep(0, Y.days)
for (j in 1:Y.days) {
if (Y.t[j] %in% T) {
N.catchable.obs[j] <- sum(Y[, j])
}
else {
N.catchable.obs[j] <- NA
}
}
N.catchable.obs <- N.catchable.obs[!is.na(N.catchable.obs)]
Y.obs <- NULL
age.obs <- NULL
for (j in 1:Y.days) {
if (Y.t[j] %in% T) {
Y.obs <- cbind(Y.obs, Y[, j])
age.obs <- cbind(age.obs, Y.t[j] - day.recruited)
}
}
y <- Y.obs
for (j in 1:J) {
for (i in 1:N) {
if (Y.obs[i, j] == 1) {
pr.capture <- 1 - exp(-E[j] * alpha.0 * exp(alpha.1 *
age.obs[i, j]))
if (runif(1) <= pr.capture) {
y[i, j] <- 1
}
else {
y[i, j] <- 0
}
}
}
}
y.tmp <- y
y <- NULL
for (i in 1:N) {
if (sum(y.tmp[i, ]) > 0) {
y <- rbind(y, y.tmp[i, ])
}
}
I <- dim(y)[1]
f <- rep(0, I)
l <- rep(0, I)
for (i in 1:I) {
found.yet <- FALSE
for (j in 1:J) {
if (y[i, j] == 1) {
if (!found.yet) {
f[i] <- j
found.yet <- TRUE
}
l[i] <- j
}
}
}
return(list(y = y, I = I, J = J, f = f, l = l, T.F = T.F,
T.L = T.L, u.True = u.True, pr.alive = pr.alive, N.catchable = N.catchable,
N.catchable.obs = N.catchable.obs))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.