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) }