F.cjs.simulate  Generation of capture histories that follow a CJS model.
Description
This function generates capture history matrices that follow open CormackJollySeber (CJS) models. A superpopulation 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
1 2  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 superpopulation of individuals.

super.s 
A matrix or vector of true survival probabilities in the superpopulation of individuals.
Number of survival probabilities in 
fit 
A previously estimated CJS object. Instead of specifying 
N1 
A scalar specifying the initial population size. I.e., 
births.per.indiv 
Either a vector of births per individual in the realized population, or
the string "constant.popln" (the default). If

R 
A scalar specifying the number of replications for the simulation. A total of 
Details
Some examples: A twogroup 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) twogroup heterogeneity model is simulated. In this example, the
realized population is sampled for 10 occasions.
Nonequal 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((1s)/s,9) ).
In this simulation, N(j)*(1s) individuals die between each occasion, but are replaced because the N(j)*s surviving individuals each have (1s)/s offspring.
Because of the superpopulation 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.
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:
hists 
The 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.n 
A vector of length NS containing the true population sizes at each sampling occasion. 
Author(s)
Trent McDonald, WEST Inc. (tmcdonald@westinc.com)
See Also
F.cjs.estim
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79  ## Not run:
## Don't run specified because these examples can take > 10 seconds.
## Simulate constant model, and analyze
ns < 10
N < 100
sim.list < F.cjs.simulate( rep(0.3,ns), rep(0.9,ns1), N1=N, R=100 )
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,ns1), 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")
## End(Not run)
