Nothing
#' Simulation tool
#'
#' Simulates movement data from an HMM.
#'
#' @param nbAnimals Number of observed individuals to simulate.
#' @param nbStates Number of behavioural states to simulate.
#' @param stepDist Name of the distribution of the step lengths (as a character string).
#' Supported distributions are: gamma, weibull, lnorm, exp. Default: gamma.
#' @param angleDist Name of the distribution of the turning angles (as a character string).
#' Supported distributions are: vm, wrpcauchy. Set to \code{"none"} if the angle distribution should
#' not be estimated. Default: vm.
#' @param stepPar Parameters of the step length distribution.
#' @param anglePar Parameters of the turning angle distribution.
#' @param beta Matrix of regression parameters for the transition probabilities (more information
#' in "Details").
#' @param covs Covariate values to include in the model, as a dataframe. Default: \code{NULL}.
#' Covariates can also be simulated according to a standard normal distribution, by setting
#' \code{covs} to \code{NULL}, and specifying \code{nbCovs>0}.
#' @param nbCovs Number of covariates to simulate (0 by default). Does not need to be specified of
#' \code{covs} is specified.
#' @param zeroInflation \code{TRUE} if the step length distribution is inflated in zero.
#' Default: \code{FALSE}. If \code{TRUE}, values for the zero-mass parameters should be
#' included in \code{stepPar}.
#' @param obsPerAnimal Either the number of the number of observations per animal (if single value),
#' or the bounds of the number of observations per animal (if vector of two values). In the latter case,
#' the numbers of obervations generated for each animal are uniformously picked from this interval.
#' Default: \code{c(500,1500)}.
#' @param model A moveHMM object. This option can be used to simulate from a fitted model. Default: NULL.
#' Note that, if this argument is specified, most other arguments will be ignored -- except for nbAnimals,
#' obsPerAnimal, covs (if covariate values different from those in the data should be specified),
#' and states.
#' @param states \code{TRUE} if the simulated states should be returned, \code{FALSE} otherwise (default).
#'
#' @return An object moveData, i.e. a dataframe of:
#' \item{ID}{The ID(s) of the observed animal(s)}
#' \item{step}{The step lengths}
#' \item{angle}{The turning angles (if any)}
#' \item{x}{Either easting or longitude}
#' \item{y}{Either northing or latitude}
#' \item{...}{Covariates (if any)}
#'
#' @details \itemize{
#' \item The matrix \code{beta} of regression coefficients for the transition probabilities has
#' one row for the intercept, plus one row for each covariate, and one column for
#' each non-diagonal element of the transition probability matrix. For example, in a 3-state
#' HMM with 2 covariates, the matrix \code{beta} has three rows (intercept + two covariates)
#' and six columns (six non-diagonal elements in the 3x3 transition probability matrix - filled in
#' row-wise).
#' In a covariate-free model (default), \code{beta} has one row, for the intercept.
#'
#' \item If the length of covariate values passed (either through 'covs', or 'model') is not the same
#' as the number of observations suggested by 'nbAnimals' and 'obsPerAnimal', then the series of
#' covariates is either shortened (removing last values - if too long) or extended (starting
#' over from the first values - if too short).
#' }
#'
#' @examples
#' # 1. Pass a fitted model to simulate from
#' # (m is a moveHMM object - as returned by fitHMM - automatically loaded with the package)
#' # We keep the default nbAnimals=1.
#' m <- example$m
#' obsPerAnimal=c(50,100)
#' data <- simData(model=m,obsPerAnimal=obsPerAnimal)
#'
#' # 2. Pass the parameters of the model to simulate from
#' stepPar <- c(1,10,1,5,0.2,0.3) # mean1, mean2, sd1, sd2, z1, z2
#' anglePar <- c(pi,0,0.5,2) # mean1, mean2, k1, k2
#' stepDist <- "gamma"
#' angleDist <- "vm"
#' data <- simData(nbAnimals=5,nbStates=2,stepDist=stepDist,angleDist=angleDist,stepPar=stepPar,
#' anglePar=anglePar,nbCovs=2,zeroInflation=TRUE,obsPerAnimal=obsPerAnimal)
#'
#' stepPar <- c(1,10,1,5) # mean1, mean2, sd1, sd2
#' anglePar <- c(pi,0,0.5,0.7) # mean1, mean2, k1, k2
#' stepDist <- "weibull"
#' angleDist <- "wrpcauchy"
#' data <- simData(nbAnimals=5,nbStates=2,stepDist=stepDist,angleDist=angleDist,stepPar=stepPar,
#' anglePar=anglePar,obsPerAnimal=obsPerAnimal)
#'
#' # step length only and zero-inflation
#' stepPar <- c(1,10,1,5,0.2,0.3) # mean1, mean2, sd1, sd2, z1, z2
#' stepDist <- "gamma"
#' data <- simData(nbAnimals=5,nbStates=2,stepDist=stepDist,angleDist="none",stepPar=stepPar,
#' nbCovs=2,zeroInflation=TRUE,obsPerAnimal=obsPerAnimal)
#'
#' # include covariates
#' # (note that it is useless to specify "nbCovs", which is determined
#' # by the number of columns of "cov")
#' cov <- data.frame(temp=rnorm(500,20,5))
#' stepPar <- c(1,10,1,5) # mean1, mean2, sd1, sd2
#' anglePar <- c(pi,0,0.5,2) # mean1, mean2, k1, k2
#' stepDist <- "gamma"
#' angleDist <- "vm"
#' data <- simData(nbAnimals=5,nbStates=2,stepDist=stepDist,angleDist=angleDist,stepPar=stepPar,
#' anglePar=anglePar,covs=cov)
#'
#' @export
#' @importFrom stats rnorm runif
simData <- function(nbAnimals=1,nbStates=2,stepDist=c("gamma","weibull","lnorm","exp"),
angleDist=c("vm","wrpcauchy","none"),stepPar=NULL,anglePar=NULL,
beta=NULL,covs=NULL,nbCovs=0,zeroInflation=FALSE,obsPerAnimal=c(500,1500),
model=NULL,states=FALSE)
{
##############################
## Check if !is.null(model) ##
##############################
if(!is.null(model)) {
# extract simulation parameters from model
nbStates <- ncol(model$mle$stepPar)
stepDist <- model$conditions$stepDist
angleDist <- model$conditions$angleDist
stepPar <- c(t(model$mle$stepPar))
anglePar <- c(t(model$mle$anglePar))
beta <- model$mle$beta
if(is.null(covs)) {
covsCol <- which(names(model$data)!="ID" & names(model$data)!="x" &
names(model$data)!="y" & names(model$data)!="step" &
names(model$data)!="angle")
covs <- model$data[,covsCol]
if(length(covsCol)>1) {
# remove intercept column, which is not expected in 'covs'
names <- colnames(covs)
covs <- data.frame(covs[,-1]) # data.frame structure is lost when only one column
colnames(covs) <- names[-1]
} else
covs <- NULL
}
# else, allow user to enter new values for covariates
zeroInflation <- model$conditions$zeroInflation
} else {
if(is.null(stepPar))
stop("'stepPar' needs to be specified")
stepDist <- match.arg(stepDist)
angleDist <- match.arg(angleDist)
}
#####################
## Check arguments ##
#####################
stepFun <- paste("r",stepDist,sep="")
angleFun <- paste("r",angleDist,sep="")
if(nbAnimals<1)
stop("nbAnimals should be at least 1.")
if(nbStates<1)
stop("nbStates should be at least 1.")
p <- parDef(stepDist,angleDist,nbStates,TRUE,zeroInflation)
if(length(stepPar)!=p$parSize[1]*nbStates | length(anglePar)!=p$parSize[2]*nbStates) {
error <- "Wrong number of parameters: there should be"
error <- paste(error,p$parSize[1]*nbStates,"step parameters and")
error <- paste(error,p$parSize[2]*nbStates,"angle parameters")
stop(error)
}
if(zeroInflation) {
stepBounds <- p$bounds[1:((p$parSize[1]-1)*nbStates),]
sp <- stepPar[1:(length(stepPar)-nbStates)]
zm <- stepPar[(length(stepPar)-nbStates+1):length(stepPar)]
if(length(which(zm<0 | zm>1))>0)
stop("The zero-mass should be in [0,1].")
} else {
stepBounds <- p$bounds[1:(p$parSize[1]*nbStates),]
sp <- stepPar
}
if(length(which(sp<=stepBounds[,1] | sp>=stepBounds[,2]))>0)
stop(paste("Check the step parameters bounds (the parameters should be",
"strictly between the bounds of their parameter space)."))
if(angleDist!="none") {
# We can't really write distribution-agnostic code here, because the bounds
# defined in parDef are not the actual bounds of the parameter space.
m <- anglePar[1:nbStates] # angle mean
k <- anglePar[(nbStates+1):length(anglePar)] # angle concentration
if(length(which(m<=(-pi) | m>pi))>0)
stop("Check the angle parameters bounds. The angle mean should be in (-pi,pi].")
if(length(which(k<=0))>0)
stop("Check the angle parameters bounds. The concentration should be strictly positive.")
if(angleDist=="wrpcauchy" & length(which(k>=1))>0)
stop("Check the angle parameters bounds. The concentration should be in (0,1).")
}
if(length(which(obsPerAnimal<1))>0)
stop("obsPerAnimal should have positive values.")
if(!is.null(covs) & nbCovs>0) {
if(ncol(covs)!=nbCovs)
warning("covs and nbCovs argument conflicting - nbCovs was set to ncol(covs)")
}
if(!is.null(covs)) {
if(!is.data.frame(covs))
stop("'covs' should be a data.frame")
}
if(!is.null(covs)) {
nbCovs <- ncol(covs)
# account for missing values of the covariates
if(length(which(is.na(covs)))>0)
warning(paste("There are",length(which(is.na(covs))),
"missing covariate values.",
"Each will be replaced by the closest available value."))
for(i in 1:nbCovs) {
if(length(which(is.na(covs[,i])))>0) { # if covariate i has missing values
if(is.na(covs[1,i])) { # if the first value of the covariate is missing
k <- 1
while(is.na(covs[k,i])) k <- k+1
for(j in k:2) covs[j-1,i] <- covs[j,i]
}
for(j in 2:nrow(trackData))
if(is.na(covs[j,i])) covs[j,i] <- covs[j-1,i]
}
}
}
if(length(obsPerAnimal)==1)
obsPerAnimal <- rep(obsPerAnimal,2)
else if(length(obsPerAnimal)!=2)
stop("obsPerAnimal should be of length 1 or 2.")
#######################################
## Prepare parameters for simulation ##
#######################################
# define number of observations for each animal
allNbObs <- rep(NA,nbAnimals)
for(zoo in 1:nbAnimals) {
if(obsPerAnimal[1]!=obsPerAnimal[2])
allNbObs[zoo] <- sample(obsPerAnimal[1]:obsPerAnimal[2],size=1)
else
allNbObs[zoo] <- obsPerAnimal[1]
}
# extend covs if not enough covariate values
if(!is.null(covs)) {
covnames <- colnames(covs)
while(sum(allNbObs)>nrow(covs))
covs <- rbind(covs,covs)
# shrink covs if too many covariate values
covs <- data.frame(covs[1:sum(allNbObs),])
colnames(covs) <- covnames
}
# generate regression parameters for transition probabilities
if(is.null(beta))
beta <- matrix(rnorm(nbStates*(nbStates-1)*(nbCovs+1)),nrow=nbCovs+1)
else if(nrow(beta)!=nbCovs+1 | ncol(beta)!=nbStates*(nbStates-1)) {
if(nbStates>1)
stop(paste("beta should have ",nbCovs+1," rows and ",nbStates*(nbStates-1)," columns.",sep=""))
else
stop("beta should be NULL")
}
# initial state distribution
delta <- rep(1,nbStates)/nbStates
# format parameters
wpar <- n2w(c(stepPar,anglePar),p$bounds,beta,delta,nbStates,estAngleMean=(angleDist!="none"))
par <- w2n(wpar,p$bounds,p$parSize,nbStates,nbCovs,estAngleMean=(angleDist!="none"),stationary=FALSE)
if(zeroInflation) {
zeroMass <- par$stepPar[nrow(par$stepPar),]
stepPar <- par$stepPar[-(nrow(par$stepPar)),,drop=FALSE]
}
else {
zeroMass <- rep(0,nbStates)
stepPar <- par$stepPar
}
anglePar <- par$anglePar # i.e. NULL if angleDist=="none"
trackData <- NULL
allCovs <- NULL
allStates <- NULL
# build the data frame to be returned
data <- data.frame(ID=character(),
step=numeric(),
angle=numeric(),
x=numeric(),
y=numeric())
###########################
## Loop over the animals ##
###########################
for (zoo in 1:nbAnimals) {
# number of observations for animal zoo
nbObs <- allNbObs[zoo]
###############################
## Simulate covariate values ##
###############################
if(nbCovs>0) {
if(is.null(covs)) {
subCovs <- data.frame(cov1=rnorm(nbObs))
if(nbCovs>1) {
for(j in 2:nbCovs) {
c <- data.frame(rnorm(nbObs))
colnames(c) <- paste("cov",j,sep="")
subCovs <- cbind(subCovs,c)
}
}
} else {
# select covariate values which concern the current animal
if(zoo<2)
ind1 <- 1
else
ind1 <- sum(allNbObs[1:(zoo-1)])+1
ind2 <- sum(allNbObs[1:zoo])
subCovs <- data.frame(covs[ind1:ind2,])
if(!is.null(covs))
colnames(subCovs) <- colnames(covs) # keep covariates names from input
}
allCovs <- rbind(allCovs,subCovs)
}
###############################
## Simulate state sequence Z ##
###############################
if(nbStates>1) {
Z <- rep(NA,nbObs)
Z[1] <- sample(1:nbStates,size=1,prob=delta)
for (k in 2:nbObs) {
gamma <- diag(nbStates)
g <- beta[1,]
if(nbCovs==1) g <- g + beta[2,]*subCovs[k,1]
if(nbCovs>1) {
for(j in 1:nbCovs)
g <- g + beta[j+1,]*subCovs[k,j]
}
gamma[!gamma] <- exp(g)
gamma <- t(gamma)
gamma <- gamma/apply(gamma,1,sum)
Z[k] <- sample(1:nbStates,size=1,prob=gamma[Z[k-1],])
}
allStates <- c(allStates,Z)
} else
Z <- rep(1,nbObs)
X <- matrix(nbObs,nrow=nbObs,ncol=2)
X[1,] <- c(0,0) # initial position of animal
phi <- 0
s <- rep(NA,nbObs)
a <- rep(NA,nbObs)
############################
## Simulate movement path ##
############################
for (k in 1:(nbObs-1)){
# prepare lists of arguments for step and angle distributions
stepArgs <- list(1) ; angleArgs <- list(1) # first argument = 1 (one random draw)
for(j in 1:nrow(stepPar))
stepArgs[[j+1]] <- stepPar[j,Z[k]]
if(angleDist!="none") {
for(j in 1:nrow(anglePar))
angleArgs[[j+1]] <- anglePar[j,Z[k]]
}
if(stepDist=="gamma") {
shape <- stepArgs[[2]]^2/stepArgs[[3]]^2
scale <- stepArgs[[3]]^2/stepArgs[[2]]
stepArgs[[2]] <- shape
stepArgs[[3]] <- 1/scale # rgamma expects rate=1/scale
}
if(runif(1)>zeroMass[Z[k]])
s[k] <- do.call(stepFun,stepArgs)
else
s[k] <- 0
if(angleDist!="none" & s[k]>0) {
a[k] <- do.call(angleFun,angleArgs)
if(a[k] > pi) a[k] <- a[k]-2*pi
if(a[k] < -pi) a[k] <- a[k]+2*pi
phi <- phi + a[k]
}
m <- s[k]*c(Re(exp(1i*phi)),Im(exp(1i*phi)))
X[k+1,] <- X[k,] + m
}
a[1] <- NA # the first angle value is arbitrary
# angle = NA if step = 0
step0 <- which(s==0)
a[c(step0,step0+1)] <- NA
d <- data.frame(ID=rep(zoo,nbObs),step=s,angle=a,x=X[,1],y=X[,2])
data <- rbind(data,d)
}
# if covs provided as argument
if(!is.null(covs) & is.null(allCovs))
allCovs <- covs
if(nbCovs>0)
data <- cbind(data,allCovs)
# include states sequence in the data
if(states)
data <- cbind(data,states=allStates)
return(moveData(data))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.