F.cjs.simulate {mra} | R Documentation |
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 probabilities are randomly 'born' into the realized population and captured. Any CJS model, including those with heterogeneous survival or capture probabilities, 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.
|
super.s |
A matrix or vector of true survival probabilities in the super-population 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 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 heterogeneous model is simulated. In this example, the
realized population is sampled for 10 occasions.
Non-equal sized heterogeneous 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, approximately 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 heterogeneous capture probabilities proportional to a vector w
can be simulated using
F.cjs.simulate( matrix( w/sum(x), length(w), 10), rep(s,9) )
.
A stochastic 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 guarantee that a certain number of individuals in the realized population have capture probabilities 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 stochastic 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@west-inc.com)
See Also
Examples
## 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,ns-1), 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 heterogeniety 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")
## End(Not run)