huggins91 {VGAM} | R Documentation |
Fits a Huggins (1991) capture-recapture model to a matrix of 0s and 1s: animals sampled on several occasions and individual animals caught at least once.
huggins91(link = "logit", earg = list(), parallel = TRUE, iprob = NULL, eim.not.oim = TRUE)
link, earg, parallel, iprob |
See |
eim.not.oim |
Logical. If |
This model operates on a response matrix of 0s and 1s. Each of
at least two columns is an occasion where animals are potentially
captured (e.g., a field trip), and each row is an individual
animal. Capture is a 1, else a 0. Each row has at least one
capture. It is well-known that animals are affected by capture,
e.g., trap-shy or trap-happy. This VGAM family function
attempts to allow the capture history to be modelled. This
involves the use of the xij
argument. Ignoring capture
history effects would mean posbinomial
could
be used by aggregating over the sampling occasions.
Huggins (1991) suggests a model involving maximizing a conditional likelihood. The form of this is a numerator divided by a denominator, where the true model has part of the linear/additive predictor modelling capture history applying to the numerator only, so that part is set to zero in the denominator. The numerator of the conditional likelihood corresponds to a sequence of Bernoulli trials, with at least one success, for each animal.
Unfortunately the Huggins model is too difficult to fit in this package, and one can only use the same linear/additive predictor in the numerator as the denominator. Hence this VGAM family function does not implement the model properly.
The number of linear/additive predictors is twice the number of sampling occasions, i.e., M = 2T, say. The first two correspond to the first sampling occasion, the next two correspond to the second sampling occasion, etc. Even-numbered linear/additive predictors should correspond to what would happen if no capture had occurred (they belong to the denominator.) Odd-numbered linear/additive predictors correspond to what actually happened (they belong to the numerator.)
The fitted value for column t is the tth numerator probability divided by the denominator.
An object of class "vglmff"
(see vglmff-class
).
The object is used by modelling functions such as vglm
,
and vgam
.
This VGAM family function is experimental and does not work properly because the linear/additive predictor in the numerator and denominator must be the same. The parameter estimates of the Huggins (1991) model ought to be similar (probably in between, in some sense) to two models: Model 1 is where the capture history variable is included, Model 2 is where the capture history variable is not included. See the example below. A third model, called Model 3, allows for 'half' the capture history to be put in both numerator and denominator. This might be thought of as a compromise between Models 1 and 2, and may be useful as a crude approximation.
Under- or over-flow may occur if the data is ill-conditioned.
The weights
argument of vglm
need not be
assigned, and the default is just a matrix of ones.
This VGAM family function is currently more complicated than it needs to be, e.g., it is possible to simplify M = T, say.
Thomas W. Yee
Huggins, R. M. (1991) Some practical aspects of a conditional likelihood approach to capture experiments. Biometrics, 47, 725–732.
vglm.control
for xij
,
dhuggins91
,
rhuggins91
.
posbinomial
.
set.seed(123); nTimePts = 5 hdata = rhuggins91(n = 1000, nTimePts = nTimePts, pvars = 2) # The truth: xcoeffs are c(-2, 1, 2) and capeffect = -1 # Model 1 is where capture history information is used model1 = vglm(cbind(y1, y2, y3, y4, y5) ~ x2 + Chistory, huggins91, data = hdata, trace = TRUE, xij = list(Chistory ~ ch0 + zch0 + ch1 + zch1 + ch2 + zch2 + ch3 + zch3 + ch4 + zch4 - 1), form2 = ~ 1 + x2 + Chistory + ch0 + ch1 + ch2 + ch3 + ch4 + zch0 + zch1 + zch2 + zch3 + zch4) coef(model1, matrix = TRUE) # Biased!! summary(model1) head(fitted(model1)) head(model.matrix(model1, type = "vlm"), 21) head(hdata) # Model 2 is where no capture history information is used model2 = vglm(cbind(y1, y2, y3, y4, y5) ~ x2, huggins91, data = hdata, trace = TRUE) coef(model2, matrix = TRUE) # Biased!! summary(model2) # Model 3 is where half the capture history is used in both # the numerator and denominator set.seed(123); nTimePts = 5 hdata2 = rhuggins91(n = 1000, nTimePts = nTimePts, pvars = 2, double.ch = TRUE) head(hdata2) # 2s have replaced the 1s in hdata model3 = vglm(cbind(y1, y2, y3, y4, y5) ~ x2 + Chistory, huggins91, data = hdata2, trace = TRUE, xij = list(Chistory ~ ch0 + zch0 + ch1 + zch1 + ch2 + zch2 + ch3 + zch3 + ch4 + zch4 - 1), form2 = ~ 1 + x2 + Chistory + ch0 + ch1 + ch2 + ch3 + ch4 + zch0 + zch1 + zch2 + zch3 + zch4) coef(model3, matrix = TRUE) # Biased!!