huggins.sim <- function(N, t, x=NULL, b=c(-3,1) ) {
# Simulate a vector of covariates and a capture history
# N = # in population
# t = number of capture occasions
# x, if provided, is a prespecified vector of covariates
# b is the vector of regression intercept and regr. slope
# result is a matrix with t+1 columns.
# First col is X; rest is capture history
if (is.null(x)) {
x <- runif(N, 0, 6) # N covariates unif between 0 and 6
}
p.capt <- 1/(1+exp(-(b[1]+b[2]*x))); # P[capture each individ.
cap.hist <- matrix(runif(N*t),nrow=N,ncol=t)
cap.hist <- cap.hist < p.capt # calculate whether individ seen each time
# result is a matrix of 0's and 1's where P[1] = p.capt[i]
# R recycles columns, so above is exactly the same as:
# for (j in 1:t) {
# cap.hist[,j] <- runif(N) < p.capt
# }
nseen <- apply(cap.hist,1,sum) # total # times seen
cbind(x,cap.hist)[nseen > 0,]
# combine x and cap hist and return rows where animal seen at least once
}
huggins.lnl <- function(beta, data) {
# calculate lnl for Huggins with individual heterogeneity dep. on 1 X
# data is a matrix of N seen x (Ntimes + 1)
# first col is the vector of covariates, rest is the capture history
x <- data[,1]
cap.hist <- data[,-1]
t <- dim(cap.hist)[2]
pcap <- 1/(1+exp(-(beta[1]+beta[2]*x)));
pseen <- 1-(1-pcap)^t
lnli <- apply(cap.hist*log(pcap) + (1-cap.hist)*log(1-pcap),1,sum) - log(pseen)
sum(lnli)
}
huggins.n <- function(beta, data) {
# estimate n and var n given estimates of beta and a data matrix
# columns of data are same as for huggins.lnl
x <- data[,1]
cap.hist <- data[,-1]
t <- dim(cap.hist)[2]
pcap <- 1/(1+exp(-(beta[1]+beta[2]*x)));
pseen <- 1-(1-pcap)^t
n <- sum(1/pseen)
var.n <- sum((1-pseen)/(pseen*pseen))
c(n=n, var=var.n)
}