posbernoulli.tb: Positive Bernoulli Family Function with Time and Behavioural...
In VGAM: Vector Generalized Linear and Additive Models

posbernoulli.tb

R Documentation

Positive Bernoulli Family Function with Time and 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). 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.

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.

Author(s)

Thomas W. Yee.

References

See posbernoulli.t.

Examples

## Not run: 
# 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)

## End(Not run)

VGAM documentation built on Sept. 18, 2024, 9:09 a.m.