VGAM (version 0.9-1)

posbernoulli.tb: Positive Bernoulli Family Function with Time and Behavioural Effects (experimental)

Description

Fits a GLM-like model to multiple (currently only two or three) Bernoulli responses where each row in the capture history matrix response has at least one success (capture). Sampling occasion effects and behavioural effects are accommodated. However, this function only handles two and three sampling occasions.

Usage

posbernoulli.tb(link = "logit", parallel.t = FALSE, parallel.b = FALSE,
                apply.parint = FALSE, imethod = 1, iprob = NULL,
                dconst = 0.1, dpower = -2)

Arguments

link, imethod, iprob, parallel.t, parallel.b, apply.parint

See CommonVGAMffArguments for information. But parallel.t and parallel.b must each be logicals only. Argument parallel.t means parallel with respec

dconst, dpower

Decay constants and power (exponent) for the ridge adjustment for the working weight matrices. At iteration $t$ of the IRLS algorithm a positive value is added to the first $\tau$ diagonal elements of the working weight matrices to make them pos

Value

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

Warning

As this model is likely to be overparameterized, probably this function should not be used (for now?). Estimation for the population size (and its SE) for the $M_{tb}$ model may be wrong. Models $M_{tbh}$ and $M_{th}$ may be wrong. But models $M_{bh}$, $M_{h}$, $M_{b}$, $M_{t}$, $M_{0}$ seem fine.

Inference, especially using standard errors, may be fraught here because the EIM is, strictly speaking, not of full rank. A similar adjustment is made by zipebcom. It is a good idea to monitor convergence. The $M_0$ model is best fitted with posbernoulli.b or posbernoulli.t or posbinomial because the standard errors are more accurate.

Details

This model (commonly known as $M_{tb}$ in the capture--recapture literature) operates on a response matrix of 0s and 1s. See posbernoulli.t for information that is in common.

This VGAM family function is experimental only. When finished, it should allow time and behavioural effects to be modelled. Evidently, the expected information matrix (EIM) is not of full rank, so dconst and dpower are used to try fix up the problem. The default link functions are $(logit \,p_{c1},\ldots,logit \,p_{c\tau},logit \,p_{r2},\ldots,logit \,p_{r\tau})^T$ where the subscript $c$ denotes capture, the subscript $r$ denotes recapture, and it is not possible to recapture the animal at sampling occasion 1. Thus $M = 2\tau - 1$. The parameters are currently prefixed by cap.prob and recap.prob for the capture and recapture probabilities.

References

See posbernoulli.t.

See Also

posbernoulli.b (including $\widehat{N}$), posbernoulli.t, posbinomial.

Examples

# Example 1: simulated data
set.seed(123)
nTimePts <- 2  # Must be 2 or 3 currently (aka tau == # of sampling occasions)
nnn <- 10000   # Number of animals
pdata <- rposbern(n = nnn, nTimePts = nTimePts, pvars = 2)
dim(pdata)
head(pdata)

clist <- list("(Intercept)" = cbind(1, c(0, 0, 1)),  # Capture effect is last coln
              x2            = rbind(1, 1, 1))
M_tbh.1 <- vglm(cbind(y1, y2) ~ x2,
                constraints = clist, trace = TRUE,
                posbernoulli.tb, data = pdata)
summary(M_tbh.1)

coef(M_tbh.1)
coef(M_tbh.1, matrix = TRUE)
constraints(M_tbh.1, matrix = TRUE)
summary(M_tbh.1)  # Standard errors are very approximate
head(fitted(M_tbh.1))
head(model.matrix(M_tbh.1, type = "vlm"), 21)
dim(depvar(M_tbh.1))


# Example 2: Perom subset data
Hlist <- list("(Intercept)" = cbind(1, c(0, 0, 0, 1, 1)),
              sex           = rbind(1, 1, 1, 1, 1),
              weight        = rbind(1, 1, 1, 1, 1))
Psubset <- subset(Perom, y1 + y2 + y3 > 0)
head(Psubset)

fit1 <- vglm(cbind(y1, y2, y3) ~ sex + weight, constraints = Hlist,
             posbernoulli.tb, data = Psubset, trace = TRUE)
coef(fit1)
coef(fit1, matrix = TRUE)
summary(fit1)  # Standard errors are very approximate

# fit1 is the same as Fit1:
Fit1 <- vglm(cbind(y1, y2, y3) ~ sex + weight, data = Psubset,
             posbernoulli.tb(parallel.t = TRUE), trace = TRUE)
constraints(Fit1)  # Same as Hlist

yyy <- depvar(fit1)
if (length(table(4 * yyy[, 1] + 2 * yyy[, 2] + 1 * yyy[, 3])) != 2^(ncol(yyy))-1)
  warning("not every combination is represented by a row in the response matrix")