Positive Bernoulli Family Function with 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). Capture history behavioural effects are accommodated.
posbernoulli.b(link = "logitlink", drop.b = FALSE ~ 1,
type.fitted = c("likelihood.cond", "mean.uncond"), I2 = FALSE,
ipcapture = NULL, iprecapture = NULL,
p.small = 1e-4, no.warning = FALSE)link, drop.b, ipcapture, iprecapture |
See |
I2 |
Logical.
This argument is used for terms that are not parallel.
If |
type.fitted |
Details at |
p.small, no.warning |
See |
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 x 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.
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm,
and vgam.
The dependent variable is not scaled to row proportions.
This is the same as posbernoulli.t
and posbernoulli.tb
but different from posbinomial
and binomialff.
Thomas W. Yee.
See posbernoulli.t.
posbernoulli.t and
posbernoulli.tb (including estimating N),
deermice,
dposbern,
rposbern,
posbinomial,
aux.posbernoulli.t,
prinia.
# 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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.