GenerateCRData: Simulate capture-recapture data according to the recruitment...

Description Usage Arguments Examples

View source: R/GenerateCRData.R

Description

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.

Usage

1
d <- GenerateCRData(T, E, N, t.bar, sigma, alpha.0, alpha.1 = 0, beta.0, beta.1)

Arguments

T

Vector of length J containing sampling days.

E

Vector of length J containing effort during each sampling day.

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.

Examples

  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))
  }

shanearichards/crrecruit documentation built on May 26, 2017, 7:01 a.m.