simHMM: Simulates data from Hidden Markov Model

In marked: Mark-Recapture Analysis for Survival and Abundance Estimation

Simulates data from Hidden Markov Model

Description

Creates a set of data from a specified HMM model for capture-recapture data.

Usage

simHMM(
  data,
  ddl = NULL,
  begin.time = 1,
  model = "hmmCJS",
  title = "",
  model.parameters = list(),
  design.parameters = list(),
  initial = NULL,
  groups = NULL,
  time.intervals = NULL,
  accumulate = TRUE,
  strata.labels = NULL
)

Arguments

data	Either the raw data which is a dataframe with at least one column named ch (a character field containing the capture history) or a processed dataframe
ddl	Design data list which contains a list element for each parameter type; if NULL it is created
begin.time	Time of first capture(release) occasion
model	Type of c-r model
title	Optional title; not used at present
model.parameters	List of model parameter specifications
design.parameters	Specification of any grouping variables for design data for each parameter
initial	Optional list (by parameter type) of initial values for beta parameters (e.g., initial=list(Phi=0.3,p=-2)
groups	Vector of names of factor variables for creating groups
time.intervals	Intervals of time between the capture occasions
accumulate	if TRUE, like capture-histories are accumulated to reduce computation
strata.labels	labels for strata used in capture history; they are converted to numeric in the order listed. Only needed to specify unobserved strata. For any unobserved strata p=0..

Details

The specification for the simulation includes a set of data with at least 2 unique ch and freq value to specify the number of ch values to simulate that start at the specified occasion. For example, 1000 50 0100 50 0010 50 would simulate 150 capture histories with 50 starting at each of occasions 1 2 and 3. The data can also contain other fields used to generate the model probabilities and each row can have freq=1 to use individual covariates. Either a dataframe (data) is provided and it is processed and the design data list are created or the processed dataframe and design data list are provided. Formula for the model parameters for generating the data are provided in model.parameters and parameter values are provided in initial.

Value

dataframe with simulated data

Author(s)

Jeff Laake

Examples

# simulate phi(.) p(.) with 1000 Females and 100 males, 3 occasions all released on first occasion
df=simHMM(data.frame(ch=c("100","110"),sex=factor(c("F","M")),freq=c(1000,100),
   stringsAsFactors=FALSE))
df=simHMM(data.frame(ch=rep("100",100),u=rnorm(100,0,1),freq=rep(1,100),
  stringsAsFactors=FALSE),
  model.parameters=list(Phi=list(formula=~u),p=list(formula=~time)),
   initial=list(Phi=c(1,1),p=c(0,1)))
df=simHMM(data.frame(ch=c("1000","0100","0010"),freq=rep(50,3),stringsAsFactors=FALSE),
  model.parameters=list(Phi=list(formula=~1),p=list(formula=~time)),
    initial=list(Phi=c(1),p=c(0,1,2)))

######################################################################## 
# Example developed by Jay Rotella, Montana State University
# example of using the 'simHMM' function in 'marked' package for
#  a multi-state model with 3 states

# simulate a single release cohort of 1000 animals with 1 release from 
#  each of the 3 states for 10 recapture occasions; 
# Note: at least 2 unique ch are needed in simHMM
simd <- data.frame(ch = c("A0000000000", "B0000000000", "C0000000000"),
	freq = c(1000, 1000, 1000),
	stringsAsFactors = FALSE)

# define simulation/fitting model; default for non-specified parameters is ~1
# but all are listed here. For the formula for Psi, the -1 is necessary to remove
# the intercept which is not needed and would be redundant. 
# The formula '~ -1 + stratum:tostratum' is often appropriate for multi-state
# transitions as it allows different values for Psi depending on 
# the current state (state at time t) and the new state (state at t+1). 
modelspec <- list(
	S = list(formula = ~ -1 + stratum),
	# by presenting the formula for S as '~-1+stratum', it is possible to 
	# present the desired survival rates directly for each group in the 
	# simulations (see how 'initial' is created below) for this simple
	# scenario where survival varies by stratum but no other covariates
	p = list(formula = ~ 1),
	Psi = list(formula =  ~ -1 + stratum:tostratum))

# process data with A,B,C strata
sd <- process.data(simd,
	model = "hmmMSCJS",
	strata.labels = c("A", "B", "C"))

# create design data
ddl <- make.design.data(sd)
# view design data (especially important for seeing which order to present
#  beta values used to set probabilities for various transitions)
head(model.matrix(~stratum, ddl$S))
head(model.matrix(~-1+stratum:tostratum, ddl$Psi), 9)

# set initial parameter values for model S=~stratum, p=~1, Psi=~stratum:tostratum
initial <- list(S = c(log(0.8/0.2), log(0.6/0.4), log(0.5/0.5)),
	p = log(0.6/0.4),
	# order of presentation is BA, CA, AB, CB, AC, BC
	# Note: values for AA, BB, CC are obtained by subtraction
	Psi = c(log(0.2/0.6), log(0.3/0.1),   # Psi(BA) & Psi(CA)
			log(0.3/0.5), log(0.6/0.1),   # Psi(AB) & Psi(CB)
			log(0.2/0.5), log(0.2/0.6)))  # Psi(AC) & Psi(BC)
# The desired probabilties for the transition matrix for Psi are as follows
#        To:
# From:    A   B    C
#    A   .5   .3   .2
#    B   .2   .6   .2
#    C   .3   .6   .4
# The values in the list above are obtained using the transitions AA, BB, CC
#  as reference values in the denominator of log-odds calculations. 
# For example, to achieve the desired probabilities for the transition 
#  from A to B use log(0.3/0.5) and from A to C use log(0.2/0.5)
# exp(log(0.3/0.5))/(1 + exp(log(0.3/0.5)) + exp(log(0.2/0.5))) = 0.3
# exp(log(0.2/0.5))/(1 + exp(log(0.3/0.5)) + exp(log(0.2/0.5))) = 0.2
# Note:           1/(1 + exp(log(0.3/0.5)) + exp(log(0.2/0.5))) = 0.5

# call simmHMM to get a single realization
realization <- simHMM(sd, ddl, model.parameters = modelspec, 
		initial = initial)

# using that realization, process data and make design data
# note that the analysis can also use model="MSCJS" in marked or "Multistrata" in RMark
# it is only the simulation that requires specification as an HMM.
sd <- process.data(realization, model = "hmmMSCJS",
	strata.labels = c("A","B","C"))
rddl <- make.design.data(sd)

# fit model
m <- crm(sd, rddl, model.parameters = modelspec,
	hessian = TRUE)

# model output
m$results$beta
initial
m$results$reals

marked documentation built on Oct. 19, 2023, 5:06 p.m.