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!!
Run the code above in your browser using DataLab