Positive Bernoulli Family Function with Time and Behavioural Effects

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). Sampling occasion effects and behavioural effects are accommodated.

Usage

posbernoulli.tb(link = "logitlink", parallel.t = FALSE ~ 1,
   parallel.b = FALSE ~ 0, drop.b = FALSE ~ 1,
   type.fitted = c("likelihood.cond", "mean.uncond"),
   imethod = 1, iprob = NULL,
   p.small = 1e-4, no.warning = FALSE,
   ridge.constant = 0.0001, ridge.power = -4)

Arguments

link, imethod, iprob

See CommonVGAMffArguments for information.

parallel.t, parallel.b, drop.b

A logical, or formula with a logical as the response. See CommonVGAMffArguments for information. The parallel.-type arguments specify whether the constraint matrices have a parallelism assumption for the temporal and behavioural effects. Argument parallel.t means parallel with respect to time, and matches the same argument name in posbernoulli.t.

Suppose the model is intercept-only. Setting parallel.t = FALSE ~ 0 results in the \(M_b\) model. Setting drop.b = FALSE ~ 0 results in the \(M_t\) model because it drops columns off the constraint matrices corresponding to any behavioural effect. Setting parallel.t = FALSE ~ 0 and setting parallel.b = FALSE ~ 0 results in the \(M_b\) model. Setting parallel.t = FALSE ~ 0, parallel.b = FALSE ~ 0 and drop.b = FALSE ~ 0 results in the \(M_0\) model. Note the default for parallel.t and parallel.b may be unsuitable for data sets which have a large \(\tau\) because of the large number of parameters; it might be too flexible. If it is desired to have the behaviour affect some of the other covariates then set drop.b = TRUE ~ 0.

The default model has a different intercept for each sampling occasion, a time-parallelism assumption for all other covariates, and a dummy variable representing a single behavioural effect (also in the intercept).

The most flexible model is to set parallel.b = TRUE ~ 0, parallel.t = TRUE ~ 0 and drop.b = TRUE ~ 0. This means that all possible temporal and behavioural effects are estimated, for the intercepts and other covariates. Such a model is not recommended; it will contain a lot of paramters.

type.fitted

Character, one of the choices for the type of fitted value returned. The default is the first one. Partial matching is okay. For "likelihood.cond": the probability defined by the conditional likelihood. For "mean.uncond": the unconditional mean, which should agree with colMeans applied to the response matrix for intercept-only models.

ridge.constant, ridge.power

Determines the ridge parameters at each IRLS iteration. They are the constant and power (exponent) for the ridge adjustment for the working weight matrices (the capture probability block matrix, hence the first \(\tau\) diagonal values). At iteration \(a\) of the IRLS algorithm a positive value is added to the first \(\tau\) diagonal elements of the working weight matrices to make them positive-definite. This adjustment is the mean of the diagonal elements of wz multipled by \(K \times a^p\) where \(K\) is ridge.constant and \(p\) is ridge.power. This is always positive but decays to zero as iterations proceed (provided \(p\) is negative etc.).

p.small, no.warning

See posbernoulli.t.

Details

This model (commonly known as \(M_{tb}\)/\(M_{tbh}\) in the capture–recapture literature) operates on a response matrix of 0s and 1s (\(n \times \tau\)). See posbernoulli.t for information that is in common. It allows time and behavioural effects to be modelled.

Evidently, the expected information matrix (EIM) seems not of full rank (especially in early iterations), so ridge.constant and ridge.power 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 pcapture and precapture for the capture and recapture probabilities. This VGAM family function may be further modified in the future.

Value

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

References

See posbernoulli.t.

Author

Thomas W. Yee.

Note

It is a good idea to apply the parallelism assumption to each sampling occasion except possibly with respect to the intercepts. Also, a simple behavioural effect such as being modelled using the intercept is recommended; if the behavioural effect is not parallel and/or allowed to apply to other covariates then there will probably be too many parameters, and hence, numerical problems. See M_tbh.1 below.

It is a good idea to monitor convergence. Simpler models such as the \(M_0\)/\(M_h\) models are best fitted with posbernoulli.t or posbernoulli.b or posbinomial.

Examples

if (FALSE) { # \dontrun{
# Example 1: simulated data
nTimePts <- 5  # (aka tau == # of sampling occasions)
nnn <- 1000   # Number of animals
pdata <- rposbern(n = nnn, nTimePts = nTimePts, pvars = 2)
dim(pdata); head(pdata)

M_tbh.1 <- vglm(cbind(y1, y2, y3, y4, y5) ~ x2,
                posbernoulli.tb, data = pdata, trace = TRUE)
coef(M_tbh.1)  # First element is the behavioural effect
coef(M_tbh.1, matrix = TRUE)
constraints(M_tbh.1, matrix = TRUE)
summary(M_tbh.1, presid = FALSE)  # Std errors are approximate
head(fitted(M_tbh.1))
head(model.matrix(M_tbh.1, type = "vlm"), 21)
dim(depvar(M_tbh.1))

M_tbh.2 <- vglm(cbind(y1, y2, y3, y4, y5) ~ x2,
                posbernoulli.tb(parallel.t = FALSE ~  0),
                data = pdata, trace = TRUE)
coef(M_tbh.2)  # First element is the behavioural effect
coef(M_tbh.2, matrix = TRUE)
constraints(M_tbh.2, matrix = TRUE)
summary(M_tbh.2, presid = FALSE)  # Std errors are approximate
head(fitted(M_tbh.2))
head(model.matrix(M_tbh.2, type = "vlm"), 21)
dim(depvar(M_tbh.2))

# Example 2: deermice subset data
fit1 <- vglm(cbind(y1, y2, y3, y4, y5, y6) ~ sex + weight,
             posbernoulli.t, data = deermice, trace = TRUE)
coef(fit1)
coef(fit1, matrix = TRUE)
constraints(fit1, matrix = TRUE)
summary(fit1, presid = FALSE)  # Standard errors are approximate

# fit1 is the same as Fit1 (a M_{th} model):
Fit1 <- vglm(cbind(y1, y2, y3, y4, y5, y6) ~ sex + weight,
             posbernoulli.tb(drop.b = TRUE ~ sex + weight,
                parallel.t = TRUE),  # But not for the intercept
             data = deermice, trace = TRUE)
constraints(Fit1)
} # }