F.cjs.simulate: F.cjs.simulate - Generation of capture histories that follow a CJS model.

Description

This function generates capture history matrices that follow open Cormack-Jolly-Seber (CJS) models. A super-population approach is taken wherein individuals with unique capture and survival probabilites are randomly 'born' into the realized population and captured. Any CJS model, including those with heterogeneity, can be simulated. Closed populations can also be simulated.

Usage

F.cjs.simulate(super.p, super.s, fit, N1 = 1000, births.per.indiv = "constant.popln", R = 100)

Arguments

super.p

A matrix or vector of true capture probabilities in the super-population of individuals.

Ifsuper.pis a VECTOR, all individuals in the realized population will have the same true capture probabilities, but c

super.s

A matrix or vector of true survival probabilities in the super-population of individuals.

Ifsuper.sis a VECTOR, all individuals in the realized population will have the same true survival probabilities aft

fit

A previously estimated CJS object. Instead of specifying super.p and super.s, a fitted CJS model can be specified. If either one of super.p or super.s is missing, the (estimated) pro

A scalar specifying the initial population size. I.e., N1 individuals will be 'born' into the realized population just before the first sampling occasion.

births.per.indiv

Either a vector of births per individual in the realized population, or the string "constant.popln" (the default). If births.per.indiv = "constant.popln", the total number of births into the realized population between ca

A scalar specifying the number of replications for the simulation. A total of R independent capture history matricies will be generated.

Value

A list of length R. Each component of this list is a list of length 2. Each of these R sublists contains the following components:
histsThe simulated capture histories for a particular iteration. This is a matrix with a random number of rows (due to the stocastic nature of captures) and NS columns.
popln.nA vector of length NS containing the true population sizes at each sampling occasion.

Details

Some examples: A two-group heterogeneous population contains one group of individuals with one common set of capture probabilities, and another group of individuals with another set of common capture probabilities. A population with one group of individuals having capture probability equal to 0.25, and another group with capture probability equal to 0.75 can be simulated using

F.cjs.simulate( rbind( rep(0.25,10),rep(0.75,10) ), rep(s,9) ).

, where s is some survival probability between 0 and 1. If s = 1, a closed (no births or deaths) two-group heterogeneity model is simulated. In this example, the realized population is sampled for 10 occasions. Non-equal sized hetergeneity groups can be simulated using

F.cjs.simulate( rbind( matrix(0.25,1,10),matrix(0.75,9,10) ), rep(1,9) ).

Using this call, approximatley 10% of individuals in the realized population will have capture probabilities equal to 0.25, while 90% will have capture probabilities equal to 0.75. Additional groups can be included by including more rows with distinct probabilities in super.p. A population with capture heterogeneity proportional to a vector w can be simulated using

F.cjs.simulate( matrix( w/sum(x), length(w), 10), rep(s,9) )

. A stocastic population that varies around a specified size of N1 = 1000 can be simulated with a statement like

F.cjs.simulate( rep(0.25,10), rep(s,9), N1=1000, births.per.indiv=rep((1-s)/s,9) ).

In this simulation, N(j)*(1-s) individuals die between each occasion, but are replaced because the N(j)*s surviving individuals each have (1-s)/s offspring. Because of the super-population approach taken here, it is not possible to specify which individuals have which survival or capture probabilities, nor to guarentee that a certain number of individuals in the realized population have capture probabilites equal to any particular value.

Examples

Run this code

## Careful: These examples can take a long time to run. 

## Simulate constant model, and analyze

ns <- 10
N <- 100
sim.list <- F.cjs.simulate( rep(0.3,ns), rep(0.9,ns-1), N1=N, R=100 )  # should do at least R = 500, but takes a while...

f.analyze <- function(x){
    fit <- F.cjs.estim( ~1, ~1, x$hists, control=mra.control(maxfn=200, cov.meth=2) )
    if( fit$exit.code == 1 ){
        return( fit$n.hat )
    } else {
        return( rep(NA,ncol(x$hists)) )
    }
}
results <- t(sapply(sim.list, f.analyze))
plot( 1:10, colMeans(results, na.rm=TRUE), xlab="Occasion", ylab="Mean population estimate", col="red", type="b")
abline( h=N )

## Plot RMSE by occasion
std <- apply(results, 2, sd, na.rm=TRUE)
bias <- apply(results - N, 2, mean, na.rm=TRUE)
plot( std, bias, type="n" )
text( std, bias, 2:10 )
abline(h=0)
title(main="RMSE by Sample Occasion")

## Show bias when heterogeneity is present
sim.list <- F.cjs.simulate( matrix(c(0.3,.7,.7,.7),4,ns), rep(0.9,ns-1), N1=N, R=100 ) 
results <- t(sapply(sim.list, f.analyze))
mean.N <- colMeans(results, na.rm=TRUE)
plot( 1:length(mean.N), mean.N, ylim=range(c(mean.N,N),na.rm=TRUE), xlab="Occasion", ylab="Mean population estimate", col="red", type="b")
abline( h=N )
abline( h=mean(mean.N), col="red", lty=2)
title(main="Heterogeniety causes negative bias")


## Simulate CJS model, first estimate one
data(dipper.histories)
ct <- as.factor( paste("T",1:ncol(dipper.histories), sep=""))
attr(ct,"nan")<-nrow(dipper.histories)
dipper.cjs <- F.cjs.estim( ~tvar(ct,drop=c(1,2)), ~tvar(ct,drop=c(1,6,7)), dipper.histories )

## Now generate histories from it.
sim.list <- F.cjs.simulate( fit=dipper.cjs, N1=100, birth.rate=rep(1,6), R=100 )

## Now analyze generated histories using lapply or sapply.  Can fit any model. 
## Here we fit the correct model.
f.analyze <- function(x){
    #   write a counter to console, this is not necessary
    i <- get("i", env=.GlobalEnv) + 1
    cat(paste("Iteration", i, "\n"))
    assign("i",i,env=.GlobalEnv)
    
    ct <- as.factor( 1:ncol(x$hists) )
    fit <- F.cjs.estim( ~tvar(ct,nan=nrow(x$hists),drop=c(1,2)), ~tvar(ct,nan=nrow(x$hists),drop=c(1,6,7)), 
        x$hists, control=mra.control(maxfn=200, cov.meth=2) )
    if( fit$exit.code == 1 ){
        return( fit$n.hat )
    } else {
        return( rep(NA,ncol(x$hists)) )
    }
}
i <- 0
results <- t(sapply(sim.list, f.analyze))
mean.N <- colMeans(results, na.rm=TRUE)
plot( 1:length(mean.N), mean.N, ylim=range(c(mean.N,N),na.rm=TRUE), xlab="Occasion", ylab="Mean population estimate", col="red", type="b")
abline( h=N )
title(main="Time varying CJS model")