posbernoulli.b: Positive Bernoulli Family Function with Behavioural Effects

Description

Fits a GLM-/GAM-like model to multiple Bernoulli responses where each row in the capture history matrix response has at least one success (capture). Capture history behavioural effects are accommodated.

Usage

posbernoulli.b(link = "logit", drop.b = FALSE ~ 1,
               type.fitted = c("likelihood.cond", "mean.uncond"), I2 = FALSE,
               ipcapture = NULL, iprecapture = NULL,
               p.small = 1e-4, no.warning = FALSE)

Arguments

link, drop.b, ipcapture, iprecapture

See CommonVGAMffArguments for information about these arguments. By default the parallelism assumption does not apply to the intercept. With an intercept-only model setting drop.b = TRUE ~ 1 results in the \(M_0\)/\(M_h\) model.

Logical. This argument is used for terms that are not parallel. If TRUE then the constraint matrix diag(2) (the general default constraint matrix in VGAM) is used, else cbind(0:1, 1). The latter means the first element/column corresponds to the behavioural effect. Consequently it and its standard error etc. can be accessed directly without subtracting two quantities.

type.fitted

Details at posbernoulli.tb.

p.small, no.warning

See posbernoulli.t.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

Details

This model (commonly known as \(M_b\)/\(M_{bh}\) in the capture--recapture literature) operates on a capture history matrix response of 0s and 1s (\(n \times \tau\)). See posbernoulli.t for details, e.g., common assumptions with other models. Once an animal is captured for the first time, it is marked/tagged so that its future capture history can be recorded. The effect of the recapture probability is modelled through a second linear/additive predictor. It is well-known that some species of animals are affected by capture, e.g., trap-shy or trap-happy. This VGAM family function does allow the capture history to be modelled via such behavioural effects. So does posbernoulli.tb but posbernoulli.t cannot.

The number of linear/additive predictors is \(M = 2\), and the default links are \((logit \,p_c, logit \,p_r)^T\) where \(p_c\) is the probability of capture and \(p_r\) is the probability of recapture. The fitted value returned is of the same dimension as the response matrix, and depends on the capture history: prior to being first captured, it is pcapture. Afterwards, it is precapture.

By default, the constraint matrices for the intercept term and the other covariates are set up so that \(p_r\) differs from \(p_c\) by a simple binary effect, on a logit scale. However, this difference (the behavioural effect) is more directly estimated by having I2 = FALSE. Then it allows an estimate of the trap-happy/trap-shy effect; these are positive/negative values respectively. If I2 = FALSE then the (nonstandard) constraint matrix used is cbind(0:1, 1), meaning the first element can be interpreted as the behavioural effect.

References

See posbernoulli.t.

Examples

Run this code

# NOT RUN {
# deermice data ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

# Fit a M_b model
M.b <- vglm(cbind(y1, y2, y3, y4, y5, y6) ~ 1,
            posbernoulli.b, data = deermice, trace = TRUE)
coef(M.b)["(Intercept):1"]  # Behavioural effect on the logit scale
coef(M.b, matrix = TRUE)
constraints(M.b, matrix = TRUE)
summary(M.b, presid = FALSE)

# Fit a M_bh model
M.bh <- vglm(cbind(y1, y2, y3, y4, y5, y6) ~ sex + weight,
             posbernoulli.b, data = deermice, trace = TRUE)
coef(M.bh, matrix = TRUE)
coef(M.bh)["(Intercept):1"]  # Behavioural effect on the logit scale
constraints(M.bh)  # (2,1) element of "(Intercept)" is for the behavioural effect
summary(M.bh, presid = FALSE)  # Significant positive (trap-happy) behavioural effect
# Approx. 95 percent confidence for the behavioural effect:
SE.M.bh <- coef(summary(M.bh))["(Intercept):1", "Std. Error"]
coef(M.bh)["(Intercept):1"] + c(-1, 1) * 1.96 * SE.M.bh

# Fit a M_h model
M.h <- vglm(cbind(y1, y2, y3, y4, y5, y6) ~ sex + weight,
            posbernoulli.b(drop.b = TRUE ~ sex + weight),
            data = deermice, trace = TRUE)
coef(M.h, matrix = TRUE)
constraints(M.h, matrix = TRUE)
summary(M.h, presid = FALSE)

# Fit a M_0 model
M.0 <- vglm(cbind(    y1 + y2 + y3 + y4 + y5 + y6,
                  6 - y1 - y2 - y3 - y4 - y5 - y6) ~ 1,
            posbinomial, data = deermice, trace = TRUE)
coef(M.0, matrix = TRUE)
summary(M.0, presid = FALSE)


# Simulated data set ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
set.seed(123); nTimePts <- 5; N <- 1000  # N is the popn size
pdata <- rposbern(n = N, nTimePts = nTimePts, pvars = 2, is.popn = TRUE)
nrow(pdata)  # Less than N (because some animals were never captured)
# The truth: xcoeffs are c(-2, 1, 2) and cap.effect = +1

M.bh.2 <- vglm(cbind(y1, y2, y3, y4, y5) ~ x2,
               posbernoulli.b, data = pdata, trace = TRUE)
coef(M.bh.2)
coef(M.bh.2, matrix = TRUE)
constraints(M.bh.2, matrix = TRUE)
summary(M.bh.2, presid = FALSE)
head(depvar(M.bh.2))    # Capture history response matrix
head(M.bh.2@extra$cap.hist1)  # Info on its capture history
head(M.bh.2@extra$cap1)  # When it was first captured
head(fitted(M.bh.2))     # Depends on capture history
(trap.effect <- coef(M.bh.2)["(Intercept):1"])  # Should be +1
head(model.matrix(M.bh.2, type = "vlm"), 21)
head(pdata)
summary(pdata)
dim(depvar(M.bh.2))
vcov(M.bh.2)

M.bh.2@extra$N.hat     # Estimate of the population size; should be about N
M.bh.2@extra$SE.N.hat  # SE of the estimate of the population size
# An approximate 95 percent confidence interval:
round(M.bh.2@extra$N.hat + c(-1, 1) * 1.96 *  M.bh.2@extra$SE.N.hat, 1)
# }

Run the code above in your browser using DataCamp Workspace