posbernoulli.b(link = "logit", drop.b = FALSE ~ 1,
type.fitted = c("likelihood.cond", "mean.uncond"), I2 = FALSE,
ipcapture = NULL, iprecapture = NULL,
p.small = 1e-4, no.warning = FALSE)
CommonVGAMffArguments
for information about
these arguments.
By default the parallelism assumption does not apply to the intercept.
With an intercept-only model
setting drop.
TRUE
then the constraint matrix diag(2)
(the general default constraint matrix in cbind(0:1, 1)
. The latter meposbernoulli.tb
.posbernoulli.t
."vglmff"
(see vglmff-class
).
The object is used by modelling functions such as vglm
,
and vgam
.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 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.
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)
Run the code above in your browser using DataLab