View source: R/family.positive.R
posbernoulli.b | R Documentation |
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 \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.
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
.
## 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 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 logit scale
# (2,1) elt is for the behavioural effect:
constraints(M.bh)[["(Intercept)"]]
summary(M.bh, presid = FALSE) # Significant trap-happy 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, nTimePts=nTimePts, pvars=2, is.popn=TRUE)
nrow(pdata) # < 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 # Population size estimate; 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)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.