Title: | Maximum Likelihood Analysis of Animal Movement Behavior Using Multivariate Hidden Markov Models |
---|---|
Description: | Extended tools for analyzing telemetry data using generalized hidden Markov models. Features of momentuHMM (pronounced ``momentum'') include data pre-processing and visualization, fitting HMMs to location and auxiliary biotelemetry or environmental data, biased and correlated random walk movement models, discrete- or continuous-time HMMs, continuous- or discrete-space movement models, approximate Langevin diffusion models, hierarchical HMMs, multiple imputation for incorporating location measurement error and missing data, user-specified design matrices and constraints for covariate modelling of parameters, random effects, decoding of the state process, visualization of fitted models, model checking and selection, and simulation. See McClintock and Michelot (2018) <doi:10.1111/2041-210X.12995>. |
Authors: | Brett McClintock, Theo Michelot |
Maintainer: | Brett McClintock <[email protected]> |
License: | GPL-3 |
Version: | 2.0.0 |
Built: | 2024-11-11 06:35:06 UTC |
Source: | https://github.com/bmcclintock/momentuHMM |
Adds smoothed gradients to momentuHMMData
or momentuHierHMMData
objects, which can help reduce bias in gradient-based habitat selection coefficients
addSmoothGradient( data, spatialCovs, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scale = NULL )
addSmoothGradient( data, spatialCovs, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scale = NULL )
data |
A |
spatialCovs |
List of |
weights |
Numeric vector indicating the weight for the points in a smoothed version of the gradient, where |
scale |
Numeric scalar indicating the distance from the observed location to each smoothing point (before adjusting by time step for continuous-time HMMs). Default: |
A momentuHMM
or momentuHierHMM
object, with the smooth gradients added for each spatial covariate.
Blackwell, P. G. and J. Matthiopoulos. 2024. Joint inference for telemetry and spatial survey data. Ecology
Akaike information criterion of momentuHMM model(s).
## S3 method for class 'momentuHMM' AIC(object, ..., k = 2, n = NULL)
## S3 method for class 'momentuHMM' AIC(object, ..., k = 2, n = NULL)
object |
A |
... |
Optional additional |
k |
Penalty per parameter. Default: 2 ; for classical AIC. |
n |
Optional sample size. If specified, the small sample correction AIC is used (i.e., |
The AIC of the model(s) provided. If several models are provided, the AICs are output in ascending order.
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m AIC(m) ## Not run: # HMM specifications nbStates <- 2 stepDist <- "gamma" angleDist <- "vm" mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar0 <- c(mu0,sigma0) anglePar0 <- c(-pi/2,pi/2,kappa0) formula <- ~cov1+cov2 # example$m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package mod1 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar0,angle=anglePar0), formula=~1,estAngleMean=list(angle=TRUE)) Par0 <- getPar0(mod1,formula=formula) mod2 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0$Par,beta0=Par0$beta, formula=formula,estAngleMean=list(angle=TRUE)) AIC(mod1,mod2) Par0nA <- getPar0(mod1,estAngleMean=list(angle=FALSE)) mod3 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0nA$Par,beta0=Par0nA$beta, formula=~1) AIC(mod1,mod2,mod3) # add'l models provided as a list using the !!! operator AIC(mod1, !!!list(mod2,mod3)) ## End(Not run)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m AIC(m) ## Not run: # HMM specifications nbStates <- 2 stepDist <- "gamma" angleDist <- "vm" mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar0 <- c(mu0,sigma0) anglePar0 <- c(-pi/2,pi/2,kappa0) formula <- ~cov1+cov2 # example$m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package mod1 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar0,angle=anglePar0), formula=~1,estAngleMean=list(angle=TRUE)) Par0 <- getPar0(mod1,formula=formula) mod2 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0$Par,beta0=Par0$beta, formula=formula,estAngleMean=list(angle=TRUE)) AIC(mod1,mod2) Par0nA <- getPar0(mod1,estAngleMean=list(angle=FALSE)) mod3 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0nA$Par,beta0=Par0nA$beta, formula=~1) AIC(mod1,mod2,mod3) # add'l models provided as a list using the !!! operator AIC(mod1, !!!list(mod2,mod3)) ## End(Not run)
Calculate Akaike information criterion model weights
AICweights(..., k = 2, n = NULL)
AICweights(..., k = 2, n = NULL)
... |
|
k |
Penalty per parameter. Default: 2 ; for classical AIC. |
n |
Optional sample size. If specified, the small sample correction AIC is used (i.e., |
Model objects must all be either of class momentuHMM
or multiple imputation model objects (of class HMMfits
and/or miHMM
).
AIC is only valid for comparing models fitted to the same data. The data for each model fit must therefore be identical. For multiple imputation model objects, respective model fits must have identical data.
The AIC weights of the models. If multiple imputation objects are provided, then the mean model weights (and standard deviations) are provided.
## Not run: # HMM specifications nbStates <- 2 stepDist <- "gamma" angleDist <- "vm" mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar0 <- c(mu0,sigma0) anglePar0 <- c(-pi/2,pi/2,kappa0) formula <- ~cov1+cov2 # example$m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package mod1 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar0,angle=anglePar0), formula=~1,estAngleMean=list(angle=TRUE)) Par0 <- getPar0(mod1,formula=formula) mod2 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0$Par,beta0=Par0$beta, formula=formula,estAngleMean=list(angle=TRUE)) AICweights(mod1,mod2) Par0nA <- getPar0(mod1,estAngleMean=list(angle=FALSE)) mod3 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0nA$Par,beta0=Par0nA$beta, formula=~1) AICweights(mod1,mod2,mod3) # add'l models provided as a list using the !!! operator AICweights(mod1, !!!list(mod2,mod3)) ## End(Not run)
## Not run: # HMM specifications nbStates <- 2 stepDist <- "gamma" angleDist <- "vm" mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar0 <- c(mu0,sigma0) anglePar0 <- c(-pi/2,pi/2,kappa0) formula <- ~cov1+cov2 # example$m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package mod1 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar0,angle=anglePar0), formula=~1,estAngleMean=list(angle=TRUE)) Par0 <- getPar0(mod1,formula=formula) mod2 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0$Par,beta0=Par0$beta, formula=formula,estAngleMean=list(angle=TRUE)) AICweights(mod1,mod2) Par0nA <- getPar0(mod1,estAngleMean=list(angle=FALSE)) mod3 <- fitHMM(example$m$data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0nA$Par,beta0=Par0nA$beta, formula=~1) AICweights(mod1,mod2,mod3) # add'l models provided as a list using the !!! operator AICweights(mod1, !!!list(mod2,mod3)) ## End(Not run)
Used in functions viterbi
, logAlpha
, logBeta
.
allProbs(m)
allProbs(m)
m |
Object |
Matrix of all probabilities.
## Not run: P <- momentuHMM:::allProbs(m=example$m) ## End(Not run)
## Not run: P <- momentuHMM:::allProbs(m=example$m) ## End(Not run)
fitHMM
(or MIfitHMM
) modelPrints parameters with labels based on DM
, formula
, and/or formulaDelta
. See fitHMM
for
further argument details.
checkPar0(data, ...) ## Default S3 method: checkPar0( data, nbStates, dist, Par0 = NULL, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, betaRef = NULL, deltaCons = NULL, stateNames = NULL, optMethod = "nlm", fixPar = NULL, prior = NULL, ... ) ## S3 method for class 'hierarchical' checkPar0( data, hierStates, hierDist, Par0 = NULL, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, optMethod = "nlm", fixPar = NULL, prior = NULL, ... )
checkPar0(data, ...) ## Default S3 method: checkPar0( data, nbStates, dist, Par0 = NULL, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, betaRef = NULL, deltaCons = NULL, stateNames = NULL, optMethod = "nlm", fixPar = NULL, prior = NULL, ... ) ## S3 method for class 'hierarchical' checkPar0( data, hierStates, hierDist, Par0 = NULL, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, optMethod = "nlm", fixPar = NULL, prior = NULL, ... )
data |
|
... |
further arguments passed to or from other methods |
nbStates |
Number of states of the HMM. |
dist |
A named list indicating the probability distributions of the data streams. |
Par0 |
Optional named list containing vectors of state-dependent probability distribution parameters for
each data stream specified in |
beta0 |
Optional matrix of regression coefficients for the transition probabilities. If |
delta0 |
Optional values or regression coefficients for the initial distribution of the HMM. If |
estAngleMean |
An optional named list indicating whether or not to estimate the angle mean for data streams with angular distributions ('vm' and 'wrpcauchy'). |
circularAngleMean |
An optional named list indicating whether to use circular-linear or circular-circular regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. |
formula |
Regression formula for the transition probability covariates. |
formulaDelta |
Regression formula for the initial distribution. |
stationary |
|
mixtures |
Number of mixtures for the state transition probabilities. |
formulaPi |
Regression formula for the mixture distribution probabilities.
Note that only the covariate values from the first row for each individual ID in |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data stream. |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. |
betaCons |
Matrix of the same dimension as |
betaRef |
Numeric vector of length |
deltaCons |
Matrix of the same dimension as |
stateNames |
Optional character vector of length nbStates indicating state names. |
optMethod |
The optimization method to be used. Can be “nlm” (the default; see |
fixPar |
An optional list of vectors indicating parameters which are assumed known prior to fitting the model. |
prior |
A function that returns the log-density of the working scale parameter prior distribution(s). |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). See |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). See |
m <- example$m checkPar0(data=m$data, nbStates=2, dist=m$conditions$dist, estAngleMean = m$conditions$estAngleMean, formula = m$conditions$formula) par <- getPar(m) checkPar0(data=m$data, nbStates=2, dist=m$conditions$dist, estAngleMean = m$conditions$estAngleMean, formula = m$conditions$formula, Par0=par$Par, beta0=par$beta, delta0=par$delta) dummyDat <- data.frame(step=0,angle=0,cov1=0,cov2=0) checkPar0(data=dummyDat, nbStates=2, dist=m$conditions$dist, estAngleMean = m$conditions$estAngleMean, formula = m$conditions$formula) ## Not run: simDat <- simData(nbStates=2, dist=m$conditions$dist, Par = par$Par, spatialCovs = list(forest=forest), centers = matrix(0,1,2), nbCovs = 2) checkPar0(data = simDat, nbStates=2, dist=m$conditions$dist, formula = ~forest, DM = list(step=list(mean=~cov1, sd=~cov2), angle=list(mean=~center1.angle,concentration=~1)), estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE)) par <- list(step=rnorm(8),angle=rnorm(4)) beta0 <- matrix(rnorm(4),2,2) delta0 <- c(0.5,0.5) checkPar0(data = simDat, nbStates=2, dist=m$conditions$dist, Par0 = par, beta0 = beta0, delta0 = delta0, formula = ~forest, DM = list(step=list(mean=~cov1, sd=~cov2), angle=list(mean=~center1.angle,concentration=~1)), estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE)) ## End(Not run)
m <- example$m checkPar0(data=m$data, nbStates=2, dist=m$conditions$dist, estAngleMean = m$conditions$estAngleMean, formula = m$conditions$formula) par <- getPar(m) checkPar0(data=m$data, nbStates=2, dist=m$conditions$dist, estAngleMean = m$conditions$estAngleMean, formula = m$conditions$formula, Par0=par$Par, beta0=par$beta, delta0=par$delta) dummyDat <- data.frame(step=0,angle=0,cov1=0,cov2=0) checkPar0(data=dummyDat, nbStates=2, dist=m$conditions$dist, estAngleMean = m$conditions$estAngleMean, formula = m$conditions$formula) ## Not run: simDat <- simData(nbStates=2, dist=m$conditions$dist, Par = par$Par, spatialCovs = list(forest=forest), centers = matrix(0,1,2), nbCovs = 2) checkPar0(data = simDat, nbStates=2, dist=m$conditions$dist, formula = ~forest, DM = list(step=list(mean=~cov1, sd=~cov2), angle=list(mean=~center1.angle,concentration=~1)), estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE)) par <- list(step=rnorm(8),angle=rnorm(4)) beta0 <- matrix(rnorm(4),2,2) delta0 <- c(0.5,0.5) checkPar0(data = simDat, nbStates=2, dist=m$conditions$dist, Par0 = par, beta0 = beta0, delta0 = delta0, formula = ~forest, DM = list(step=list(mean=~cov1, sd=~cov2), angle=list(mean=~center1.angle,concentration=~1)), estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE)) ## End(Not run)
Computes the standard errors and confidence intervals on the beta (i.e., working) scale of the data stream probability distribution parameters,
as well as for the transition probabilities regression parameters. Working scale depends on the real (i.e., natural) scale of the parameters. For
non-circular distributions or for circular distributions with estAngleMean
=FALSE:
CIbeta(m, alpha = 0.95)
CIbeta(m, alpha = 0.95)
m |
A |
alpha |
Significance level of the confidence intervals. Default: 0.95 (i.e. 95% CIs). |
1) if both lower and upper bounds are finite then logit is the working scale; 2) if lower bound is finite and upper bound is infinite then log is the working scale.
For circular distributions with estAngleMean
=TRUE and no constraints imposed by a design matrix (DM) or bounds (userBounds), then the working parameters
are complex functions of both the angle mean and concentrations/sd natural parameters (in this case, it's probably best just to focus on the real parameter
estimates!). However, if constraints are imposed by DM or userBounds on circular distribution parameters with estAngleMean
=TRUE and circularAngleMean
=FALSE:
1) if the natural bounds are (-pi,pi] then tangent is the working scale, otherwise if both lower and upper bounds are finite then logit is the working scale; 2) if lower bound is finite and upper bound is infinite then log is the working scale.
When circular-circular regression is specified using circularAngleMean
, the working scale
for the mean turning angle is not as easily interpretable, but the
link function is atan2(sin(X)*B,1+cos(X)*B), where X are the angle covariates and B the angle coefficients.
Under this formulation, the reference turning angle is 0 (i.e., movement in the same direction as the previous time step).
In other words, the mean turning angle is zero when the coefficient(s) B=0.
A list of the following objects:
... |
List(s) of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the working parameters of the data streams |
beta |
List of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the working parameters of the transition probabilities |
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m CIbeta(m)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m CIbeta(m)
This function can be used to convert angular covariates (e.g., ocean currents, wind direction) measured in radians relative to the x-axis to turning angle
covariates sutiable for circular-circular regression in fitHMM
or MIfitHMM
.
circAngles(refAngle, data, coordNames = c("x", "y"))
circAngles(refAngle, data, coordNames = c("x", "y"))
refAngle |
Numeric vector of standard direction angles (in radians) relative to the x-axis, where 0 = east, pi/2 = north, pi = west, -pi/2 = south |
data |
data frame containing fields for the x- and y-coordinates (identified by |
coordNames |
Names of the columns of coordinates in |
A vector of turning angles between the movement direction at time step t-1 and refAngle
at time t
# extract data from momentuHMM example data<-example$m$data # generate fake angle covariates u <- rnorm(nrow(data)) # horizontal component v <- rnorm(nrow(data)) # vertical component refAngle <- atan2(v,u) # add turning angle covariate to data data$cov3 <- circAngles(refAngle=refAngle,data=data)
# extract data from momentuHMM example data<-example$m$data # generate fake angle covariates u <- rnorm(nrow(data)) # horizontal component v <- rnorm(nrow(data)) # vertical component refAngle <- atan2(v,u) # add turning angle covariate to data data$cov3 <- circAngles(refAngle=refAngle,data=data)
Computes the standard errors and confidence intervals on the real (i.e., natural) scale of the data stream probability distribution parameters, as well as for the transition probabilities parameters. If covariates are included in the probability distributions or TPM formula, the mean values of non-factor covariates are used for calculating the natural parameters. For any covariate(s) of class 'factor', then the value(s) from the first observation in the data are used.
CIreal(m, alpha = 0.95, covs = NULL, parms = NULL) ## Default S3 method: CIreal(m, alpha = 0.95, covs = NULL, parms = NULL) ## S3 method for class 'hierarchical' CIreal(m, alpha = 0.95, covs = NULL, parms = NULL)
CIreal(m, alpha = 0.95, covs = NULL, parms = NULL) ## Default S3 method: CIreal(m, alpha = 0.95, covs = NULL, parms = NULL) ## S3 method for class 'hierarchical' CIreal(m, alpha = 0.95, covs = NULL, parms = NULL)
m |
A |
alpha |
Significance level of the confidence intervals. Default: 0.95 (i.e. 95% CIs). |
covs |
Data frame consisting of a single row indicating the covariate values to be used in the calculations. By default, no covariates are specified. |
parms |
Optional character vector indicating which groups of real parameters to calculate confidence intervals for (e.g., 'step', 'angle', 'gamma', 'delta', etc.). Default: NULL, in which case confidence intervals are calculated for all groups of parameters in the model. |
For any covariates that are not specified using covs
, the means of the covariate(s) are used
(unless the covariate is a factor, in which case the first factor in the data is used).
A list of the following objects:
... |
List(s) of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the natural parameters of the data streams |
gamma |
List of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the transition probabilities |
delta |
List of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the initial state probabilities |
hierGamma |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m ci1<-CIreal(m) # specify 'covs' ci2<-CIreal(m,covs=data.frame(cov1=mean(m$data$cov1),cov2=mean(m$data$cov2))) all.equal(ci1,ci2)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m ci1<-CIreal(m) # specify 'covs' ci2<-CIreal(m,covs=data.frame(cov1=mean(m$data$cov1),cov2=mean(m$data$cov2))) all.equal(ci1,ci2)
Transform a raster into a (x,y,z) list
collapseRaster(rast)
collapseRaster(rast)
rast |
|
List of three elements: x (grid of x values), y (grid of y values), and z (matrix of values).
This function can be used to merge crwData
or crwHierData
objects (as returned by crawlWrap
) with additional data streams
and/or covariates that are unrelated to location.
crawlMerge(crwData, data, Time.name)
crawlMerge(crwData, data, Time.name)
crwData |
A |
data |
A data frame containing required columns |
Time.name |
Character string indicating name of the time column to be used for merging |
Specifically, the function merges the crwData$crwPredict
data frame with data
based on the ID
, Time.name
, and, if crwData
is hierarchical, level
columns. Thus data
must contain ID
, Time.name
, and, if crwData
is hierarchical, level
columns.
Only rows of data
with ID
, Time.name
, and, if crwData
is hierarchical, level
values that exactly match crwData$crwPredict
are merged. Typically, the Time.name
column
in data
should match predicted times of locations in crwData$crwPredict
(i.e. those corresponding to crwData$crwPredict$locType=="p"
)
A crwData
object
## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crwMLE models to obsData and predict locations # at default intervals for both individuals crwOut <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model,attempts=100) # create data frame with fake data stream data <- data.frame(ID=rep(factor(c(1,2)),times=c(753,652)), time=c(1:753,1:652), fake=rpois(753+652,5)) # merge fake data stream with crwOut crwOut <- crawlMerge(crwOut,data,"time") ## End(Not run)
## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crwMLE models to obsData and predict locations # at default intervals for both individuals crwOut <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model,attempts=100) # create data frame with fake data stream data <- data.frame(ID=rep(factor(c(1,2)),times=c(753,652)), time=c(1:753,1:652), fake=rpois(753+652,5)) # merge fake data stream with crwOut crwOut <- crawlMerge(crwOut,data,"time") ## End(Not run)
Wrapper function for fitting crawl::crwMLE models and predicting locations with crawl::crwPredict for multiple individuals.
crawlWrap( obsData, timeStep = 1, ncores = 1, retryFits = 0, retrySD = 1, retryParallel = FALSE, mov.model = ~1, err.model = NULL, activity = NULL, drift = NULL, coord = c("x", "y"), proj = NULL, Time.name = "time", time.scale = "hours", theta, fixPar, method = "L-BFGS-B", control = NULL, constr = NULL, prior = NULL, need.hess = TRUE, initialSANN = list(maxit = 200), attempts = 1, predTime = NULL, fillCols = FALSE, coordLevel = NULL, ... )
crawlWrap( obsData, timeStep = 1, ncores = 1, retryFits = 0, retrySD = 1, retryParallel = FALSE, mov.model = ~1, err.model = NULL, activity = NULL, drift = NULL, coord = c("x", "y"), proj = NULL, Time.name = "time", time.scale = "hours", theta, fixPar, method = "L-BFGS-B", control = NULL, constr = NULL, prior = NULL, need.hess = TRUE, initialSANN = list(maxit = 200), attempts = 1, predTime = NULL, fillCols = FALSE, coordLevel = NULL, ... )
obsData |
data.frame object containing fields for animal ID ('ID'), time of observation (identified by |
timeStep |
Length of the time step at which to predict regular locations from the fitted model. Unless |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
retryFits |
Number of times to attempt to achieve convergence and valid (i.e., not NaN) variance estimates after the initial model fit. |
retrySD |
An optional list of scalars or vectors for each individual indicating the standard deviation to use for normal perturbations of |
retryParallel |
Logical indicating whether or not to perform |
mov.model |
List of mov.model objects (see |
err.model |
List of err.model objects (see |
activity |
List of activity objects (see |
drift |
List of drift objects (see |
coord |
A 2-vector of character values giving the names of the "x" and
"y" coordinates in |
proj |
A list of valid epsg integer codes or proj4string for |
Time.name |
Character indicating name of the location time column. See |
time.scale |
character. Scale for conversion of POSIX time to numeric for modeling. Defaults to "hours". |
theta |
List of theta objects (see |
fixPar |
List of fixPar objects (see |
method |
Optimization method that is passed to |
control |
Control list which is passed to |
constr |
List of constr objects (see |
prior |
List of prior objects (see |
need.hess |
A logical value which decides whether or not to evaluate the Hessian for parameter standard errors |
initialSANN |
Control list for |
attempts |
The number of times likelihood optimization will be
attempted in cases where the fit does not converge or is otherwise non-valid. Note this is not the same as |
predTime |
List of predTime objects (see |
fillCols |
Logical indicating whether or not to use the crawl:: |
coordLevel |
Character string indicating the level of the hierarchy for the location data. Ignored unless |
... |
Additional arguments that are ignored. |
Consult crwMLE
and crwPredict
for futher details about model fitting and prediction.
Note that the names of the list elements corresponding to each individual in mov.model
, err.model
, activity
, drift
,
theta
, fixPar
, constr
, prior
, and predTime
must match the individual IDs in obsData
. If only one element is provided
for any of these arguments, then the same element will be applied to all individuals.
A crwData
or crwHierData
object, i.e. a list of:
crwFits |
A list of |
crwPredict |
A |
The crwData
object is used in MIfitHMM
analyses that account for temporal irregularity or location measurement error.
## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crwMLE models to obsData and predict locations # at default intervals for both individuals crwOut1 <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model,attempts=100) # Fit the same crwMLE models and predict locations # at same intervals but specify for each individual using lists crwOut2 <- crawlWrap(obsData=obsData, theta=list(c(4,0),c(4,0)), fixPar=list(c(1,1,NA,NA),c(1,1,NA,NA)), err.model=list(err.model,err.model), predTime=list('1'=seq(1,633),'2'=seq(1,686))) ## End(Not run)
## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crwMLE models to obsData and predict locations # at default intervals for both individuals crwOut1 <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model,attempts=100) # Fit the same crwMLE models and predict locations # at same intervals but specify for each individual using lists crwOut2 <- crawlWrap(obsData=obsData, theta=list(c(4,0),c(4,0)), fixPar=list(c(1,1,NA,NA),c(1,1,NA,NA)), err.model=list(err.model,err.model), predTime=list('1'=seq(1,633),'2'=seq(1,686))) ## End(Not run)
crwData
objectsConstructor of crwData
objects
crwData(m)
crwData(m)
m |
A list of attributes of crawl output: |
An object crwData
.
crwHierData
objectsConstructor of crwHierData
objects
crwHierData(m)
crwHierData(m)
m |
A list of attributes of crawl output: |
An object crwHierData
.
crwHierSim
objectsConstructor of crwHierSim
objects
crwHierSim(m)
crwHierSim(m)
m |
A list of attributes required for multiple imputation data generated from a
|
An object crwHierSim
.
crwSim
objectsConstructor of crwSim
objects
crwSim(m)
crwSim(m)
m |
A list of attributes required for multiple imputation data generated from a
|
An object crwSim
.
ctds
fitted model objectThis is essentially a wrapper function for get.rate.matrix
and get.UD
that has been extended for momentuHMM
.
ctdsUD( ctds, spatialCovs, spatialCovs.grad, zero.idx = integer(), method = "lu", maxiter, start, tol, log = FALSE )
ctdsUD( ctds, spatialCovs, spatialCovs.grad, zero.idx = integer(), method = "lu", maxiter, start, tol, log = FALSE )
ctds |
A |
spatialCovs |
List of |
spatialCovs.grad |
List of |
zero.idx |
Integer vector of the indices of raster cells that are not passable and should be excluded. These are cells where movement should be impossible. Default is zero.idx=integer(). |
method |
Either "lu" (default) or "limit". See |
maxiter |
Total number of iterations for limit method if tolerance not reached first. Defaults to 100. Ignored for method='lu'. |
start |
A value for the starting distribution for the 'limit' method. Defaults to 1/num. cells. Ignored for method='lu'. |
tol |
Value used to assess convergence for limit method. If max(abs(pi1-pi0))<tol, limit method has converged. Defaults to sqrt(.Machine$double.eps). |
log |
Logical indicating whether or not to return the log UD. Default: |
A list of raster
objects.
Probability density function of the Bernoulli distribution (written in C++)
dbern_rcpp(x, prob, foo)
dbern_rcpp(x, prob, foo)
x |
Vector of quantiles |
prob |
success probability |
foo |
Unused (for compatibility with template) |
Vector of densities
Probability density function of the beta distribution (written in C++)
dbeta_rcpp(x, shape1, shape2)
dbeta_rcpp(x, shape1, shape2)
x |
Vector of quantiles |
shape1 |
Shape1 |
shape2 |
Shape2 |
Vector of densities
Probability density function of the categorical distribution (written in C++)
dcat_rcpp(x, prob, foo)
dcat_rcpp(x, prob, foo)
x |
Vector of quantiles |
prob |
success probability |
foo |
Unused (for compatibility with template) |
Vector of densities
Probability density function of Rice distribution under the CTCRW model of Johnson et al. (2008) (written in C++)
dcrwrice_rcpp(x, beta, sigma)
dcrwrice_rcpp(x, beta, sigma)
x |
numeric data vector of length |
beta |
correlation parameter |
sigma |
speed parameter |
Vector of densities
Probability density function of the Von Mises distribution under the CTCRW model of Johnson et al. (2008) (written in C++)
dcrwvm_rcpp(x, beta, sigma)
dcrwvm_rcpp(x, beta, sigma)
x |
numeric data vector of length |
beta |
correlation parameter |
sigma |
speed parameter |
Vector of densities
Probability density function of the exponential distribution (written in C++)
dexp_rcpp(x, rate, foo)
dexp_rcpp(x, rate, foo)
x |
Vector of quantiles |
rate |
Rate |
foo |
Unused (for compatibility with template) |
Vector of densities
Probability density function of the gamma distribution (written in C++)
dgamma_rcpp(x, mu, sigma)
dgamma_rcpp(x, mu, sigma)
x |
Vector of quantiles |
mu |
Mean |
sigma |
Standard deviation |
Vector of densities
Calculate distance between points y and z and turning angle between points x, y, and z
distAngle(x, y, z, type = "UTM", angleCov = TRUE)
distAngle(x, y, z, type = "UTM", angleCov = TRUE)
x |
location 1 |
y |
location 2 |
z |
location 3 |
type |
|
angleCov |
logical indicating to not return NA when x=y or y=z. Default: TRUE (i.e. NA is not returned if x=y or y=z). |
Used in prepData
and simData
to get distance and turning angle covariates
between locations (x1,x2), (y1,y2) and activity center (z1,z2).
If type='LL'
then distance is calculated as great circle distance using spDistsN1
, and turning angle is calculated based on initial bearings using bearing
.
2-vector with first element the distance between y and z and second element the turning angle between (x,y) and (y,z).
Probability density function of the log-normal distribution (written in C++)
dlnorm_rcpp(x, meanlog, sdlog)
dlnorm_rcpp(x, meanlog, sdlog)
x |
Vector of quantiles |
meanlog |
Mean of the distribution on the log-scale |
sdlog |
Standard deviation of the distribution on the log-scale |
Vector of densities
Probability density function of the logistic distribution (written in C++)
dlogis_rcpp(x, location, scale)
dlogis_rcpp(x, location, scale)
x |
Vector of quantiles |
location |
mean of the distribution |
scale |
Dispersion parameter |
Vector of densities
C++ implementation of multivariate Normal probability density function for multiple inputs
dmvnorm_rcpp(x, mean, varcovM)
dmvnorm_rcpp(x, mean, varcovM)
x |
data matrix of dimension |
mean |
mean vectors matrix of dimension |
varcovM |
(nbStates * (p + (p * p - p)/2 )) x n matrix of all sds (sd.x.state1, sd.x.state2, ..., sd.y.state1, ...) and then correlations (corr.xy.state1, corr.xy.state2, ..., corr.xz.state1, corr.xz.state2, ...) from lower triangle of variance-covariance matrix. |
matrix of densities of dimension K x n
.
Probability density function of the negative binomial distribution (written in C++)
dnbinom_rcpp(x, mu, size)
dnbinom_rcpp(x, mu, size)
x |
Vector of quantiles |
mu |
Mean of the distribution |
size |
Dispersion parameter |
Vector of densities
Probability density function of the normal distribution (written in C++)
dnorm_rcpp(x, mean, sd)
dnorm_rcpp(x, mean, sd)
x |
Vector of quantiles |
mean |
Mean of the distribution |
sd |
Standard deviation of the distribution |
Vector of densities
Probability density function of the Poisson distribution (written in C++)
dpois_rcpp(x, rate, foo)
dpois_rcpp(x, rate, foo)
x |
Vector of quantiles |
rate |
Rate |
foo |
Unused (for compatibility with template) |
Vector of densities
Probability density function of non-central student t (written in C++)
dt_rcpp(x, df, ncp)
dt_rcpp(x, df, ncp)
x |
Vector of quantiles |
df |
degrees of freedom |
ncp |
non-centrality parameter |
Vector of densities
Probability density function of the Von Mises distribution, defined as a function of the modified Bessel function of order 0 (written in C++)
dvm_rcpp(x, mu, kappa)
dvm_rcpp(x, mu, kappa)
x |
Vector of quantiles |
mu |
Mean |
kappa |
Concentration |
Vector of densities
Probability density function of the Weibull distribution (written in C++)
dweibull_rcpp(x, shape, scale)
dweibull_rcpp(x, shape, scale)
x |
Vector of quantiles |
shape |
Shape |
scale |
Scale |
Vector of densities
Probability density function of the wrapped Cauchy distribution (written in C++)
dwrpcauchy_rcpp(x, mu, rho)
dwrpcauchy_rcpp(x, mu, rho)
x |
Vector of quantiles |
mu |
Mean |
rho |
Concentration |
Vector of densities
These data are used in the examples and tests of functions to keep them as short as possible.
example miExample forest
example miExample forest
example
is a list of the following objects for demonstrating fitHMM
:
m
A momentuHMM
object
simPar
The parameters used to simulate data
par0
The initial parameters in the optimization to fit m
miExample
is a list of the following objects for demonstrating crawlWrap
, MIfitHMM
, and MIpool
:
obsData
Simulated observation data with measurement error and temporal irregularity (generated by simData
)
bPar
initial parameter estimates for MIfitHMM
examples
forest
is a simulated spatial covariate raster
object of the RasterLayer
class
Expand vector of free working parameters to vector of all working parameters including any fixed parameters (used in fitHMM.R and nLogLike.R)
expandPar( optPar, optInd, fixPar, wparIndex, betaCons, deltaCons, nbStates, nbCovsDelta, stationary, nbCovs, nbRecovs = 0, mixtures = 1, nbCovsPi = 0, TMB = FALSE, Par0, beta0, delta0 )
expandPar( optPar, optInd, fixPar, wparIndex, betaCons, deltaCons, nbStates, nbCovsDelta, stationary, nbCovs, nbRecovs = 0, mixtures = 1, nbCovsPi = 0, TMB = FALSE, Par0, beta0, delta0 )
optPar |
vector of free working parameters |
optInd |
indices of constrained parameters |
fixPar |
Vector of working parameters which are assumed known prior to fitting the model (NA indicates parameters is to be estimated) |
wparIndex |
Vector of indices for the elements of |
betaCons |
Matrix of the same dimension as |
deltaCons |
Matrix of the same dimension as |
nbStates |
Number of states of the HMM |
nbCovsDelta |
Number of initial distribution covariates |
stationary |
|
nbCovs |
Number of t.p.m. covariates |
nbRecovs |
Number of recharge covariates |
mixtures |
Number of mixtures for the state transition probabilities |
nbCovsPi |
Number of mixture probability covariates |
TMB |
Logical indicating whether or not optMethod='TMB'. Default: |
Par0 |
list of initial data stream parameter values. Ignored unless |
beta0 |
matrix of initial state transition parameter values. Ignored unless |
delta0 |
initial distribution parameter values. Ignored unless |
A vector of all working parameters including any fixed parameters
## Not run: nbStates <- 2 stepDist <- "gamma" # step distribution angleDist <- "vm" # turning angle distribution # extract data from momentuHMM example data <- example$m$data ### 1. fit the model to the simulated data # define initial values for the parameters mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar <- c(mu0,sigma0) # no zero-inflation, so no zero-mass included anglePar <- kappa0 # not estimating angle mean, so not included formula <- ~cov1+cos(cov2) # constrain cov1 effect to state 1 -> 2 and cov2 effect to state 2 -> 1 fixPar <- list(beta=c(NA,NA,0,NA,0,NA)) m <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=anglePar),formula=formula,fixPar=fixPar) # convert free parameter vector (m$mod$wpar) to full set of working parameters (m$mod$estimate) est <- momentuHMM:::expandPar(m$mod$wpar,m$conditions$optInd,unlist(m$conditions$fixPar), m$conditions$wparIndex,m$conditions$betaCons,m$conditions$deltaCons, nbStates, ncol(m$covsDelta)-1,m$conditions$stationary,nrow(m$mle$beta)-1) all(est==m$mod$estimate) ## End(Not run)
## Not run: nbStates <- 2 stepDist <- "gamma" # step distribution angleDist <- "vm" # turning angle distribution # extract data from momentuHMM example data <- example$m$data ### 1. fit the model to the simulated data # define initial values for the parameters mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar <- c(mu0,sigma0) # no zero-inflation, so no zero-mass included anglePar <- kappa0 # not estimating angle mean, so not included formula <- ~cov1+cos(cov2) # constrain cov1 effect to state 1 -> 2 and cov2 effect to state 2 -> 1 fixPar <- list(beta=c(NA,NA,0,NA,0,NA)) m <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=anglePar),formula=formula,fixPar=fixPar) # convert free parameter vector (m$mod$wpar) to full set of working parameters (m$mod$estimate) est <- momentuHMM:::expandPar(m$mod$wpar,m$conditions$optInd,unlist(m$conditions$fixPar), m$conditions$wparIndex,m$conditions$betaCons,m$conditions$deltaCons, nbStates, ncol(m$covsDelta)-1,m$conditions$stationary,nrow(m$mle$beta)-1) all(est==m$mod$estimate) ## End(Not run)
This function computes the exponential of a square matrix A
, defined as the sum from r=0
to infinity of A^r/r!
, using the default method of expm
.
expmatrix_rcpp(x, kappa = NA_real_, check = FALSE)
expmatrix_rcpp(x, kappa = NA_real_, check = FALSE)
x |
a square matrix. |
kappa |
maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is |
check |
logical indicating whether or not to check transition probability matrix for issues |
The matrix exponential of x.
Fit an approximate continuous-time (multivariate) hidden Markov model to the data provided, using numerical optimization of the log-likelihood function. The discrete-time approximation of the continuous-time model improves as the time between observations decreases. Note that any time-varying covariates are assumed piece-wise constant between observations.
fitCTHMM(data, ...) ## S3 method for class 'momentuHMMData' fitCTHMM( data, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, ... ) ## S3 method for class 'momentuHierHMMData' fitCTHMM( data, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, ... )
fitCTHMM(data, ...) ## S3 method for class 'momentuHMMData' fitCTHMM( data, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, ... ) ## S3 method for class 'momentuHierHMMData' fitCTHMM( data, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, ... )
data |
A |
... |
further arguments passed to or from other methods |
nbStates |
Number of states of the HMM. |
dist |
A named list indicating the probability distributions of the data streams. Currently
supported distributions are 'bern', 'beta', 'cat', 'ctds', 'exp', 'gamma', 'lnorm', 'logis', 'negbinom', 'norm', 'mvnorm2' (bivariate normal distribution), 'mvnorm3' (trivariate normal distribution),
'pois', 'rw_norm' (normal random walk), 'rw_mvnorm2' (bivariate normal random walk), 'rw_mvnorm3' (trivariate normal random walk), 't', and 'weibull'. See For continuous-time HMMs, the (multivariate) normal random walk ('rw_norm', 'rw_mvnorm2', 'rw_mvnorm3') and Poisson ('pois') distributions are modeled as a function of the time interval between observations |
Par0 |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in If |
beta0 |
Initial matrix of regression coefficients for the transition probabilities (more
information in 'Details').
Default: |
delta0 |
Initial value for the initial distribution of the HMM. Default: |
estAngleMean |
An optional named list indicating whether or not to estimate the angle mean for data streams with angular
distributions ('vm' and 'wrpcauchy'). For example, |
circularAngleMean |
An optional named list indicating whether to use circular-linear (FALSE) or circular-circular (TRUE)
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. For example,
Alternatively, |
formula |
Regression formula for the transition probability covariates. Default: |
formulaDelta |
Regression formula for the initial distribution. Default for |
stationary |
|
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. Default: |
nlmPar |
List of parameters to pass to the optimization function |
fit |
|
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of Design matrices specified using formulas allow standard functions in R formulas
(e.g., |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability
distribution parameters for each data stream. For each matrix, the first column pertains to the lower bound and the second column the upper bound. For example, for a 2-state model using the wrapped Cauchy ('wrpcauchy') distribution for
a data stream named 'angle' with |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
betaCons |
Matrix of the same dimension as |
deltaCons |
Matrix of the same dimension as |
mvnCoords |
Character string indicating the name of location data that are to be modeled using a multivariate normal distribution. For example, if |
stateNames |
Optional character vector of length nbStates indicating state names. |
knownStates |
Vector of values of the state process which are known prior to fitting the model (if any). Default: NULL (states are not known). This should be a vector with length the number of rows of 'data'; each element should either be an integer (the value of the known states) or NA if the state is not known. |
fixPar |
An optional list of vectors indicating parameters which are assumed known prior to fitting the model. Default: NULL
(no parameters are fixed). For data streams, each element of If |
retryFits |
Non-negative integer indicating the number of times to attempt to iteratively fit the model using random perturbations of the current parameter estimates as the
initial values for likelihood optimization. Normal(0, |
retrySD |
An optional list of scalars or vectors indicating the standard deviation to use for normal perturbations of each working scale parameter when |
ncores |
Number of cores to use for parallel processing when |
optMethod |
The optimization method to be used. Can be “nlm” (the default; see |
control |
A list of control parameters to be passed to |
prior |
A function that returns the log-density of the working scale parameter prior distribution(s). See 'Details'. |
modelName |
An optional character string providing a name for the fitted model. If provided, |
kappa |
maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). Default: |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
fitCTHMM
assumes the snapshot property applies to all data stream distributions (i.e. observations are “instantaneous”) except for the (multivariate) normal random walk (rw_norm
, rw_mvnorm2
, rw_mvnorm3
) and Poisson (pois
) distributions. For these particular distributions, the observed data are not “instantaneous”; they depend on the time interval between observations and, hence, the state sequence during the entire interval.
When fitting with
fitCTHMM
(or MIfitCTHMM
), it is critical that the frequency of observations is high relative to the serial correlation in the hidden state process in order for the discrete-time approximation of fitCTHMM
to be reasonably accurate for these distributions.
fitCTHMM.momentuHierHMMData
is very similar to fitCTHMM.momentuHMMData
except that instead of simply specifying the number of states (nbStates
), distributions (dist
), and a single t.p.m. formula (formula
), the hierStates
argument specifies the hierarchical nature of the states,
the hierDist
argument specifies the hierarchical nature of the data streams, and the hierFormula
argument specifies a t.p.m. formula for each level of the hierarchy. All are specified as
Node
objects from the data.tree
package.
A momentuHMM
or momentuHierHMM
object, i.e. a list of:
mle |
A named list of the maximum likelihood estimates of the parameters of the model (if the numerical algorithm
has indeed identified the global maximum of the likelihood function). Elements are included for the parameters of each
data strea, as well as |
CIreal |
Standard errors and 95% confidence intervals on the real (i.e., natural) scale of parameters |
CIbeta |
Standard errors and 95% confidence intervals on the beta (i.e., working) scale of parameters |
data |
The momentuHMMData or momentuHierHMMData object |
mod |
List object returned by the numerical optimizer |
conditions |
Conditions used to fit the model, e.g., |
rawCovs |
Raw covariate values for transition probabilities, as found in the data (if any). Used in |
stateNames |
The names of the states. |
knownStates |
Vector of values of the state process which are known. |
covsDelta |
Design matrix for initial distribution. |
Fit a (multivariate) hidden Markov model to the data provided, using numerical optimization of the log-likelihood function.
fitHMM(data, ...) ## S3 method for class 'momentuHMMData' fitHMM( data, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, betaRef = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, ... ) ## S3 method for class 'momentuHierHMMData' fitHMM( data, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, ... )
fitHMM(data, ...) ## S3 method for class 'momentuHMMData' fitHMM( data, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, betaRef = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, ... ) ## S3 method for class 'momentuHierHMMData' fitHMM( data, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = list(), fit = TRUE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, ncores = 1, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, ... )
data |
A |
... |
further arguments passed to or from other methods |
nbStates |
Number of states of the HMM. |
dist |
A named list indicating the probability distributions of the data streams. Currently
supported distributions are 'bern', 'beta', 'cat', exp', 'gamma', 'lnorm', 'logis', 'negbinom', 'norm', 'mvnorm2' (bivariate normal distribution), 'mvnorm3' (trivariate normal distribution),
'pois', 'rw_norm' (normal random walk), 'rw_mvnorm2' (bivariate normal random walk), 'rw_mvnorm3' (trivariate normal random walk), 'vm', 'vmConsensus', 'weibull', and 'wrpcauchy'. For example,
|
Par0 |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in If |
beta0 |
Initial matrix of regression coefficients for the transition probabilities (more
information in 'Details').
Default: |
delta0 |
Initial value for the initial distribution of the HMM. Default: |
estAngleMean |
An optional named list indicating whether or not to estimate the angle mean for data streams with angular
distributions ('vm' and 'wrpcauchy'). For example, |
circularAngleMean |
An optional named list indicating whether to use circular-linear (FALSE) or circular-circular (TRUE)
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. For example,
Alternatively, |
formula |
Regression formula for the transition probability covariates. Default: |
formulaDelta |
Regression formula for the initial distribution. Default for |
stationary |
|
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. Default: |
nlmPar |
List of parameters to pass to the optimization function |
fit |
|
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of Design matrices specified using formulas allow standard functions in R formulas
(e.g., |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability
distribution parameters for each data stream. For each matrix, the first column pertains to the lower bound and the second column the upper bound. For example, for a 2-state model using the wrapped Cauchy ('wrpcauchy') distribution for
a data stream named 'angle' with |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
betaCons |
Matrix of the same dimension as |
betaRef |
Numeric vector of length |
deltaCons |
Matrix of the same dimension as |
mvnCoords |
Character string indicating the name of location data that are to be modeled using a multivariate normal distribution. For example, if |
stateNames |
Optional character vector of length nbStates indicating state names. |
knownStates |
Vector of values of the state process which are known prior to fitting the model (if any). Default: NULL (states are not known). This should be a vector with length the number of rows of 'data'; each element should either be an integer (the value of the known states) or NA if the state is not known. |
fixPar |
An optional list of vectors indicating parameters which are assumed known prior to fitting the model. Default: NULL
(no parameters are fixed). For data streams, each element of If |
retryFits |
Non-negative integer indicating the number of times to attempt to iteratively fit the model using random perturbations of the current parameter estimates as the
initial values for likelihood optimization. Normal(0, |
retrySD |
An optional list of scalars or vectors indicating the standard deviation to use for normal perturbations of each working scale parameter when |
ncores |
Number of cores to use for parallel processing when |
optMethod |
The optimization method to be used. Can be “nlm” (the default; see |
control |
A list of control parameters to be passed to |
prior |
Either a function that returns the log-density of the working scale parameter prior distribution(s) (if |
modelName |
An optional character string providing a name for the fitted model. If provided, |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). Default: |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
By default the matrix beta0
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 formula
covariates, the matrix 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), beta0
has one row, for the intercept. While the diagonal elements are by default the reference elements, other elements can serve as the reference
using the betaRef
argument. For example, in a 3-state model, setting betaRef=c(3,2,3)
changes the reference elements to state transition 1 -> 3 for state 1 (instead of 1 -> 1), state
transition 2 -> 2 for state 2 (same as default), and state transition 3 -> 3 for state 3 (same as default).
When covariates are not included in formulaDelta
(i.e. formulaDelta=NULL
), then delta0
(and fixPar$delta
) are specified as a vector of length nbStates
that
sums to 1. When any formula is specified for formulaDelta
(e.g. formulaDelta=~1
, formulaDelta=~cov1
), then delta0
(and fixPar$delta
) must be specified
as a k x (nbStates
-1) matrix of working parameters, where k is the number of regression coefficients and the columns correspond to states 2:nbStates
. For example, in a 3-state
HMM with formulaDelta=~cov1+cov2
, the matrix delta0
has three rows (intercept + two covariates)
and 2 columns (corresponding to states 2 and 3). The initial distribution working parameters are transformed to the real scale as exp(covsDelta*Delta)/rowSums(exp(covsDelta*Delta))
, where covsDelta
is the N x k design matrix, Delta=cbind(rep(0,k),delta0)
is a k x nbStates
matrix of working parameters,
and N=length(unique(data$ID))
.
The choice of initial parameters (particularly Par0
and beta0
) is crucial to fit a model. The algorithm might not find
the global optimum of the likelihood function if the initial parameters are poorly chosen.
If DM
is specified for a particular data stream, then the initial values are specified on
the working (i.e., beta) scale of the parameters. The working scale of each parameter is determined by the link function used.
If a parameter P is bound by (0,Inf) then the working scale is the log(P) scale. If the parameter bounds are (-pi,pi) then the working
scale is tan(P/2) unless circular-circular regression is used. Otherwise if the parameter bounds are finite then logit(P) is the working scale. However, when both
zero- and one-inflation are included, then a multinomial logit link is used because the sum of the zeromass and onemass probability parameters cannot exceed 1.
The function getParDM
is intended to help with obtaining initial values on the working scale when specifying a design matrix and other
parameter constraints (see example below). When circular-circular regression is specified using circularAngleMean
, the working scale
for the mean turning angle is not as easily interpretable, but the
link function is atan2(sin(X)*B,1+cos(X)*B), where X are the angle covariates and B the angle coefficients (see Duchesne et al. 2015).
Under this formulation, the reference turning angle is 0 (i.e., movement in the same direction as the previous time step).
In other words, the mean turning angle is zero when the coefficient(s) B=0.
Circular-circular regression in momentuHMM
is designed for turning angles (not bearings) as computed by simData
and prepData
.
Any circular-circular regression angle covariates for time step t should therefore be relative to the previous
direction of movement for time step t-1. In other words, circular-circular regression covariates for time step t should be the turning angle
between the direction of movement for time step t-1 and the standard direction of the covariate relative to the x-axis for time step t. If provided standard directions in radians relative to the x-axis
(where 0 = east, pi/2 = north, pi = west, and -pi/2 = south), circAngles
or prepData
can perform this calculation for you.
When the circular-circular regression model is used, the special function angleFormula(cov,strength,by)
can be used in DM
for the mean of angular distributions (i.e. 'vm', 'vmConsensus', and 'wrpcauchy'),
where cov
is an angle covariate (e.g. wind direction), strength
is an optional positive real covariate (e.g. wind speed), and by
is an optional factor variable for individual- or group-level effects (e.g. ID, sex). The strength
argument allows angle covariates to be weighted based on their relative strength or importance at time step t as in
Rivest et al. (2016). In this case, the link function for the mean angle is atan2((Z * sin(X)) %*% B,1+(Z * cos(X)) %*% B), where X are the angle covariates, Z the strength covariates, and B the angle coefficients (see Rivest et al. 2016).
State-specific formulas can be specified in DM
using special formula functions. These special functions can take
the names paste0("state",1:nbStates)
(where the integer indicates the state-specific formula). For example,
DM=list(step=list(mean=~cov1+state1(cov2),sd=~cov2+state2(cov1)))
includes cov1
on the mean parameter for all states, cov2
on the mean parameter for state 1, cov2
on the sd parameter for all states, and cov1
on the sd parameter for state 2.
State- and parameter-specific formulas can be specified for transition probabilities in formula
using special formula functions.
These special functions can take the names paste0("state",1:nbStates)
(where the integer indicates the current state from which transitions occur),
paste0("toState",1:nbStates)
(where the integer indicates the state to which transitions occur),
or paste0("betaCol",nbStates*(nbStates-1))
(where the integer indicates the column of the beta
matrix). For example with nbStates=3
,
formula=~cov1+betaCol1(cov2)+state3(cov3)+toState1(cov4)
includes cov1
on all transition probability parameters, cov2
on the beta
column corresponding
to the transition from state 1->2, cov3
on transition probabilities from state 3 (i.e., beta
columns corresponding to state transitions 3->1 and 3->2),
and cov4
on transition probabilities to state 1 (i.e., beta
columns corresponding to state transitions 2->1 and 3->1).
betaCons
can be used to impose equality constraints among the t.p.m. parameters. It must be a matrix of the same dimension as beta0
and be composed of integers, where each beta parameter is sequentially indexed in a column-wise fashion (see checkPar0
). Parameter indices in betaCons
must therefore be integers between 1
and nbStates*(nbStates-1)
.
Use of betaCons
is perhaps best demonstrated by example. If no constraints are imposed (the default), then betaCons=matrix(1:length(beta0),nrow(beta0),ncol(beta0))
such that
each beta parameter is (column-wise) sequentially identified by a unique integer. Suppose we wish to fit a model with nbStates=3
states and a covariate (‘cov1’) on the t.p.m. With no constraints on the t.p.m., we would have
betaCons=matrix(1:(2*(nbStates*(nbStates-1))),nrow=2,ncol=nbStates*(nbStates-1),dimnames=list(c("(Intercept)","cov1"),c("1 -> 2","1 -> 3","2 -> 1","2 -> 3","3 -> 1","3 -> 2")))
. If we then wanted to constrain the t.p.m. such that the covariate effect is identical for transitions from state 1 to states 2 and 3 (and vice versa), we have
betaCons=matrix(c(1,2,3,2,5,6,7,8,9,6,11,12),nrow=2,ncol=nbStates*(nbStates-1),dimnames=list(c("(Intercept)","cov1"),c("1 -> 2","1 -> 3","2 -> 1","2 -> 3","3 -> 1","3 -> 2")))
; this results in 10 estimated beta parameters (instead of 12), the “cov1” effects indexed by a “2” (“1 -> 2” and “1 -> 3”) constrained to be equal, and
the “cov1” effects indexed by a “6” (“2 -> 1” and “3 -> 1”) constrained to be equal.
Now suppose we instead wish to constrain these sets of state transition probabilities to be equal, i.e., Pr(1 -> 2) = Pr(1 -> 3) and Pr(2 -> 1) = Pr(3 -> 1); then we have betaCons=matrix(c(1,2,1,2,5,6,7,8,5,6,11,12),nrow=2,ncol=nbStates*(nbStates-1),dimnames=list(c("(Intercept)","cov1"),c("1 -> 2","1 -> 3","2 -> 1","2 -> 3","3 -> 1","3 -> 2")))
Cyclical relationships (e.g., hourly, monthly) may be modeled in DM
or formula
using the cosinor(x,period)
special formula function for covariate x
and sine curve period of time length period
. For example, if
the data are hourly, a 24-hour cycle can be modeled using ~cosinor(cov1,24)
, where the covariate cov1
is a repeating sequential series
of integers indicating the hour of day (0,1,...,23,0,1,...,23,0,1,...
) (note that fitHMM
will not do this for you, the appropriate covariate must be included in data
; see example below).
The cosinor(x,period)
function converts x
to 2 covariates cosinorCos(x)=cos(2*pi*x/period)
and cosinorSin(x)=sin(2*pi*x/period
for inclusion in the model (i.e., 2 additional parameters per state). The amplitude of the sine wave
is thus sqrt(B_cos^2 + B_sin^2)
, where B_cos
and B_sin
are the working parameters correponding to cosinorCos(x)
and cosinorSin(x)
, respectively (e.g., see Cornelissen 2014).
Similar to that used in crawlWrap
, the prior
argument is a user-specified function that returns the log-density of the working scale parameter prior distribution(s). In addition to including prior information about parameters, one area where priors can be particularly useful is for handling numerical issues that can arise when parameters are near a boundary.
When parameters are near boundaries, they can wander into the “nether regions” of the parameter space during optimization. For example, setting prior=function(par) {sum(dnorm(par,0,sd,log=TRUE))}
with a reasonably large sd
(e.g. 100 or 1000) can help prevent working parameters
from straying too far along the real line. Here par
is the vector of working scale parameters (as returned by fitHMM
, e.g., see example$m$mod$estimate
) in the following order: data stream working parameters (in order names(dist)
), beta working parameters, and delta working parameters. Instead of specifying the same prior on all parameters, different priors could be specified on different parameters (and not all parameters must have user-specified priors). For example,
prior=function(par){dnorm(par[3],0,100,log=TRUE)}
would only specify a prior for the third working parameter. Note that the prior
function must return a scalar on the log scale. See 'harbourSealExample.R' in the “vignettes” source directory for an example using the prior
argument.
fitHMM.momentuHierHMMData
is very similar to fitHMM.momentuHMMData
except that instead of simply specifying the number of states (nbStates
), distributions (dist
), and a single t.p.m. formula (formula
), the hierStates
argument specifies the hierarchical nature of the states,
the hierDist
argument specifies the hierarchical nature of the data streams, and the hierFormula
argument specifies a t.p.m. formula for each level of the hierarchy. All are specified as
Node
objects from the data.tree
package.
A momentuHMM
or momentuHierHMM
object, i.e. a list of:
mle |
A named list of the maximum likelihood estimates of the parameters of the model (if the numerical algorithm
has indeed identified the global maximum of the likelihood function). Elements are included for the parameters of each
data strea, as well as |
CIreal |
Standard errors and 95% confidence intervals on the real (i.e., natural) scale of parameters |
CIbeta |
Standard errors and 95% confidence intervals on the beta (i.e., working) scale of parameters |
data |
The momentuHMMData or momentuHierHMMData object |
mod |
List object returned by the numerical optimizer |
conditions |
Conditions used to fit the model, e.g., |
rawCovs |
Raw covariate values for transition probabilities, as found in the data (if any). Used in |
stateNames |
The names of the states. |
knownStates |
Vector of values of the state process which are known. |
covsDelta |
Design matrix for initial distribution. |
Cornelissen, G. 2014. Cosinor-based rhythmometry. Theoretical Biology and Medical Modelling 11:16.
Duchesne, T., Fortin, D., Rivest L-P. 2015. Equivalence between step selection functions and biased correlated random walks for statistical inference on animal movement. PLoS ONE 10 (4): e0122947.
Langrock R., King R., Matthiopoulos J., Thomas L., Fortin D., Morales J.M. 2012. Flexible and practical modeling of animal telemetry data: hidden Markov models and extensions. Ecology, 93 (11), 2336-2342.
Leos-Barajas, V., Gangloff, E.J., Adam, T., Langrock, R., van Beest, F.M., Nabe-Nielsen, J. and Morales, J.M. 2017. Multi-scale modeling of animal movement and general behavior data using hidden Markov models with hierarchical structures. Journal of Agricultural, Biological and Environmental Statistics, 22 (3), 232-248.
Maruotti, A., and T. Ryden. 2009. A semiparametric approach to hidden Markov models under longitudinal observations. Statistics and Computing 19: 381-393.
McClintock B.T., King R., Thomas L., Matthiopoulos J., McConnell B.J., Morales J.M. 2012. A general discrete-time modeling framework for animal movement using multistate random walks. Ecological Monographs, 82 (3), 335-349.
McClintock B.T., Russell D.J., Matthiopoulos J., King R. 2013. Combining individual animal movement and ancillary biotelemetry data to investigate population-level activity budgets. Ecology, 94 (4), 838-849.
Patterson T.A., Basson M., Bravington M.V., Gunn J.S. 2009. Classifying movement behaviour in relation to environmental conditions using hidden Markov models. Journal of Animal Ecology, 78 (6), 1113-1123.
Rivest, LP, Duchesne, T, Nicosia, A, Fortin, D, 2016. A general angular regression model for the analysis of data on animal movement in ecology. Journal of the Royal Statistical Society: Series C (Applied Statistics), 65(3):445-463.
nbStates <- 2 stepDist <- "gamma" # step distribution angleDist <- "vm" # turning angle distribution # extract data from momentuHMM example data <- example$m$data ### 1. fit the model to the simulated data # define initial values for the parameters mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar <- c(mu0,sigma0) # no zero-inflation, so no zero-mass included anglePar <- kappa0 # not estimating angle mean, so not included formula <- ~cov1+cos(cov2) m <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=anglePar),formula=formula) print(m) ## Not run: ### 2. fit the exact same model to the simulated data using DM formulas # Get initial values for the parameters on working scale Par0 <- getParDM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=list(mean=~1,sd=~1),angle=list(concentration=~1))) mDMf <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0,formula=formula, DM=list(step=list(mean=~1,sd=~1),angle=list(concentration=~1))) print(mDMf) ### 3. fit the exact same model to the simulated data using DM matrices # define DM DMm <- list(step=diag(4),angle=diag(2)) # user-specified dimnames not required but are recommended dimnames(DMm$step) <- list(c("mean_1","mean_2","sd_1","sd_2"), c("mean_1:(Intercept)","mean_2:(Intercept)", "sd_1:(Intercept)","sd_2:(Intercept)")) dimnames(DMm$angle) <- list(c("concentration_1","concentration_2"), c("concentration_1:(Intercept)","concentration_2:(Intercept)")) mDMm <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0,formula=formula, DM=DMm) print(mDMm) ### 4. fit step mean parameter covariate model to the simulated data using DM stepDMf <- list(mean=~cov1,sd=~1) Par0 <- getParDM(data,nbStates,list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=stepDMf,angle=DMm$angle)) mDMfcov <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0, formula=formula, DM=list(step=stepDMf,angle=DMm$angle)) print(mDMfcov) ### 5. fit the exact same step mean parameter covariate model using DM matrix stepDMm <- matrix(c(1,0,0,0,"cov1",0,0,0,0,1,0,0,0,"cov1",0,0, 0,0,1,0,0,0,0,1),4,6,dimnames=list(c("mean_1","mean_2","sd_1","sd_2"), c("mean_1:(Intercept)","mean_1:cov1","mean_2:(Intercept)","mean_2:cov1", "sd_1:(Intercept)","sd_2:(Intercept)"))) Par0 <- getParDM(data,nbStates,list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=stepDMm,angle=DMm$angle)) mDMmcov <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0, formula=formula, DM=list(step=stepDMm,angle=DMm$angle)) print(mDMmcov) ### 6. fit circular-circular angle mean covariate model to the simulated data using DM # Generate fake circular covariate using circAngles data$cov3 <- circAngles(refAngle=2*atan(rnorm(nrow(data))),data) # Fit circular-circular regression model for angle mean # Note no intercepts are estimated for angle means because these are by default # the previous movement direction (i.e., a turning angle of zero) mDMcircf <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=c(0,0,Par0$angle)), formula=formula, estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE), DM=list(angle=list(mean=~cov3,concentration=~1))) print(mDMcircf) ### 7. fit the exact same circular-circular angle mean model using DM matrices # Note no intercept terms are included in DM for angle means because the intercept is # by default the previous movement direction (i.e., a turning angle of zero) mDMcircm <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=c(0,0,Par0$angle)), formula=formula, estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE), DM=list(angle=matrix(c("cov3",0,0,0,0,"cov3",0,0,0,0,1,0,0,0,0,1),4,4))) print(mDMcircm) ### 8. Cosinor and state-dependent formulas nbStates<-2 dist<-list(step="gamma") Par<-list(step=c(100,1000,50,100)) # include 24-hour cycle on all transition probabilities # include 12-hour cycle on transitions from state 2 formula=~cosinor(hour24,24)+state2(cosinor(hour12,12)) # specify appropriate covariates covs<-data.frame(hour24=0:23,hour12=0:11) beta<-matrix(c(-1.5,1,1,NA,NA,-1.5,-1,-1,1,1),5,2) # row names for beta not required but can be helpful rownames(beta)<-c("(Intercept)", "cosinorCos(hour24, 24)", "cosinorSin(hour24, 24)", "cosinorCos(hour12, 12)", "cosinorSin(hour12, 12)") data.cos<-simData(nbStates=nbStates,dist=dist,Par=Par, beta=beta,formula=formula,covs=covs) m.cosinor<-fitHMM(data.cos,nbStates=nbStates,dist=dist,Par0=Par,formula=formula) m.cosinor ### 9. Piecewise constant B-spline on step length mean and angle concentration nObs <- 1000 # length of simulated track cov <- data.frame(time=1:nObs) # time covariate for splines dist <- list(step="gamma",angle="vm") stepDM <- list(mean=~splines2::bSpline(time,df=2,degree=0),sd=~1) angleDM <- list(mean=~1,concentration=~splines2::bSpline(time,df=2,degree=0)) DM <- list(step=stepDM,angle=angleDM) Par <- list(step=c(log(1000),1,-1,log(100)),angle=c(0,log(10),2,-5)) data.spline<-simData(obsPerAnimal=nObs,nbStates=1,dist=dist,Par=Par,DM=DM,covs=cov) Par0 <- list(step=Par$step,angle=Par$angle[-1]) m.spline<-fitHMM(data.spline,nbStates=1,dist=dist,Par0=Par0, DM=list(step=stepDM, angle=angleDM["concentration"])) ### 10. Initial state (delta) based on covariate nObs <- 100 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.01,0.75)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs))) # sex covariate formulaDelta <- ~ sex + 0 # Female begins in state 1, male begins in state 2 delta <- matrix(c(-100,100),2,1,dimnames=list(c("sexF","sexM"),"state 2")) data.delta<-simData(nbAnimals=2,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, delta=delta,formulaDelta=formulaDelta,covs=cov) Par0 <- list(step=Par$step, angle=Par$angle[3:4]) m.delta <- fitHMM(data.delta, nbStates=2, dist=dist, Par0 = Par0, formulaDelta=formulaDelta) ### 11. Two mixtures based on covariate nObs <- 100 nbAnimals <- 20 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.1,2)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs*nbAnimals/2))) formulaPi <- ~ sex + 0 # Females more likely in mixture 1, males more likely in mixture 2 beta <- list(beta=matrix(c(-1.5,-0.5,-1.5,-3),2,2), pi=matrix(c(-2,2),2,1,dimnames=list(c("sexF","sexM"),"mix2"))) data.mix<-simData(nbAnimals=nbAnimals,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, beta=beta,formulaPi=formulaPi,mixtures=2,covs=cov) Par0 <- list(step=Par$step, angle=Par$angle[3:4]) m.mix <- fitHMM(data.mix, nbStates=2, dist=dist, Par0 = Par0, beta0=beta,formulaPi=formulaPi,mixtures=2) ## End(Not run)
nbStates <- 2 stepDist <- "gamma" # step distribution angleDist <- "vm" # turning angle distribution # extract data from momentuHMM example data <- example$m$data ### 1. fit the model to the simulated data # define initial values for the parameters mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar <- c(mu0,sigma0) # no zero-inflation, so no zero-mass included anglePar <- kappa0 # not estimating angle mean, so not included formula <- ~cov1+cos(cov2) m <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=anglePar),formula=formula) print(m) ## Not run: ### 2. fit the exact same model to the simulated data using DM formulas # Get initial values for the parameters on working scale Par0 <- getParDM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=list(mean=~1,sd=~1),angle=list(concentration=~1))) mDMf <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0,formula=formula, DM=list(step=list(mean=~1,sd=~1),angle=list(concentration=~1))) print(mDMf) ### 3. fit the exact same model to the simulated data using DM matrices # define DM DMm <- list(step=diag(4),angle=diag(2)) # user-specified dimnames not required but are recommended dimnames(DMm$step) <- list(c("mean_1","mean_2","sd_1","sd_2"), c("mean_1:(Intercept)","mean_2:(Intercept)", "sd_1:(Intercept)","sd_2:(Intercept)")) dimnames(DMm$angle) <- list(c("concentration_1","concentration_2"), c("concentration_1:(Intercept)","concentration_2:(Intercept)")) mDMm <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0,formula=formula, DM=DMm) print(mDMm) ### 4. fit step mean parameter covariate model to the simulated data using DM stepDMf <- list(mean=~cov1,sd=~1) Par0 <- getParDM(data,nbStates,list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=stepDMf,angle=DMm$angle)) mDMfcov <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0, formula=formula, DM=list(step=stepDMf,angle=DMm$angle)) print(mDMfcov) ### 5. fit the exact same step mean parameter covariate model using DM matrix stepDMm <- matrix(c(1,0,0,0,"cov1",0,0,0,0,1,0,0,0,"cov1",0,0, 0,0,1,0,0,0,0,1),4,6,dimnames=list(c("mean_1","mean_2","sd_1","sd_2"), c("mean_1:(Intercept)","mean_1:cov1","mean_2:(Intercept)","mean_2:cov1", "sd_1:(Intercept)","sd_2:(Intercept)"))) Par0 <- getParDM(data,nbStates,list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=stepDMm,angle=DMm$angle)) mDMmcov <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=Par0, formula=formula, DM=list(step=stepDMm,angle=DMm$angle)) print(mDMmcov) ### 6. fit circular-circular angle mean covariate model to the simulated data using DM # Generate fake circular covariate using circAngles data$cov3 <- circAngles(refAngle=2*atan(rnorm(nrow(data))),data) # Fit circular-circular regression model for angle mean # Note no intercepts are estimated for angle means because these are by default # the previous movement direction (i.e., a turning angle of zero) mDMcircf <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=c(0,0,Par0$angle)), formula=formula, estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE), DM=list(angle=list(mean=~cov3,concentration=~1))) print(mDMcircf) ### 7. fit the exact same circular-circular angle mean model using DM matrices # Note no intercept terms are included in DM for angle means because the intercept is # by default the previous movement direction (i.e., a turning angle of zero) mDMcircm <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar,angle=c(0,0,Par0$angle)), formula=formula, estAngleMean=list(angle=TRUE), circularAngleMean=list(angle=TRUE), DM=list(angle=matrix(c("cov3",0,0,0,0,"cov3",0,0,0,0,1,0,0,0,0,1),4,4))) print(mDMcircm) ### 8. Cosinor and state-dependent formulas nbStates<-2 dist<-list(step="gamma") Par<-list(step=c(100,1000,50,100)) # include 24-hour cycle on all transition probabilities # include 12-hour cycle on transitions from state 2 formula=~cosinor(hour24,24)+state2(cosinor(hour12,12)) # specify appropriate covariates covs<-data.frame(hour24=0:23,hour12=0:11) beta<-matrix(c(-1.5,1,1,NA,NA,-1.5,-1,-1,1,1),5,2) # row names for beta not required but can be helpful rownames(beta)<-c("(Intercept)", "cosinorCos(hour24, 24)", "cosinorSin(hour24, 24)", "cosinorCos(hour12, 12)", "cosinorSin(hour12, 12)") data.cos<-simData(nbStates=nbStates,dist=dist,Par=Par, beta=beta,formula=formula,covs=covs) m.cosinor<-fitHMM(data.cos,nbStates=nbStates,dist=dist,Par0=Par,formula=formula) m.cosinor ### 9. Piecewise constant B-spline on step length mean and angle concentration nObs <- 1000 # length of simulated track cov <- data.frame(time=1:nObs) # time covariate for splines dist <- list(step="gamma",angle="vm") stepDM <- list(mean=~splines2::bSpline(time,df=2,degree=0),sd=~1) angleDM <- list(mean=~1,concentration=~splines2::bSpline(time,df=2,degree=0)) DM <- list(step=stepDM,angle=angleDM) Par <- list(step=c(log(1000),1,-1,log(100)),angle=c(0,log(10),2,-5)) data.spline<-simData(obsPerAnimal=nObs,nbStates=1,dist=dist,Par=Par,DM=DM,covs=cov) Par0 <- list(step=Par$step,angle=Par$angle[-1]) m.spline<-fitHMM(data.spline,nbStates=1,dist=dist,Par0=Par0, DM=list(step=stepDM, angle=angleDM["concentration"])) ### 10. Initial state (delta) based on covariate nObs <- 100 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.01,0.75)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs))) # sex covariate formulaDelta <- ~ sex + 0 # Female begins in state 1, male begins in state 2 delta <- matrix(c(-100,100),2,1,dimnames=list(c("sexF","sexM"),"state 2")) data.delta<-simData(nbAnimals=2,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, delta=delta,formulaDelta=formulaDelta,covs=cov) Par0 <- list(step=Par$step, angle=Par$angle[3:4]) m.delta <- fitHMM(data.delta, nbStates=2, dist=dist, Par0 = Par0, formulaDelta=formulaDelta) ### 11. Two mixtures based on covariate nObs <- 100 nbAnimals <- 20 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.1,2)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs*nbAnimals/2))) formulaPi <- ~ sex + 0 # Females more likely in mixture 1, males more likely in mixture 2 beta <- list(beta=matrix(c(-1.5,-0.5,-1.5,-3),2,2), pi=matrix(c(-2,2),2,1,dimnames=list(c("sexF","sexM"),"mix2"))) data.mix<-simData(nbAnimals=nbAnimals,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, beta=beta,formulaPi=formulaPi,mixtures=2,covs=cov) Par0 <- list(step=Par$step, angle=Par$angle[3:4]) m.mix <- fitHMM(data.mix, nbStates=2, dist=dist, Par0 = Par0, beta0=beta,formulaPi=formulaPi,mixtures=2) ## End(Not run)
Convert hierarchical HMM structure to a conventional HMM
formatHierHMM( data, hierStates, hierDist, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, workBounds = NULL, betaCons = NULL, deltaCons = NULL, fixPar = NULL, checkData = TRUE, ... )
formatHierHMM( data, hierStates, hierDist, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, workBounds = NULL, betaCons = NULL, deltaCons = NULL, fixPar = NULL, checkData = TRUE, ... )
data |
|
hierStates |
A hierarchical data structure |
hierDist |
A hierarchical data structure |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). See |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). See |
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). See |
workBounds |
A list with elements named |
betaCons |
A hierarchical data structure |
deltaCons |
A hierarchical data structure |
fixPar |
A list with elements named |
checkData |
logical indicating whether or not to check the suitability of |
... |
further arguments passed to or from other methods |
A list of arguments needed for specifying a hierarchical HMM as a conventional HMM in fitHMM
or MIfitHMM
, including:
nbStates |
See |
dist |
See |
formula |
See |
formulaDelta |
See |
beta0 |
See |
delta0 |
See |
betaRef |
See |
betaCons |
See |
deltaCons |
See |
fixPar |
See |
workBounds |
See |
stateNames |
See |
Get names of any covariates used in probability distribution parameters
getCovNames(m, p, distname)
getCovNames(m, p, distname)
m |
|
p |
list returned by |
distname |
Name of the data stream |
A list of:
DMterms |
Names of all covariates included in the design matrix for the data stream |
DMpartems |
A list of the names of all covariates for each of the probability distribution parameters |
Loop for creating full design matrix (X) from pseudo-design matrix (DM). Written in C++. Used in getDM
.
getDM_rcpp(X, covs, DM, nr, nc, cov, nbObs)
getDM_rcpp(X, covs, DM, nr, nc, cov, nbObs)
X |
full design matrix |
covs |
matrix of covariates |
DM |
pseudo design matrix |
nr |
number of rows in design matrix |
nc |
number of column in design matrix |
cov |
covariate names |
nbObs |
number of observations |
full design matrix (X)
Calculate gradient of spatial covariates using bilinear interpolation
getGradients(data, spatialCovs, collapseRast, coordNames = c("x", "y"))
getGradients(data, spatialCovs, collapseRast, coordNames = c("x", "y"))
data |
Data frame of data streams. At a minimum, it must contain fields matching |
spatialCovs |
List of |
collapseRast |
List of collapsed |
coordNames |
Names of the coordinates in |
The gradients are appended to data
with “.x
” (easting gradient) and “.y
” (northing gradient) suffixes added to the names of spatialCovs
. For example, if cov1
is the name of a spatial covariate, then the returned data
object will include the fields “cov1.x
” and “cov1.y
”.
Get starting values from momentuHMM, miHMM, or miSum object returned by fitHMM, MIfitHMM, or MIpool
getPar(m)
getPar(m)
m |
A |
A list of parameter values (Par, beta, delta) that can be used as starting values in fitHMM
or MIfitHMM
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m Par <- getPar(m)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m Par <- getPar(m)
momentuHMM
or momentuHierHMM
model fitFor nested models, this function will extract starting parameter values (i.e., Par0
in fitHMM
or MIfitHMM
) from an existing momentuHMM
or momentuHierHMM
model fit based on the provided arguments for the new model.
Any parameters that are not in common between model
and the new model (as specified by the arguments) are set to 0
. This function is intended to help users incrementally build and fit more complicated models from simpler nested models (and vice versa).
getPar0(model, ...) ## Default S3 method: getPar0( model, nbStates = length(model$stateNames), estAngleMean = model$conditions$estAngleMean, circularAngleMean = model$conditions$circularAngleMean, formula = model$conditions$formula, formulaDelta = model$conditions$formulaDelta, stationary = model$conditions$stationary, mixtures = model$conditions$mixtures, formulaPi = model$conditions$formulaPi, DM = model$conditions$DM, betaRef = model$conditions$betaRef, stateNames = model$stateNames, ... ) ## S3 method for class 'hierarchical' getPar0( model, hierStates = model$conditions$hierStates, estAngleMean = model$conditions$estAngleMean, circularAngleMean = model$conditions$circularAngleMean, hierFormula = model$conditions$hierFormula, hierFormulaDelta = model$conditions$hierFormulaDelta, mixtures = model$conditions$mixtures, formulaPi = model$conditions$formulaPi, DM = model$conditions$DM, ... )
getPar0(model, ...) ## Default S3 method: getPar0( model, nbStates = length(model$stateNames), estAngleMean = model$conditions$estAngleMean, circularAngleMean = model$conditions$circularAngleMean, formula = model$conditions$formula, formulaDelta = model$conditions$formulaDelta, stationary = model$conditions$stationary, mixtures = model$conditions$mixtures, formulaPi = model$conditions$formulaPi, DM = model$conditions$DM, betaRef = model$conditions$betaRef, stateNames = model$stateNames, ... ) ## S3 method for class 'hierarchical' getPar0( model, hierStates = model$conditions$hierStates, estAngleMean = model$conditions$estAngleMean, circularAngleMean = model$conditions$circularAngleMean, hierFormula = model$conditions$hierFormula, hierFormulaDelta = model$conditions$hierFormulaDelta, mixtures = model$conditions$mixtures, formulaPi = model$conditions$formulaPi, DM = model$conditions$DM, ... )
model |
A |
... |
further arguments passed to or from other methods |
nbStates |
Number of states in the new model. Default: |
estAngleMean |
Named list indicating whether or not the angle mean for data streams with angular
distributions ('vm' and 'wrpcauchy') are to be estimated in the new model. Default: |
circularAngleMean |
Named list indicating whether circular-linear or circular-circular
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles are to be used in the new model. See |
formula |
Regression formula for the transition probability covariates of the new model (see |
formulaDelta |
Regression formula for the initial distribution covariates of the new model (see |
stationary |
|
mixtures |
Number of mixtures for the state transition probabilities (see |
formulaPi |
Regression formula for the mixture distribution probabilities (see |
DM |
Named list indicating the design matrices to be used for the probability distribution parameters of each data stream in the new model (see |
betaRef |
Numeric vector of length |
stateNames |
Character vector of length |
hierStates |
A hierarchical model structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy (see |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
All other fitHMM
(or MIfitHMM
) model specifications (e.g., dist
, hierDist
, userBounds
, workBounds
, etc.) and data
are assumed to be the same
for model
and the new model (as specified by nbStates
, hierStates
, estAngleMean
, circularAngleMean
, formula
, hierFormula
, formulaDelta
, hierFormulaDelta
, DM
, etc.).
Note that for hierarchical models, getPar0
will return hierarchical data.tree objects (hierBeta
and hierDelta
) with the default values for fixPar
, betaCons
, and deltaCons
;
if hierarchical t.p.m. or initial distribution parameters are subject to constraints, then these must be set manually on the list object returned by getPar0
.
A named list containing starting values suitable for Par0
and beta0
arguments in fitHMM
or MIfitHMM
:
Par |
A list of vectors of state-dependent probability distribution parameters for
each data stream specified in |
beta |
Matrix of regression coefficients for the transition probabilities |
delta |
Initial distribution of the HMM. Only returned if |
getPar
, getParDM
, fitHMM
, MIfitHMM
# model is a momentuHMM object, automatically loaded with the package model <- example$m data <- model$data dist <- model$conditions$dist nbStates <- length(model$stateNames) estAngleMean <- model$conditions$estAngleMean newformula <- ~cov1+cov2 Par0 <- getPar0(model,formula=newformula) ## Not run: newModel <- fitHMM(model$data,dist=dist,nbStates=nbStates, Par0=Par0$Par,beta0=Par0$beta, formula=newformula, estAngleMean=estAngleMean) ## End(Not run) newDM1 <- list(step=list(mean=~cov1,sd=~cov1)) Par0 <- getPar0(model,DM=newDM1) ## Not run: newModel1 <- fitHMM(model$data,dist=dist,nbStates=nbStates, Par0=Par0$Par,beta0=Par0$beta, formula=model$conditions$formula, estAngleMean=estAngleMean, DM=newDM1) ## End(Not run) # same model but specify DM for step using matrices newDM2 <- list(step=matrix(c(1,0,0,0, "cov1",0,0,0, 0,1,0,0, 0,"cov1",0,0, 0,0,1,0, 0,0,"cov1",0, 0,0,0,1, 0,0,0,"cov1"),nrow=nbStates*2)) # to be extracted, new design matrix column names must match # column names of model$conditions$fullDM colnames(newDM2$step)<-paste0(rep(c("mean_","sd_"),each=2*nbStates), rep(1:nbStates,each=2), rep(c(":(Intercept)",":cov1"),2*nbStates)) Par0 <- getPar0(model,DM=newDM2) ## Not run: newModel2 <- fitHMM(model$data,dist=dist,nbStates=nbStates, Par0=Par0$Par,beta0=Par0$beta, formula=model$conditions$formula, estAngleMean=estAngleMean, DM=newDM2) ## End(Not run)
# model is a momentuHMM object, automatically loaded with the package model <- example$m data <- model$data dist <- model$conditions$dist nbStates <- length(model$stateNames) estAngleMean <- model$conditions$estAngleMean newformula <- ~cov1+cov2 Par0 <- getPar0(model,formula=newformula) ## Not run: newModel <- fitHMM(model$data,dist=dist,nbStates=nbStates, Par0=Par0$Par,beta0=Par0$beta, formula=newformula, estAngleMean=estAngleMean) ## End(Not run) newDM1 <- list(step=list(mean=~cov1,sd=~cov1)) Par0 <- getPar0(model,DM=newDM1) ## Not run: newModel1 <- fitHMM(model$data,dist=dist,nbStates=nbStates, Par0=Par0$Par,beta0=Par0$beta, formula=model$conditions$formula, estAngleMean=estAngleMean, DM=newDM1) ## End(Not run) # same model but specify DM for step using matrices newDM2 <- list(step=matrix(c(1,0,0,0, "cov1",0,0,0, 0,1,0,0, 0,"cov1",0,0, 0,0,1,0, 0,0,"cov1",0, 0,0,0,1, 0,0,0,"cov1"),nrow=nbStates*2)) # to be extracted, new design matrix column names must match # column names of model$conditions$fullDM colnames(newDM2$step)<-paste0(rep(c("mean_","sd_"),each=2*nbStates), rep(1:nbStates,each=2), rep(c(":(Intercept)",":cov1"),2*nbStates)) Par0 <- getPar0(model,DM=newDM2) ## Not run: newModel2 <- fitHMM(model$data,dist=dist,nbStates=nbStates, Par0=Par0$Par,beta0=Par0$beta, formula=model$conditions$formula, estAngleMean=estAngleMean, DM=newDM2) ## End(Not run)
Convert starting values on the natural scale of data stream probability distributions to a feasible set of working scale parameters based on a design matrix and other parameter constraints.
getParDM(data, ...) ## Default S3 method: getParDM( data = data.frame(), nbStates, dist, Par, zeroInflation = NULL, oneInflation = NULL, estAngleMean = NULL, circularAngleMean = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, optMethod = "nlm", ... ) ## S3 method for class 'hierarchical' getParDM( data = data.frame(), hierStates, hierDist, Par, zeroInflation = NULL, oneInflation = NULL, estAngleMean = NULL, circularAngleMean = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, ... )
getParDM(data, ...) ## Default S3 method: getParDM( data = data.frame(), nbStates, dist, Par, zeroInflation = NULL, oneInflation = NULL, estAngleMean = NULL, circularAngleMean = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, optMethod = "nlm", ... ) ## S3 method for class 'hierarchical' getParDM( data = data.frame(), hierStates, hierDist, Par, zeroInflation = NULL, oneInflation = NULL, estAngleMean = NULL, circularAngleMean = NULL, DM = NULL, userBounds = NULL, workBounds = NULL, ... )
data |
Optional If a data frame is provided, then either |
... |
further arguments passed to or from other methods |
nbStates |
Number of states of the HMM. |
dist |
A named list indicating the probability distributions of the data streams. Currently
supported distributions are 'bern', 'beta', 'exp', 'gamma', 'lnorm', 'norm', 'mvnorm2' (bivariate normal distribution), 'mvnorm3' (trivariate normal distribution),
'pois', 'rw_norm' (normal random walk), 'rw_mvnorm2' (bivariate normal random walk), 'rw_mvnorm3' (trivariate normal random walk), 'vm', 'vmConsensus', 'weibull', and 'wrpcauchy'. For example,
|
Par |
A named list containing vectors of state-dependent probability distribution parameters for
each data stream specified in |
zeroInflation |
A named list of logicals indicating whether the probability distributions of the data streams should be zero-inflated. If |
oneInflation |
Named list of logicals indicating whether the probability distributions of the data streams are one-inflated. If |
estAngleMean |
An optional named list indicating whether or not to estimate the angle mean for data streams with angular
distributions ('vm' and 'wrpcauchy'). Any |
circularAngleMean |
An optional named list indicating whether to use circular-linear or circular-circular
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. See |
DM |
A named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability
distribution parameters for each data stream. For example, for a 2-state model using the wrapped Cauchy ('wrpcauchy') distribution for
a data stream named 'angle' with |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
optMethod |
The optimization method to be used. Can be “nlm” (the default; see |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
If design matrix includes non-factor covariates, then natural scale parameters are assumed to correspond to the
mean value(s) for the covariate(s) (if nrow(data)>1
) and getParDM
simply returns one possible solution to the
system of linear equations defined by Par
, DM
, and any other constraints using singular value decomposition.
This can be helpful for exploring relationships between the natural and working scale parameters when covariates are included, but getParDM
will not necessarily return “good” starting values (i.e., Par0
) for fitHMM
or MIfitHMM
.
A list of parameter values that can be used as starting values (Par0
) in fitHMM
or MIfitHMM
getPar
, getPar0
, fitHMM
, MIfitHMM
# data is a momentuHMMData object, automatically loaded with the package data <- example$m$data stepDist <- "gamma" angleDist <- "vm" nbStates <- 2 stepPar0 <- c(15,50,10,20) # natural scale mean_1, mean_2, sd_1, sd_2 anglePar0 <- c(0.7,1.5) # natural scale conentration_1, concentration_2 # get working parameters for 'DM' that constrains step length mean_1 < mean_2 stepDM <- matrix(c(1,1,0,0,0,1,0,0,0,0,1,0,0,0,0,1),4,4, dimnames=list(NULL,c("mean:(Intercept)","mean_2", "sd_1:(Intercept)","sd_2:(Intercept)"))) stepworkBounds <- matrix(c(-Inf,Inf),4,2,byrow=TRUE, dimnames=list(colnames(stepDM),c("lower","upper"))) stepworkBounds["mean_2","lower"] <- 0 #coefficient for 'mean_2' constrained to be positive wPar0 <- getParDM(nbStates=2,dist=list(step=stepDist), Par=list(step=stepPar0), DM=list(step=stepDM),workBounds=list(step=stepworkBounds)) ## Not run: # Fit HMM using wPar0 as initial values for the step data stream mPar <- fitHMM(data,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par0=list(step=wPar0$step,angle=anglePar0), DM=list(step=stepDM),workBounds=list(step=stepworkBounds)) ## End(Not run) # get working parameters for 'DM' using 'cov1' and 'cov2' covariates stepDM2 <- list(mean=~cov1,sd=~cov2) wPar20 <- getParDM(data,nbStates=2,dist=list(step=stepDist), Par=list(step=stepPar0), DM=list(step=stepDM2)) ## Not run: # Fit HMM using wPar20 as initial values for the step data stream mPar2 <- fitHMM(data,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par0=list(step=wPar20$step,angle=anglePar0), DM=list(step=stepDM2)) ## End(Not run)
# data is a momentuHMMData object, automatically loaded with the package data <- example$m$data stepDist <- "gamma" angleDist <- "vm" nbStates <- 2 stepPar0 <- c(15,50,10,20) # natural scale mean_1, mean_2, sd_1, sd_2 anglePar0 <- c(0.7,1.5) # natural scale conentration_1, concentration_2 # get working parameters for 'DM' that constrains step length mean_1 < mean_2 stepDM <- matrix(c(1,1,0,0,0,1,0,0,0,0,1,0,0,0,0,1),4,4, dimnames=list(NULL,c("mean:(Intercept)","mean_2", "sd_1:(Intercept)","sd_2:(Intercept)"))) stepworkBounds <- matrix(c(-Inf,Inf),4,2,byrow=TRUE, dimnames=list(colnames(stepDM),c("lower","upper"))) stepworkBounds["mean_2","lower"] <- 0 #coefficient for 'mean_2' constrained to be positive wPar0 <- getParDM(nbStates=2,dist=list(step=stepDist), Par=list(step=stepPar0), DM=list(step=stepDM),workBounds=list(step=stepworkBounds)) ## Not run: # Fit HMM using wPar0 as initial values for the step data stream mPar <- fitHMM(data,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par0=list(step=wPar0$step,angle=anglePar0), DM=list(step=stepDM),workBounds=list(step=stepworkBounds)) ## End(Not run) # get working parameters for 'DM' using 'cov1' and 'cov2' covariates stepDM2 <- list(mean=~cov1,sd=~cov2) wPar20 <- getParDM(data,nbStates=2,dist=list(step=stepDist), Par=list(step=stepPar0), DM=list(step=stepDM2)) ## Not run: # Fit HMM using wPar20 as initial values for the step data stream mPar2 <- fitHMM(data,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par0=list(step=wPar20$step,angle=anglePar0), DM=list(step=stepDM2)) ## End(Not run)
Computation of the transition probability matrix for each time step as a function of the covariates and the regression parameters.
getTrProbs(data, ...) ## Default S3 method: getTrProbs( data, nbStates, beta, workBounds = NULL, formula = ~1, mixtures = 1, betaRef = NULL, stateNames = NULL, getCI = FALSE, covIndex = NULL, alpha = 0.95, Time.name = NULL, Time.unit = "auto", rateMatrix = FALSE, kappa = Inf, ... ) ## S3 method for class 'hierarchical' getTrProbs( data, hierStates, hierBeta, workBounds = NULL, hierFormula = NULL, mixtures = 1, hierDist, getCI = FALSE, covIndex = NULL, alpha = 0.95, Time.name = NULL, Time.unit = "auto", rateMatrix = FALSE, kappa = Inf, ... )
getTrProbs(data, ...) ## Default S3 method: getTrProbs( data, nbStates, beta, workBounds = NULL, formula = ~1, mixtures = 1, betaRef = NULL, stateNames = NULL, getCI = FALSE, covIndex = NULL, alpha = 0.95, Time.name = NULL, Time.unit = "auto", rateMatrix = FALSE, kappa = Inf, ... ) ## S3 method for class 'hierarchical' getTrProbs( data, hierStates, hierBeta, workBounds = NULL, hierFormula = NULL, mixtures = 1, hierDist, getCI = FALSE, covIndex = NULL, alpha = 0.95, Time.name = NULL, Time.unit = "auto", rateMatrix = FALSE, kappa = Inf, ... )
data |
If a data frame is provided, then either |
... |
further arguments passed to or from other methods; most are ignored if |
nbStates |
Number of states. Ignored unless |
beta |
Matrix of regression coefficients for the transition probabilities |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the transition probability parameters ('beta' and, for recharge models, 'g0' and 'theta'). |
formula |
Regression formula for the transition probability covariates. Ignored unless |
mixtures |
Number of mixtures for the state transition probabilities. Ignored unless |
betaRef |
Indices of reference elements for t.p.m. multinomial logit link. Ignored unless |
stateNames |
Optional character vector of length nbStates indicating state names. Ignored unless |
getCI |
Logical indicating whether to calculate standard errors and logit-transformed confidence intervals based on fitted |
covIndex |
Integer vector indicating specific rows of the data to be used in the calculations. This can be useful for reducing unnecessarily long computation times (paricularly when |
alpha |
Significance level of the confidence intervals (if |
Time.name |
Character string indicating name of the time column in |
Time.unit |
Character string indicating units for time difference between observations (e.g. 'auto', 'secs', 'mins', 'hours', 'days', 'weeks'). Ignored unless |
rateMatrix |
Logical indicating whether to return the transition rate matrix. Default: |
kappa |
maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is |
hierStates |
A hierarchical model structure |
hierBeta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). See |
hierDist |
A hierarchical data structure |
If mixtures=1
, an array of dimension nbStates
x nbStates
x nrow(data)
containing the t.p.m for each observation in data
.
If mixtures>1
, a list of length mixtures
, where each element is an array of dimension nbStates
x nbStates
x nrow(data)
containing the t.p.m for each observation in data
.
If getCI=TRUE
then a list of arrays is returned (with elements est
, se
, lower
, and upper
).
If a hierarchical HMM structure is provided, then a hierarchical data structure containing the state transition probabilities for each time step at each level of the hierarchy ('gamma') is returned.
m <- example$m trProbs <- getTrProbs(m) # equivalent trProbs <- getTrProbs(m$data,nbStates=2,beta=m$mle$beta,formula=m$conditions$formula) ## Not run: # calculate SEs and 95% CIs trProbsSE <- getTrProbs(m, getCI=TRUE) # plot estimates and CIs for each state transition par(mfrow=c(2,2)) for(i in 1:2){ for(j in 1:2){ plot(trProbsSE$est[i,j,],type="l", ylim=c(0,1), ylab=paste(i,"->",j)) arrows(1:dim(trProbsSE$est)[3], trProbsSE$lower[i,j,], 1:dim(trProbsSE$est)[3], trProbsSE$upper[i,j,], length=0.025, angle=90, code=3, col=gray(.5), lwd=1.3) } } # limit calculations to first 10 observations trProbsSE_10 <- getTrProbs(m, getCI=TRUE, covIndex=1:10) ## End(Not run)
m <- example$m trProbs <- getTrProbs(m) # equivalent trProbs <- getTrProbs(m$data,nbStates=2,beta=m$mle$beta,formula=m$conditions$formula) ## Not run: # calculate SEs and 95% CIs trProbsSE <- getTrProbs(m, getCI=TRUE) # plot estimates and CIs for each state transition par(mfrow=c(2,2)) for(i in 1:2){ for(j in 1:2){ plot(trProbsSE$est[i,j,],type="l", ylim=c(0,1), ylab=paste(i,"->",j)) arrows(1:dim(trProbsSE$est)[3], trProbsSE$lower[i,j,], 1:dim(trProbsSE$est)[3], trProbsSE$upper[i,j,], length=0.025, angle=90, code=3, col=gray(.5), lwd=1.3) } } # limit calculations to first 10 observations trProbsSE_10 <- getTrProbs(m, getCI=TRUE, covIndex=1:10) ## End(Not run)
HMMfits
objectsConstructor of HMMfits
objects
HMMfits(m)
HMMfits(m)
m |
A list of
|
An object HMMfits
.
Check that an object is of class crwData
. Used in MIfitHMM
.
is.crwData(x)
is.crwData(x)
x |
An R object |
TRUE
if x
is of class crwData
, FALSE
otherwise.
Check that an object is of class crwHierData
. Used in MIfitHMM
.
is.crwHierData(x)
is.crwHierData(x)
x |
An R object |
TRUE
if x
is of class crwHierData
, FALSE
otherwise.
Check that an object is of class crwHierSim
.
is.crwHierSim(x)
is.crwHierSim(x)
x |
An R object |
TRUE
if x
is of class crwHierSim
, FALSE
otherwise.
Check that an object is of class HMMfits
.
is.HMMfits(x)
is.HMMfits(x)
x |
An R object |
TRUE
if x
is of class HMMfits
, FALSE
otherwise.
Check that an object is of class momentuHierHMM
. Used in CIreal
, CIbeta
,
plotPR
, plotStates
, pseudoRes
, stateProbs
,
and viterbi
.
is.momentuHierHMM(x)
is.momentuHierHMM(x)
x |
An R object |
TRUE
if x
is of class momentuHierHMM
, FALSE
otherwise.
Check that an object is of class momentuHierHMMData
. Used in fitHMM
.
is.momentuHierHMMData(x)
is.momentuHierHMMData(x)
x |
An R object |
TRUE
if x
is of class momentuHierHMMData
, FALSE
otherwise.
Check that an object is of class momentuHMM
. Used in CIreal
, CIbeta
,
plotPR
, plotStates
, pseudoRes
, stateProbs
,
and viterbi
.
is.momentuHMM(x)
is.momentuHMM(x)
x |
An R object |
TRUE
if x
is of class momentuHMM
, FALSE
otherwise.
Check that an object is of class momentuHMMData
. Used in fitHMM
.
is.momentuHMMData(x)
is.momentuHMMData(x)
x |
An R object |
TRUE
if x
is of class momentuHMMData
, FALSE
otherwise.
Used in stateProbs
and pseudoRes
.
logAlpha(m)
logAlpha(m)
m |
A |
A list of length model$conditions$mixtures
where each element is a matrix of forward log-probabilities for each mixture.
## Not run: # m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m la <- momentuHMM:::logAlpha(m) ## End(Not run)
## Not run: # m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m la <- momentuHMM:::logAlpha(m) ## End(Not run)
Used in stateProbs
.
logBeta(m)
logBeta(m)
m |
A |
A list of length model$conditions$mixtures
where each element is a matrix of backward log-probabilities for each mixture.
## Not run: # m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m lb <- momentuHMM:::logBeta(m) ## End(Not run)
## Not run: # m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m lb <- momentuHMM:::logBeta(m) ## End(Not run)
Fit an approximate continuous-time (multivariate) hidden Markov model to multiple imputation data. Multiple imputation is a method for accommodating missing data, temporal-irregularity, or location measurement error in hidden Markov models, where pooled parameter estimates reflect uncertainty attributable to observation error. The discrete-time approximation of the continuous-time model improves as the time between observations decreases. Note that any time-varying covariates are assumed piece-wise constant between observations.
MIfitCTHMM(miData, ...) ## Default S3 method: MIfitCTHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, Time.name = "time", Time.unit = "auto", CTDS = FALSE, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, rast, directions = 4, zero.idx = integer(), interpMethod = "ShortestPath", spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, ... ) ## S3 method for class 'hierarchical' MIfitCTHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, Time.name = "time", Time.unit = "auto", CTDS = FALSE, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, rast, directions = 4, zero.idx = integer(), interpMethod = "ShortestPath", spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, ... )
MIfitCTHMM(miData, ...) ## Default S3 method: MIfitCTHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, Time.name = "time", Time.unit = "auto", CTDS = FALSE, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, rast, directions = 4, zero.idx = integer(), interpMethod = "ShortestPath", spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, ... ) ## S3 method for class 'hierarchical' MIfitCTHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, Time.name = "time", Time.unit = "auto", CTDS = FALSE, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, kappa = Inf, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, rast, directions = 4, zero.idx = integer(), interpMethod = "ShortestPath", spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, ... )
miData |
A |
... |
further arguments passed to or from other methods |
nSims |
Number of imputations in which to fit the HMM using |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
poolEstimates |
Logical indicating whether or not to calculate pooled parameter estimates across the |
alpha |
Significance level for calculating confidence intervals of pooled estimates when |
na.rm |
Logical indicating whether or not to exclude model fits with |
Time.name |
Character string indicating name of the time column. Default: 'time'. Ignored if |
Time.unit |
Character string indicating units for time difference between observations (e.g. 'auto', 'secs', 'mins', 'hours', 'days', 'weeks'). Ignored unless |
CTDS |
Logical indicating whether multiple imputation location data should be prepared for a continuous-time, discrete-space analysis. See |
nbStates |
Number of states of the HMM. See |
dist |
A named list indicating the probability distributions of the data streams. See |
Par0 |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in |
beta0 |
Initial matrix of regression coefficients for the transition probabilities. See |
delta0 |
Initial values for the initial distribution of the HMM. See |
estAngleMean |
An optional named list indicating whether or not to estimate the angle mean for data streams with angular
distributions ('vm' and 'wrpcauchy'). See |
circularAngleMean |
An optional named list indicating whether to use circular-linear or circular-circular
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. See |
formula |
Regression formula for the transition probability covariates. See |
formulaDelta |
Regression formula for the initial distribution. See |
stationary |
|
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. See |
nlmPar |
List of parameters to pass to the optimization function |
fit |
|
useInitial |
Logical indicating whether or not to use parameter estimates for the first model fit as initial values for all subsequent model fits.
If |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. See |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability
distribution parameters for each data stream. See |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. See |
betaCons |
Matrix of the same dimension as |
deltaCons |
Matrix of the same dimension as |
mvnCoords |
Character string indicating the name of location data that are to be modeled using a multivariate normal distribution. For example, if |
stateNames |
Optional character vector of length nbStates indicating state names. |
knownStates |
Vector of values of the state process which are known prior to fitting the
model (if any). See |
fixPar |
An optional list of vectors indicating parameters which are assumed known prior to fitting the model. See |
retryFits |
Non-negative integer indicating the number of times to attempt to iteratively fit the model using random perturbations of the current parameter estimates as the
initial values for likelihood optimization. See |
retrySD |
An optional list of scalars or vectors indicating the standard deviation to use for normal perturbations of each working scale parameter when |
optMethod |
The optimization method to be used. Can be “nlm” (the default; see |
control |
A list of control parameters to be passed to |
prior |
A function that returns the log-density of the working scale parameter prior distribution(s). See |
modelName |
An optional character string providing a name for the fitted model. If provided, |
kappa |
maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is |
covNames |
Names of any covariates in |
spatialCovs |
List of raster layer(s) for any spatial covariates. See |
centers |
2-column matrix providing the x-coordinates (column 1) and y-coordinates (column 2) for any activity centers (e.g., potential
centers of attraction or repulsion) from which distance and angle covariates will be calculated based on realizations of the position process.
See |
centroids |
List where each element is a data frame containing the x-coordinates ('x'), y-coordinates ('y'), and times (with user-specified name, e.g., ‘time’) for centroids (i.e., dynamic activity centers where the coordinates can change over time)
from which distance and angle covariates will be calculated based on the location data. See |
angleCovs |
Character vector indicating the names of any circular-circular regression angular covariates in |
altCoordNames |
Character string indicating an alternative name for the returned location data. See |
gradient |
Logical indicating whether or not to calculate gradients of |
smoothGradient |
Logical indicating whether or not to calculate smoothed gradients. See |
weights |
Numeric vector indicating the weight for the points in a smoothed version of the gradient. See |
scaleGrad |
Numeric scalar indicating the distance from the observed location to each smoothing point (before adjusting by time step). See |
method |
Method for obtaining weights for movement parameter samples. See |
parIS |
Size of the parameter importance sample. See |
dfSim |
Degrees of freedom for the t approximation to the parameter posterior. See 'df' argument in |
grid.eps |
Grid size for |
crit |
Criterion for deciding "significance" of quadrature points
(difference in log-likelihood). See |
scaleSim |
Scale multiplier for the covariance matrix of the t approximation. See 'scale' argument in |
quad.ask |
Logical, for method='quadrature'. Whether or not the sampler should ask if quadrature sampling should take place. It is used to stop the sampling if the number of likelihood evaluations would be extreme. Default: FALSE. Ignored if |
force.quad |
A logical indicating whether or not to force the execution
of the quadrature method for large parameter vectors. See |
fullPost |
Logical indicating whether to draw parameter values as well to simulate full posterior. See |
dfPostIS |
Degrees of freedom for multivariate t distribution approximation to parameter posterior. See 'df' argument in |
scalePostIS |
Extra scaling factor for t distribution approximation. See 'scale' argument in |
thetaSamp |
If multiple parameter samples are available in |
export |
Character vector of the names of any additional objects or functions in the global environment that are used in |
rast |
A raster object or raster stack object that will define the discrete-space grid cells for the CTDS movement path. See |
directions |
Integer. See |
zero.idx |
Integer vector of the indices of raster cells that are not passable and should be excluded. See |
interpMethod |
Specifies interpolation method. See |
spatialCovs.grad |
List of |
crw |
Logical. See |
normalize.gradients |
Logical. See |
grad.point.decreasing |
Logical. See |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy. See |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
miData
can either be a crwData
or crwHierData
object (as returned by crawlWrap
), a crwSim
or crwHierSim
object (as returned by MIfitCTHMM
when fit=FALSE
),
or a list of momentuHMMData
or momentuHierHMMData
objects (e.g., each element of the list as returned by prepData
).
If miData
is a crwData
(or crwHierData
) object, MIfitCTHMM
uses a combination of
crwSimulator
, crwPostIS
, prepData
, and fitCTHMM
to draw nSims
realizations of the position process
and fit the specified HMM to each imputation of the data. The vast majority of MIfitCTHMM
arguments are identical to the corresponding arguments from these functions.
If miData
is a crwData
or crwHierData
object, nSims
determines both the number of realizations of the position process to draw
(using crwSimulator
and crwPostIS
) as well as the number of HMM fits.
If miData
is a crwSim
(or crwHierSim
) object or a list of momentuHMMData
(or momentuHierHMMData
) object(s), the specified HMM will simply be fitted to each of the momentuHMMData
(or momentuHierHMMData
) objects
and all arguments related to crwSimulator
, crwPostIS
, or prepData
are ignored.
fitCTHMM
assumes the snapshot property applies to all data stream distributions (i.e. observations are “instantaneous”) except for the (multivariate) normal random walk (rw_norm
, rw_mvnorm2
, rw_mvnorm3
) and Poisson (pois
) distributions. For these particular distributions, the observed data are not “instantaneous”; they depend on the time interval between observations and, hence, the state sequence during the entire interval.
When fitting with
fitCTHMM
(or MIfitCTHMM
), it is critical that the frequency of observations is high relative to the serial correlation in the hidden state process in order for the discrete-time approximation of fitCTHMM
to be reasonably accurate for these distributions.
If nSims>1
, poolEstimates=TRUE
, and fit=TRUE
, a miHMM
object, i.e., a list consisting of:
miSum |
|
HMMfits |
List of length |
If poolEstimates=FALSE
and fit=TRUE
, a list of length nSims
consisting of momentuHMM
objects is returned.
However, if fit=FALSE
and miData
is a crwData
object, then MIfitCTHMM
returns a crwSim
object, i.e., a list containing miData
(a list of
momentuHMMData
objects) and crwSimulator
(a list of crwSimulator
objects),and most other arguments related to fitCTHMM
are ignored.
Hooten M.B., Johnson D.S., McClintock B.T., Morales J.M. 2017. Animal Movement: Statistical Models for Telemetry Data. CRC Press, Boca Raton.
McClintock B.T. 2017. Incorporating telemetry error into hidden Markov movement models using multiple imputation. Journal of Agricultural, Biological, and Environmental Statistics.
crawlWrap
, crwPostIS
, crwSimulator
, fitCTHMM
, getParDM
, MIpool
, prepData
, prepCTDS
Fit a (multivariate) hidden Markov model to multiple imputation data. Multiple imputation is a method for accommodating missing data, temporal-irregularity, or location measurement error in hidden Markov models, where pooled parameter estimates reflect uncertainty attributable to observation error.
MIfitHMM(miData, ...) ## Default S3 method: MIfitHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, betaRef = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, ... ) ## S3 method for class 'hierarchical' MIfitHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, ... )
MIfitHMM(miData, ...) ## Default S3 method: MIfitHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, nbStates, dist, Par0, beta0 = NULL, delta0 = NULL, estAngleMean = NULL, circularAngleMean = NULL, formula = ~1, formulaDelta = NULL, stationary = FALSE, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, betaRef = NULL, deltaCons = NULL, mvnCoords = NULL, stateNames = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, ... ) ## S3 method for class 'hierarchical' MIfitHMM( miData, nSims, ncores = 1, poolEstimates = TRUE, alpha = 0.95, na.rm = FALSE, hierStates, hierDist, Par0, hierBeta = NULL, hierDelta = NULL, estAngleMean = NULL, circularAngleMean = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, nlmPar = NULL, fit = TRUE, useInitial = FALSE, DM = NULL, userBounds = NULL, workBounds = NULL, betaCons = NULL, deltaCons = NULL, mvnCoords = NULL, knownStates = NULL, fixPar = NULL, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), prior = NULL, modelName = NULL, covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, smoothGradient = FALSE, weights = c(1/2, 1/8, 1/8, 1/8, 1/8), scaleGrad = NULL, method = "IS", parIS = 1000, dfSim = Inf, grid.eps = 1, crit = 2.5, scaleSim = 1, quad.ask = FALSE, force.quad = TRUE, fullPost = TRUE, dfPostIS = Inf, scalePostIS = 1, thetaSamp = NULL, export = NULL, ... )
miData |
A |
... |
further arguments passed to or from other methods |
nSims |
Number of imputations in which to fit the HMM using |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
poolEstimates |
Logical indicating whether or not to calculate pooled parameter estimates across the |
alpha |
Significance level for calculating confidence intervals of pooled estimates when |
na.rm |
Logical indicating whether or not to exclude model fits with |
nbStates |
Number of states of the HMM. See |
dist |
A named list indicating the probability distributions of the data streams. See |
Par0 |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in |
beta0 |
Initial matrix of regression coefficients for the transition probabilities. See |
delta0 |
Initial values for the initial distribution of the HMM. See |
estAngleMean |
An optional named list indicating whether or not to estimate the angle mean for data streams with angular
distributions ('vm' and 'wrpcauchy'). See |
circularAngleMean |
An optional named list indicating whether to use circular-linear or circular-circular
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. See |
formula |
Regression formula for the transition probability covariates. See |
formulaDelta |
Regression formula for the initial distribution. See |
stationary |
|
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. See |
nlmPar |
List of parameters to pass to the optimization function |
fit |
|
useInitial |
Logical indicating whether or not to use parameter estimates for the first model fit as initial values for all subsequent model fits.
If |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. See |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability
distribution parameters for each data stream. See |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. See |
betaCons |
Matrix of the same dimension as |
betaRef |
Numeric vector of length |
deltaCons |
Matrix of the same dimension as |
mvnCoords |
Character string indicating the name of location data that are to be modeled using a multivariate normal distribution. For example, if |
stateNames |
Optional character vector of length nbStates indicating state names. |
knownStates |
Vector of values of the state process which are known prior to fitting the
model (if any). See |
fixPar |
An optional list of vectors indicating parameters which are assumed known prior to fitting the model. See |
retryFits |
Non-negative integer indicating the number of times to attempt to iteratively fit the model using random perturbations of the current parameter estimates as the
initial values for likelihood optimization. See |
retrySD |
An optional list of scalars or vectors indicating the standard deviation to use for normal perturbations of each working scale parameter when |
optMethod |
The optimization method to be used. Can be “nlm” (the default; see |
control |
A list of control parameters to be passed to |
prior |
Either a function that returns the log-density of the working scale parameter prior distribution(s) (if |
modelName |
An optional character string providing a name for the fitted model. If provided, |
covNames |
Names of any covariates in |
spatialCovs |
List of raster layer(s) for any spatial covariates. See |
centers |
2-column matrix providing the x-coordinates (column 1) and y-coordinates (column 2) for any activity centers (e.g., potential
centers of attraction or repulsion) from which distance and angle covariates will be calculated based on realizations of the position process.
See |
centroids |
List where each element is a data frame containing the x-coordinates ('x'), y-coordinates ('y'), and times (with user-specified name, e.g., ‘time’) for centroids (i.e., dynamic activity centers where the coordinates can change over time)
from which distance and angle covariates will be calculated based on the location data. See |
angleCovs |
Character vector indicating the names of any circular-circular regression angular covariates in |
altCoordNames |
Character string indicating an alternative name for the returned location data. See |
gradient |
Logical indicating whether or not to calculate gradients of |
smoothGradient |
Logical indicating whether or not to calculate smoothed gradients. See |
weights |
Numeric vector indicating the weight for the points in a smoothed version of the gradient. See |
scaleGrad |
Numeric scalar indicating the distance from the observed location to each smoothing point. See |
method |
Method for obtaining weights for movement parameter samples. See |
parIS |
Size of the parameter importance sample. See |
dfSim |
Degrees of freedom for the t approximation to the parameter posterior. See 'df' argument in |
grid.eps |
Grid size for |
crit |
Criterion for deciding "significance" of quadrature points
(difference in log-likelihood). See |
scaleSim |
Scale multiplier for the covariance matrix of the t approximation. See 'scale' argument in |
quad.ask |
Logical, for method='quadrature'. Whether or not the sampler should ask if quadrature sampling should take place. It is used to stop the sampling if the number of likelihood evaluations would be extreme. Default: FALSE. Ignored if |
force.quad |
A logical indicating whether or not to force the execution
of the quadrature method for large parameter vectors. See |
fullPost |
Logical indicating whether to draw parameter values as well to simulate full posterior. See |
dfPostIS |
Degrees of freedom for multivariate t distribution approximation to parameter posterior. See 'df' argument in |
scalePostIS |
Extra scaling factor for t distribution approximation. See 'scale' argument in |
thetaSamp |
If multiple parameter samples are available in |
export |
Character vector of the names of any additional objects or functions in the global environment that are used in |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy. See |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
miData
can either be a crwData
or crwHierData
object (as returned by crawlWrap
), a crwSim
or crwHierSim
object (as returned by MIfitHMM
when fit=FALSE
),
or a list of momentuHMMData
or momentuHierHMMData
objects (e.g., each element of the list as returned by prepData
).
If miData
is a crwData
(or crwHierData
) object, MIfitHMM
uses a combination of
crwSimulator
, crwPostIS
, prepData
, and fitHMM
to draw nSims
realizations of the position process
and fit the specified HMM to each imputation of the data. The vast majority of MIfitHMM
arguments are identical to the corresponding arguments from these functions.
If miData
is a crwData
or crwHierData
object, nSims
determines both the number of realizations of the position process to draw
(using crwSimulator
and crwPostIS
) as well as the number of HMM fits.
If miData
is a crwSim
(or crwHierSim
) object or a list of momentuHMMData
(or momentuHierHMMData
) object(s), the specified HMM will simply be fitted to each of the momentuHMMData
(or momentuHierHMMData
) objects
and all arguments related to crwSimulator
, crwPostIS
, or prepData
are ignored.
If nSims>1
, poolEstimates=TRUE
, and fit=TRUE
, a miHMM
object, i.e., a list consisting of:
miSum |
|
HMMfits |
List of length |
If poolEstimates=FALSE
and fit=TRUE
, a list of length nSims
consisting of momentuHMM
objects is returned.
However, if fit=FALSE
and miData
is a crwData
object, then MIfitHMM
returns a crwSim
object, i.e., a list containing miData
(a list of
momentuHMMData
objects) and crwSimulator
(a list of crwSimulator
objects),and most other arguments related to fitHMM
are ignored.
Hooten M.B., Johnson D.S., McClintock B.T., Morales J.M. 2017. Animal Movement: Statistical Models for Telemetry Data. CRC Press, Boca Raton.
McClintock B.T. 2017. Incorporating telemetry error into hidden Markov movement models using multiple imputation. Journal of Agricultural, Biological, and Environmental Statistics.
crawlWrap
, crwPostIS
, crwSimulator
, fitHMM
, getParDM
, MIpool
, prepData
# Don't run because it takes too long on a single core ## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # create crwData object by fitting crwMLE models to obsData and predict locations # at default intervals for both individuals crwOut <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # HMM specifications nbStates <- 2 stepDist <- "gamma" angleDist <- "vm" mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar0 <- c(mu0,sigma0) anglePar0 <- c(-pi/2,pi/2,kappa0) formula <- ~cov1+cos(cov2) nbCovs <- 2 beta0 <- matrix(c(rep(-1.5,nbStates*(nbStates-1)),rep(0,nbStates*(nbStates-1)*nbCovs)), nrow=nbCovs+1,byrow=TRUE) # first fit HMM to best predicted position process bestData<-prepData(crwOut,covNames=c("cov1","cov2")) bestFit<-fitHMM(bestData, nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar0,angle=anglePar0),beta0=beta0, formula=formula,estAngleMean=list(angle=TRUE)) print(bestFit) # extract estimates from 'bestFit' bPar0 <- getPar(bestFit) # Fit nSims=5 imputations of the position process miFits<-MIfitHMM(miData=crwOut,nSims=5, nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=bPar0$Par,beta0=bPar0$beta,delta0=bPar0$delta, formula=formula,estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # print pooled estimates print(miFits) ## End(Not run)
# Don't run because it takes too long on a single core ## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # create crwData object by fitting crwMLE models to obsData and predict locations # at default intervals for both individuals crwOut <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # HMM specifications nbStates <- 2 stepDist <- "gamma" angleDist <- "vm" mu0 <- c(20,70) sigma0 <- c(10,30) kappa0 <- c(1,1) stepPar0 <- c(mu0,sigma0) anglePar0 <- c(-pi/2,pi/2,kappa0) formula <- ~cov1+cos(cov2) nbCovs <- 2 beta0 <- matrix(c(rep(-1.5,nbStates*(nbStates-1)),rep(0,nbStates*(nbStates-1)*nbCovs)), nrow=nbCovs+1,byrow=TRUE) # first fit HMM to best predicted position process bestData<-prepData(crwOut,covNames=c("cov1","cov2")) bestFit<-fitHMM(bestData, nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=list(step=stepPar0,angle=anglePar0),beta0=beta0, formula=formula,estAngleMean=list(angle=TRUE)) print(bestFit) # extract estimates from 'bestFit' bPar0 <- getPar(bestFit) # Fit nSims=5 imputations of the position process miFits<-MIfitHMM(miData=crwOut,nSims=5, nbStates=nbStates,dist=list(step=stepDist,angle=angleDist), Par0=bPar0$Par,beta0=bPar0$beta,delta0=bPar0$delta, formula=formula,estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # print pooled estimates print(miFits) ## End(Not run)
miHMM
objectsConstructor of miHMM
objects
miHMM(m)
miHMM(m)
m |
A list with attributes
|
An object miHMM
.
Calculate pooled parameter estimates and states across multiple imputations
MIpool(im, alpha = 0.95, ncores = 1, covs = NULL, na.rm = FALSE)
MIpool(im, alpha = 0.95, ncores = 1, covs = NULL, na.rm = FALSE)
im |
List comprised of |
alpha |
Significance level for calculating confidence intervals of pooled estimates (including location error ellipses). Default: 0.95. |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
covs |
Data frame consisting of a single row indicating the covariate values to be used in the calculation of pooled natural parameters.
For any covariates that are not specified using |
na.rm |
Logical indicating whether or not to exclude model fits with |
Pooled estimates, standard errors, and confidence intervals are calculated using standard multiple imputation formulas. Working scale parameters are pooled
using MIcombine
and t-distributed confidence intervals. Natural scale parameters and normally-distributed confidence intervals are calculated by transforming the pooled working scale parameters
and, if applicable, are based on covariate means across all imputations (and/or values specified in covs
).
The calculation of pooled error ellipses uses dataEllipse
from the car
package. The suggested package car
is not automatically imported by momentuHMM
and must be installed in order to calculate error ellipses. A warning will be triggered if the car
package is required but not installed.
Note that pooled estimates for timeInStates
and stateProbs
do not include within-model uncertainty and are based entirely on across-model variability.
A miSum
object, i.e., a list comprised of model and pooled parameter summaries, including data
(averaged across imputations), conditions
, Par
, and MIcombine
(as returned by MIcombine
for working parameters).
miSum$Par
is a list comprised of:
beta |
Pooled estimates for the working parameters |
real |
Estimates for the natural parameters based on pooled working parameters and covariate means (or |
timeInStates |
The proportion of time steps assigned to each state |
states |
The most freqent state assignment for each time step based on the |
stateProbs |
Pooled state probability estimates for each time step |
mixtureProbs |
Pooled mixture probabilities for each individual (only applies if |
hierStateProbs |
Pooled state probability estimates for each time step at each level of the hierarchy (only applies if |
## Not run: # Extract data and crawl inputs from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # Pool estimates miSum <- MIpool(HMMfits) print(miSum) ## End(Not run)
## Not run: # Extract data and crawl inputs from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # Pool estimates miSum <- MIpool(HMMfits) print(miSum) ## End(Not run)
miSum
objectsConstructor of miSum
objects
miSum(m)
miSum(m)
m |
A list of attributes required for multiple imputation summaries: |
An object miSum
.
For a fitted model, this function computes the probability of each individual being in a particular mixture
mixtureProbs(m, getCI = FALSE, alpha = 0.95)
mixtureProbs(m, getCI = FALSE, alpha = 0.95)
m |
|
getCI |
Logical indicating whether to calculate standard errors and logit-transformed confidence intervals for fitted |
alpha |
Significance level of the confidence intervals (if |
When getCI=TRUE
, it can take a while for large data sets and/or a large number of mixtures because the model likelihood for each individual must be repeatedly evaluated in order to numerically approximate the SEs.
The matrix of individual mixture probabilities, with element [i,j] the probability of individual i being in mixture j
Maruotti, A., and T. Ryden. 2009. A semiparametric approach to hidden Markov models under longitudinal observations. Statistics and Computing 19: 381-393.
## Not run: nObs <- 100 nbAnimals <- 20 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.1,2)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs*nbAnimals/2))) formulaPi <- ~ sex + 0 # Females more likely in mixture 1, males more likely in mixture 2 beta <- list(beta=matrix(c(-1.5,-0.5,-1.5,-3),2,2), pi=matrix(c(-2,2),2,1,dimnames=list(c("sexF","sexM"),"mix2"))) data.mix<-simData(nbAnimals=nbAnimals,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, beta=beta,formulaPi=formulaPi,mixtures=2,covs=cov) Par0 <- list(step=Par$step, angle=Par$angle[3:4]) m.mix <- fitHMM(data.mix, nbStates=2, dist=dist, Par0 = Par0, beta0=beta,formulaPi=formulaPi,mixtures=2) mixProbs <- mixtureProbs(m.mix, getCI=TRUE) ## End(Not run)
## Not run: nObs <- 100 nbAnimals <- 20 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.1,2)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs*nbAnimals/2))) formulaPi <- ~ sex + 0 # Females more likely in mixture 1, males more likely in mixture 2 beta <- list(beta=matrix(c(-1.5,-0.5,-1.5,-3),2,2), pi=matrix(c(-2,2),2,1,dimnames=list(c("sexF","sexM"),"mix2"))) data.mix<-simData(nbAnimals=nbAnimals,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, beta=beta,formulaPi=formulaPi,mixtures=2,covs=cov) Par0 <- list(step=Par$step, angle=Par$angle[3:4]) m.mix <- fitHMM(data.mix, nbStates=2, dist=dist, Par0 = Par0, beta0=beta,formulaPi=formulaPi,mixtures=2) mixProbs <- mixtureProbs(m.mix, getCI=TRUE) ## End(Not run)
momentuHierHMM
objectsConstructor of momentuHierHMM
objects
momentuHierHMM(m)
momentuHierHMM(m)
m |
A list of attributes of the fitted model: |
An object momentuHierHMM
.
momentuHierHMMData
objectsConstructor of momentuHierHMMData
objects
momentuHierHMMData(data)
momentuHierHMMData(data)
data |
A dataframe containing: |
An object momentuHierHMMData
.
momentuHMM
objectsConstructor of momentuHMM
objects
momentuHMM(m)
momentuHMM(m)
m |
A list of attributes of the fitted model: |
An object momentuHMM
.
momentuHMMData
objectsConstructor of momentuHMMData
objects
momentuHMMData(data)
momentuHMMData(data)
data |
A dataframe containing: |
An object momentuHMMData
.
Scales each data stream probability distribution parameter from its natural interval to the set of real numbers, to allow for unconstrained optimization. Used during the optimization of the log-likelihood. Parameters of any data streams for which a design matrix is specified are ignored.
n2w( par, bounds, beta, delta = NULL, nbStates, estAngleMean, DM, Bndind, dist, TMB = FALSE )
n2w( par, bounds, beta, delta = NULL, nbStates, estAngleMean, DM, Bndind, dist, TMB = FALSE )
par |
Named list of vectors containing the initial parameter values for each data stream. |
bounds |
Named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
beta |
List of regression coefficients for the transition probabilities. |
delta |
Initial distribution. Default: |
nbStates |
The number of states of the HMM. |
estAngleMean |
Named list indicating whether or not to estimate the angle mean for data streams with angular distributions ('vm' and 'wrpcauchy'). |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of |
Bndind |
Named list indicating whether |
dist |
A named list indicating the probability distributions of the data streams. |
TMB |
logical indicating whether or not |
A vector of unconstrained parameters.
## Not run: m<-example$m nbStates <- 2 nbCovs <- 2 parSize <- list(step=2,angle=2) par <- list(step=c(t(m$mle$step)),angle=c(t(m$mle$angle))) bounds <- m$conditions$bounds beta <- matrix(rnorm(6),ncol=2,nrow=3) delta <- c(0.6,0.4) #working parameters wpar <- momentuHMM:::n2w(par,bounds,list(beta=beta),log(delta[-1]/delta[1]),nbStates, m$conditions$estAngleMean,NULL,m$conditions$Bndind, m$conditions$dist) #natural parameter p <- momentuHMM:::w2n(wpar,bounds,parSize,nbStates,nbCovs,m$conditions$estAngleMean, m$conditions$circularAngleMean,lapply(m$conditions$dist,function(x) x=="vmConsensus"), m$conditions$stationary,m$conditions$fullDM, m$conditions$DMind,1,m$conditions$dist,m$conditions$Bndind, matrix(1,nrow=length(unique(m$data$ID)),ncol=1),covsDelta=m$covsDelta, workBounds=m$conditions$workBounds) ## End(Not run)
## Not run: m<-example$m nbStates <- 2 nbCovs <- 2 parSize <- list(step=2,angle=2) par <- list(step=c(t(m$mle$step)),angle=c(t(m$mle$angle))) bounds <- m$conditions$bounds beta <- matrix(rnorm(6),ncol=2,nrow=3) delta <- c(0.6,0.4) #working parameters wpar <- momentuHMM:::n2w(par,bounds,list(beta=beta),log(delta[-1]/delta[1]),nbStates, m$conditions$estAngleMean,NULL,m$conditions$Bndind, m$conditions$dist) #natural parameter p <- momentuHMM:::w2n(wpar,bounds,parSize,nbStates,nbCovs,m$conditions$estAngleMean, m$conditions$circularAngleMean,lapply(m$conditions$dist,function(x) x=="vmConsensus"), m$conditions$stationary,m$conditions$fullDM, m$conditions$DMind,1,m$conditions$dist,m$conditions$Bndind, matrix(1,nrow=length(unique(m$data$ID)),ncol=1),covsDelta=m$covsDelta, workBounds=m$conditions$workBounds) ## End(Not run)
Negative log-likelihood function
nLogLike( optPar, nbStates, formula, bounds, parSize, data, dist, covs, estAngleMean, circularAngleMean, consensus, zeroInflation, oneInflation, stationary = FALSE, fullDM, DMind, Bndind, knownStates, fixPar, wparIndex, nc, meanind, covsDelta, workBounds, prior = NULL, betaCons = NULL, betaRef, deltaCons = NULL, optInd = NULL, recovs = NULL, g0covs = NULL, mixtures = 1, covsPi, recharge = NULL, aInd = aInd, CT = FALSE, dtIndex = NULL, kappa = Inf, crwST = FALSE )
nLogLike( optPar, nbStates, formula, bounds, parSize, data, dist, covs, estAngleMean, circularAngleMean, consensus, zeroInflation, oneInflation, stationary = FALSE, fullDM, DMind, Bndind, knownStates, fixPar, wparIndex, nc, meanind, covsDelta, workBounds, prior = NULL, betaCons = NULL, betaRef, deltaCons = NULL, optInd = NULL, recovs = NULL, g0covs = NULL, mixtures = 1, covsPi, recharge = NULL, aInd = aInd, CT = FALSE, dtIndex = NULL, kappa = Inf, crwST = FALSE )
optPar |
Vector of working parameters. |
nbStates |
Number of states of the HMM. |
formula |
Regression formula for the transition probability covariates. |
bounds |
Named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
parSize |
Named list indicating the number of natural parameters of the data stream probability distributions |
data |
An object |
dist |
Named list indicating the probability distributions of the data streams. |
covs |
data frame containing the beta model covariates (if any) |
estAngleMean |
Named list indicating whether or not to estimate the angle mean for data streams with angular distributions ('vm' and 'wrpcauchy'). |
circularAngleMean |
Named list indicating whether to use circular-linear or circular-circular
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. See |
consensus |
Named list indicating whether to use the circular-circular regression consensus model |
zeroInflation |
Named list of logicals indicating whether the probability distributions of the data streams are zero-inflated. |
oneInflation |
Named list of logicals indicating whether the probability distributions of the data streams are one-inflated. |
stationary |
|
fullDM |
Named list containing the full (i.e. not shorthand) design matrix for each data stream. |
DMind |
Named list indicating whether |
Bndind |
Named list indicating whether |
knownStates |
Vector of values of the state process which are known prior to fitting the model (if any). |
fixPar |
Vector of working parameters which are assumed known prior to fitting the model (NA indicates parameters is to be estimated). |
wparIndex |
Vector of indices for the elements of |
nc |
indicator for zeros in fullDM |
meanind |
index for circular-circular regression mean angles with at least one non-zero entry in fullDM |
covsDelta |
data frame containing the delta model covariates (if any) |
workBounds |
named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters |
prior |
A function that returns the log-density of the working scale parameter prior distribution(s) |
betaCons |
Matrix of the same dimension as |
betaRef |
Indices of reference elements for t.p.m. multinomial logit link. |
deltaCons |
Matrix of the same dimension as |
optInd |
indices of constrained parameters |
recovs |
data frame containing the recharge model theta covariates (if any) |
g0covs |
data frame containing the recharge model g0 covariates (if any) |
mixtures |
Number of mixtures for the state transition probabilities |
covsPi |
data frame containing the pi model covariates |
recharge |
recharge model specification (only used for hierarchical models) |
aInd |
vector of indices of first observation for each animal |
CT |
logical indicating whether to fit discrete-time approximation of a continuous-time model |
dtIndex |
time difference index for calculating transition probabilities of hierarchical continuous-time models |
kappa |
maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is |
crwST |
logical indicating whether to fit a correlated step and turn random walk model |
The negative log-likelihood of the parameters given the data.
## Not run: # data is a momentuHMMData object (as returned by prepData), automatically loaded with the package data <- example$m$data m<-example$m Par <- getPar(m) nbStates <- length(m$stateNames) inputs <- momentuHMM:::checkInputs(nbStates,m$conditions$dist,Par$Par,m$conditions$estAngleMean, m$conditions$circularAngleMean,m$conditions$zeroInflation,m$conditions$oneInflation, m$conditions$DM,m$conditions$userBounds, m$stateNames) wpar <- momentuHMM:::n2w(Par$Par,m$conditions$bounds,list(beta=Par$beta), log(Par$delta[-1]/Par$delta[1]),nbStates,m$conditions$estAngleMean, m$conditions$DM,m$conditions$Bndind, m$conditions$dist) l <- momentuHMM:::nLogLike(wpar,nbStates,m$conditions$formula,m$conditions$bounds, inputs$p$parSize,data,inputs$dist,model.matrix(m$conditions$formula,data), m$conditions$estAngleMean,m$conditions$circularAngleMean,inputs$consensus, m$conditions$zeroInflation,m$conditions$oneInflation,m$conditions$stationary, m$conditions$fullDM,m$conditions$DMind, m$conditions$Bndind,m$knownStates,unlist(m$conditions$fixPar), m$conditions$wparIndex,covsDelta=m$covsDelta,workBounds=m$conditions$workBounds, betaRef=m$conditions$betaRef,covsPi=m$covsPi) ## End(Not run)
## Not run: # data is a momentuHMMData object (as returned by prepData), automatically loaded with the package data <- example$m$data m<-example$m Par <- getPar(m) nbStates <- length(m$stateNames) inputs <- momentuHMM:::checkInputs(nbStates,m$conditions$dist,Par$Par,m$conditions$estAngleMean, m$conditions$circularAngleMean,m$conditions$zeroInflation,m$conditions$oneInflation, m$conditions$DM,m$conditions$userBounds, m$stateNames) wpar <- momentuHMM:::n2w(Par$Par,m$conditions$bounds,list(beta=Par$beta), log(Par$delta[-1]/Par$delta[1]),nbStates,m$conditions$estAngleMean, m$conditions$DM,m$conditions$Bndind, m$conditions$dist) l <- momentuHMM:::nLogLike(wpar,nbStates,m$conditions$formula,m$conditions$bounds, inputs$p$parSize,data,inputs$dist,model.matrix(m$conditions$formula,data), m$conditions$estAngleMean,m$conditions$circularAngleMean,inputs$consensus, m$conditions$zeroInflation,m$conditions$oneInflation,m$conditions$stationary, m$conditions$fullDM,m$conditions$DMind, m$conditions$Bndind,m$knownStates,unlist(m$conditions$fixPar), m$conditions$wparIndex,covsDelta=m$covsDelta,workBounds=m$conditions$workBounds, betaRef=m$conditions$betaRef,covsPi=m$covsPi) ## End(Not run)
Computation of the negative log-likelihood (forward algorithm - written in C++)
nLogLike_rcpp( nbStates, covs, data, dataNames, dist, Par, aInd, zeroInflation, oneInflation, stationary, knownStates, betaRef, mixtures, dtIndex, CT = FALSE, kappa = NA_real_ )
nLogLike_rcpp( nbStates, covs, data, dataNames, dist, Par, aInd, zeroInflation, oneInflation, stationary, knownStates, betaRef, mixtures, dtIndex, CT = FALSE, kappa = NA_real_ )
nbStates |
Number of states, |
covs |
Covariates, |
data |
A |
dataNames |
Character vector containing the names of the data streams, |
dist |
Named list indicating the probability distributions of the data streams. |
Par |
Named list containing the state-dependent parameters of the data streams, matrix of regression coefficients for the transition probabilities ('beta'), and initial distribution ('delta'). |
aInd |
Vector of indices of the rows at which the data switches to another animal |
zeroInflation |
Named list of logicals indicating whether the probability distributions of the data streams are zero-inflated. |
oneInflation |
Named list of logicals indicating whether the probability distributions of the data streams are one-inflated. |
stationary |
|
knownStates |
Vector of values of the state process which are known prior to fitting the model (if any). Default: NULL (states are not known). This should be a vector with length the number of rows of 'data'; each element should either be an integer (the value of the known states) or NA if the state is not known. |
betaRef |
Indices of reference elements for t.p.m. multinomial logit link. |
mixtures |
Number of mixtures for the state transition probabilities |
dtIndex |
time difference index for calculating transition probabilities of hierarchical continuous-time models |
CT |
logical indicating whether to fit discrete-time approximation of a continuous-time model |
kappa |
maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is |
Negative log-likelihood
Parameters definition
parDef( dist, nbStates, estAngleMean, zeroInflation, oneInflation, DM, userBounds = NULL )
parDef( dist, nbStates, estAngleMean, zeroInflation, oneInflation, DM, userBounds = NULL )
dist |
Named list indicating the probability distributions of the data streams. |
nbStates |
Number of states of the HMM. |
estAngleMean |
Named list indicating whether or not to estimate the angle mean for data streams with angular distributions ('vm' and 'wrpcauchy'). |
zeroInflation |
Named list of logicals indicating whether the probability distributions of the data streams should be zero-inflated. |
oneInflation |
Named list of logicals indicating whether the probability distributions of the data streams are one-inflated. |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability
distribution parameters for each data stream. For example, for a 2-state model using the wrapped Cauchy ('wrpcauchy') distribution for
a data stream named 'angle' with |
A list of:
parSize |
Named list indicating the number of natural parameters of the data stream probability distributions. |
bounds |
Named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
parNames |
Names of parameters of the probability distribution for each data stream. |
Bndind |
Named list indicating whether |
## Not run: pD<-momentuHMM:::parDef(list(step="gamma",angle="wrpcauchy"), nbStates=2,list(step=FALSE,angle=FALSE),list(step=FALSE,angle=FALSE), list(step=FALSE,angle=FALSE),NULL,NULL) ## End(Not run)
## Not run: pD<-momentuHMM:::parDef(list(step="gamma",angle="wrpcauchy"), nbStates=2,list(step=FALSE,angle=FALSE),list(step=FALSE,angle=FALSE), list(step=FALSE,angle=FALSE),NULL,NULL) ## End(Not run)
crwData
Plot observed locations, error ellipses (if applicable), predicted locations, and prediction intervals from crwData
or crwHierData
object.
## S3 method for class 'crwData' plot( x, animals = NULL, compact = FALSE, ask = TRUE, plotEllipse = TRUE, crawlPlot = FALSE, ... )
## S3 method for class 'crwData' plot( x, animals = NULL, compact = FALSE, ask = TRUE, plotEllipse = TRUE, crawlPlot = FALSE, ... )
x |
An object |
animals |
Vector of indices or IDs of animals for which information will be plotted.
Default: |
compact |
|
ask |
If |
plotEllipse |
If |
crawlPlot |
Logical indicating whether or not to create individual plots using |
... |
Further arguments for passing to |
In order for error ellipses to be plotted, the names for the semi-major axis, semi-minor axis, and
orientation in x$crwPredict
must respectively be error_semimajor_axis
, error_semiminor_axis
,
and error_ellipse_orientation
.
If the crwData
(or crwHierData
) object was created using data generated by
simData
(or simHierData
) or simObsData
, then the true locations (mux
,muy
) are also plotted.
## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # create crwData object crwOut <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) plot(crwOut,compact=TRUE,ask=FALSE,plotEllipse=FALSE) ## End(Not run)
## Not run: # extract simulated obsData from example data obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # create crwData object crwOut <- crawlWrap(obsData=obsData, theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) plot(crwOut,compact=TRUE,ask=FALSE,plotEllipse=FALSE) ## End(Not run)
miHMM
For multiple imputation analyses, plot the pooled data stream densities over histograms of the data, probability distribution parameters and transition probabilities as functions of the covariates, and maps of the animals' tracks colored by the decoded states.
## S3 method for class 'miHMM' plot( x, animals = NULL, covs = NULL, ask = TRUE, breaks = "Sturges", hist.ylim = NULL, sepAnimals = FALSE, sepStates = FALSE, col = NULL, cumul = TRUE, plotTracks = TRUE, plotCI = FALSE, alpha = 0.95, plotStationary = FALSE, plotEllipse = TRUE, return = FALSE, ... )
## S3 method for class 'miHMM' plot( x, animals = NULL, covs = NULL, ask = TRUE, breaks = "Sturges", hist.ylim = NULL, sepAnimals = FALSE, sepStates = FALSE, col = NULL, cumul = TRUE, plotTracks = TRUE, plotCI = FALSE, alpha = 0.95, plotStationary = FALSE, plotEllipse = TRUE, return = FALSE, ... )
x |
Object |
animals |
Vector of indices or IDs of animals for which information will be plotted.
Default: |
covs |
Data frame consisting of a single row indicating the covariate values to be used in plots. If none are specified, the means of any covariates appearing in the model are used (unless covariate is a factor, in which case the first factor appearing in the data is used). |
ask |
If |
breaks |
Histogram parameter. See |
hist.ylim |
Parameter |
sepAnimals |
If |
sepStates |
If |
col |
Vector or colors for the states (one color per state). |
cumul |
If TRUE, the sum of weighted densities is plotted (default). |
plotTracks |
If TRUE, the Viterbi-decoded tracks are plotted (default). |
plotCI |
Logical indicating whether to include confidence intervals in natural parameter plots (default: FALSE) |
alpha |
Significance level of the confidence intervals (if |
plotStationary |
Logical indicating whether to plot the stationary state probabilities as a function of any covariates (default: FALSE) |
plotEllipse |
Logical indicating whether to plot error ellipses around imputed location means. Default: TRUE. |
return |
Logical indicating whether to return a list containing estimates, SEs, CIs, and covariate values used to create the plots for each mixture and state. Ignored if |
... |
Additional arguments passed to |
The state-dependent densities are weighted by the frequency of each state in the most
probable state sequence (decoded with the function viterbi
for each imputation). For example, if the
most probable state sequence indicates that one third of observations correspond to the first
state, and two thirds to the second state, the plots of the densities in the first state are
weighted by a factor 1/3, and in the second state by a factor 2/3.
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) miHMM <- momentuHMM:::miHMM(list(miSum=MIpool(HMMfits),HMMfits=HMMfits)) plot(miHMM) ## End(Not run)
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) miHMM <- momentuHMM:::miHMM(list(miSum=MIpool(HMMfits),HMMfits=HMMfits)) plot(miHMM) ## End(Not run)
miSum
Plot the fitted step and angle densities over histograms of the data, transition probabilities as functions of the covariates, and maps of the animals' tracks colored by the decoded states.
## S3 method for class 'miSum' plot( x, animals = NULL, covs = NULL, ask = TRUE, breaks = "Sturges", hist.ylim = NULL, sepAnimals = FALSE, sepStates = FALSE, col = NULL, cumul = TRUE, plotTracks = TRUE, plotCI = FALSE, alpha = 0.95, plotStationary = FALSE, plotEllipse = TRUE, return = FALSE, ... )
## S3 method for class 'miSum' plot( x, animals = NULL, covs = NULL, ask = TRUE, breaks = "Sturges", hist.ylim = NULL, sepAnimals = FALSE, sepStates = FALSE, col = NULL, cumul = TRUE, plotTracks = TRUE, plotCI = FALSE, alpha = 0.95, plotStationary = FALSE, plotEllipse = TRUE, return = FALSE, ... )
x |
Object |
animals |
Vector of indices or IDs of animals for which information will be plotted.
Default: |
covs |
Data frame consisting of a single row indicating the covariate values to be used in plots. If none are specified, the means of any covariates appearing in the model are used (unless covariate is a factor, in which case the first factor appearing in the data is used). |
ask |
If |
breaks |
Histogram parameter. See |
hist.ylim |
Parameter |
sepAnimals |
If |
sepStates |
If |
col |
Vector or colors for the states (one color per state). |
cumul |
If TRUE, the sum of weighted densities is plotted (default). |
plotTracks |
If TRUE, the Viterbi-decoded tracks are plotted (default). |
plotCI |
Logical indicating whether to include confidence intervals in natural parameter plots (default: FALSE) |
alpha |
Significance level of the confidence intervals (if |
plotStationary |
Logical indicating whether to plot the stationary state probabilities as a function of any covariates (default: FALSE) |
plotEllipse |
Logical indicating whether to plot error ellipses around imputed location means. Default: TRUE. |
return |
Logical indicating whether to return a list containing estimates, SEs, CIs, and covariate values used to create the plots for each mixture and state. Ignored if |
... |
Additional arguments passed to |
The state-dependent densities are weighted by the frequency of each state in the most
probable state sequence (decoded with the function viterbi
for each imputation). For example, if the
most probable state sequence indicates that one third of observations correspond to the first
state, and two thirds to the second state, the plots of the densities in the first state are
weighted by a factor 1/3, and in the second state by a factor 2/3.
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # Pool estimates miSum <- MIpool(HMMfits) plot(miSum) ## End(Not run)
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # Pool estimates miSum <- MIpool(HMMfits) plot(miSum) ## End(Not run)
momentuHMM
Plot the fitted step and angle densities over histograms of the data, transition probabilities as functions of the covariates, and maps of the animals' tracks colored by the decoded states.
## S3 method for class 'momentuHMM' plot( x, animals = NULL, covs = NULL, ask = TRUE, breaks = "Sturges", hist.ylim = NULL, sepAnimals = FALSE, sepStates = FALSE, col = NULL, cumul = TRUE, plotTracks = TRUE, plotCI = FALSE, alpha = 0.95, plotStationary = FALSE, return = FALSE, ... )
## S3 method for class 'momentuHMM' plot( x, animals = NULL, covs = NULL, ask = TRUE, breaks = "Sturges", hist.ylim = NULL, sepAnimals = FALSE, sepStates = FALSE, col = NULL, cumul = TRUE, plotTracks = TRUE, plotCI = FALSE, alpha = 0.95, plotStationary = FALSE, return = FALSE, ... )
x |
Object |
animals |
Vector of indices or IDs of animals for which information will be plotted.
Default: |
covs |
Data frame consisting of a single row indicating the covariate values to be used in plots. If none are specified, the means of any covariates appearing in the model are used (unless covariate is a factor, in which case the first factor in the data is used). |
ask |
If |
breaks |
Histogram parameter. See |
hist.ylim |
An optional named list of vectors specifying |
sepAnimals |
If |
sepStates |
If |
col |
Vector or colors for the states (one color per state). |
cumul |
If TRUE, the sum of weighted densities is plotted (default). |
plotTracks |
If TRUE, the Viterbi-decoded tracks are plotted (default). |
plotCI |
Logical indicating whether to include confidence intervals in natural parameter plots (default: FALSE) |
alpha |
Significance level of the confidence intervals (if |
plotStationary |
Logical indicating whether to plot the stationary state probabilities as a function of any covariates (default: FALSE). Ignored unless covariate are included in |
return |
Logical indicating whether to return a list containing: 1) estimates, SEs, CIs, and covariate values used to create the plots for each mixture and state ( |
... |
Additional arguments passed to |
The state-dependent densities are weighted by the frequency of each state in the most
probable state sequence (decoded with the function viterbi
). For example, if the
most probable state sequence indicates that one third of observations correspond to the first
state, and two thirds to the second state, the plots of the densities in the first state are
weighted by a factor 1/3, and in the second state by a factor 2/3.
Confidence intervals for natural parameters are calculated from the working parameter point and covariance estimates
using finite-difference approximations of the first derivative for the transformation (see grad
).
For example, if dN
is the numerical approximation of the first derivative of the transformation N = exp(x_1 * B_1 + x_2 * B_2)
for covariates (x_1, x_2) and working parameters (B_1, B_2), then
var(N)=dN %*% Sigma %*% dN
, where Sigma=cov(B_1,B_2)
, and normal confidence intervals can be
constructed as N +/- qnorm(1-(1-alpha)/2) * se(N)
.
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m plot(m,ask=TRUE,animals=1,breaks=20,plotCI=TRUE)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m plot(m,ask=TRUE,animals=1,breaks=20,plotCI=TRUE)
momentuHMMData
or momentuHierHMMData
Plot momentuHMMData
or momentuHierHMMData
## S3 method for class 'momentuHMMData' plot( x, dataNames = c("step", "angle"), animals = NULL, compact = FALSE, ask = TRUE, breaks = "Sturges", ... )
## S3 method for class 'momentuHMMData' plot( x, dataNames = c("step", "angle"), animals = NULL, compact = FALSE, ask = TRUE, breaks = "Sturges", ... )
x |
An object |
dataNames |
Names of the variables to plot. Default is |
animals |
Vector of indices or IDs of animals for which information will be plotted.
Default: |
compact |
|
ask |
If |
breaks |
Histogram parameter. See |
... |
Currently unused. For compatibility with generic method. |
# data is a momentuHMMData object (as returned by prepData), automatically loaded with the package data <- example$m$data plot(data,dataNames=c("step","angle","cov1","cov2"), compact=TRUE,breaks=20,ask=FALSE)
# data is a momentuHMMData object (as returned by prepData), automatically loaded with the package data <- example$m$data plot(data,dataNames=c("step","angle","cov1","cov2"), compact=TRUE,breaks=20,ask=FALSE)
Plots time series, qq-plots (against the standard normal distribution) using qqPlot
, and sample
ACF functions of the pseudo-residuals for each data stream
plotPR(m, lag.max = NULL, ncores = 1)
plotPR(m, lag.max = NULL, ncores = 1)
m |
A |
lag.max |
maximum lag at which to calculate the acf. See |
ncores |
number of cores to use for parallel processing |
If some turning angles in the data are equal to pi, the corresponding pseudo-residuals will not be included. Indeed, given that the turning angles are defined on (-pi,pi], an angle of pi results in a pseudo-residual of +Inf (check Section 6.2 of reference for more information on the computation of pseudo-residuals).
If some data streams are zero-inflated and/or one-inflated, the corresponding pseudo- residuals are shown as segments, because pseudo-residuals for discrete data are defined as segments (see Zucchini and MacDonald, 2009, Section 6.2).
For multiple imputation analyses, if m
is a miHMM
object or a list of momentuHMM
objects, then
the pseudo-residuals are individually calculated and plotted for each model fit. Note that pseudo-residuals for miSum
objects (as returned by MIpool
) are based on pooled parameter
estimates and the means of the data values across all imputations (and therefore may not be particularly meaningful).
Zucchini, W. and MacDonald, I.L. 2009. Hidden Markov Models for Time Series: An Introduction Using R. Chapman & Hall (London).
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m plotPR(m)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m plotPR(m)
Plot tracking data on a satellite map. This function plots coordinates in longitude
and latitude (not UTM), so if data
coordinates are not provided in longitude and latitude, then the coordinate reference system must be provided using the projargs
argument. This function uses the package ggmap
to fetch a satellite image from Google. An Internet connection is required to use
this function.
plotSat( data, zoom = NULL, location = NULL, segments = TRUE, compact = TRUE, col = NULL, alpha = 1, size = 1, shape = 16, states = NULL, animals = NULL, ask = TRUE, return = FALSE, stateNames = NULL, projargs = NULL )
plotSat( data, zoom = NULL, location = NULL, segments = TRUE, compact = TRUE, col = NULL, alpha = 1, size = 1, shape = 16, states = NULL, animals = NULL, ask = TRUE, return = FALSE, stateNames = NULL, projargs = NULL )
data |
Data frame or |
zoom |
The zoom level, as defined for |
location |
Location of the center of the map to be plotted (this must be in the same coordinate reference system as |
segments |
|
compact |
|
col |
Palette of colours to use for the dots and segments. If not specified, uses default palette. |
alpha |
Transparency argument for |
size |
Size argument for |
shape |
Shape argument for |
states |
A sequence of integers, corresponding to the decoded states for these data (such that the observations are colored by states). |
animals |
Vector of indices or IDs of animals/tracks to be plotted.
Default: |
ask |
If |
return |
If |
stateNames |
Optional character vector of length |
projargs |
A character string of PROJ.4 projection arguments indicating the coordinate reference system for |
If the plot displays the message "Sorry, we have no imagery here", try a lower level of zoom.
D. Kahle and H. Wickham. ggmap: Spatial Visualization with ggplot2. The R Journal, 5(1), 144-161. URL: http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf
Plot tracking data over a raster layer.
plotSpatialCov( data, spatialCov, segments = TRUE, compact = TRUE, col = NULL, alpha = 1, size = 1, shape = 16, states = NULL, animals = NULL, ask = TRUE, return = FALSE, stateNames = NULL, colors = c("darkblue", "lightblue") )
plotSpatialCov( data, spatialCov, segments = TRUE, compact = TRUE, col = NULL, alpha = 1, size = 1, shape = 16, states = NULL, animals = NULL, ask = TRUE, return = FALSE, stateNames = NULL, colors = c("darkblue", "lightblue") )
data |
Data frame or |
spatialCov |
|
segments |
|
compact |
|
col |
Palette of colours to use for the dots and segments. If not specified, uses default palette. |
alpha |
Transparency argument for |
size |
Size argument for |
shape |
Shape argument for |
states |
A sequence of integers, corresponding to the decoded states for these data. If specified, the observations are colored by states. |
animals |
Vector of indices or IDs of animals/tracks to be plotted.
Default: |
ask |
If |
return |
If |
stateNames |
Optional character vector of length |
colors |
Colours argument for |
## Not run: stepDist <- "gamma" angleDist <- "vm" # plot simulated data over forest raster automatically loaded with the packge spatialCov<-list(forest=forest) data <- simData(nbAnimals=2,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par=list(step=c(100,1000,50,100),angle=c(0,0,0.1,5)), beta=matrix(c(5,-10,-25,50),nrow=2,ncol=2,byrow=TRUE), formula=~forest,spatialCovs=spatialCov, obsPerAnimal=225,states=TRUE) plotSpatialCov(data,forest,states=data$states) ## End(Not run)
## Not run: stepDist <- "gamma" angleDist <- "vm" # plot simulated data over forest raster automatically loaded with the packge spatialCov<-list(forest=forest) data <- simData(nbAnimals=2,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par=list(step=c(100,1000,50,100),angle=c(0,0,0.1,5)), beta=matrix(c(5,-10,-25,50),nrow=2,ncol=2,byrow=TRUE), formula=~forest,spatialCovs=spatialCov, obsPerAnimal=225,states=TRUE) plotSpatialCov(data,forest,states=data$states) ## End(Not run)
Plot the states and states probabilities.
plotStates(m, animals = NULL, ask = TRUE)
plotStates(m, animals = NULL, ask = TRUE)
m |
A |
animals |
Vector of indices or IDs of animals for which states will be plotted. |
ask |
If |
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m # plot states for first and second animals plotStates(m,animals=c(1,2))
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m # plot states for first and second animals plotStates(m,animals=c(1,2))
Plot stationary state probabilities
plotStationary( model, covs = NULL, col = NULL, plotCI = FALSE, alpha = 0.95, return = FALSE, ... )
plotStationary( model, covs = NULL, col = NULL, plotCI = FALSE, alpha = 0.95, return = FALSE, ... )
model |
|
covs |
Optional data frame consisting of a single row indicating the covariate values to be used in plots. If none are specified, the means of any covariates appearing in the model are used (unless covariate is a factor, in which case the first factor in the data is used). |
col |
Vector or colors for the states (one color per state). |
plotCI |
Logical indicating whether to include confidence intervals in plots (default: FALSE) |
alpha |
Significance level of the confidence intervals (if |
return |
Logical indicating whether to return a list containing estimates, SEs, CIs, and covariate values used to create the plots for each mixture and state. Ignored if |
... |
Additional arguments passed to |
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m plotStationary(m)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m plotStationary(m)
This wrapper function for path2ctmc
and ctmc2glm
converts a data.frame
of coordinates, other data streams, and non-spatial covariates to a momentuHMMData
object that can be passed directly to fitCTHMM
(or as a list to MIfitCTHMM
).
prepCTDS(data, ...) ## Default S3 method: prepCTDS( data, Time.unit = "auto", rast, directions = 4, zero.idx = integer(), print.iter = FALSE, interpMethod = "ShortestPath", spatialCovs = NULL, spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, covNames = NULL, ncores = 1, ... ) ## S3 method for class 'hierarchical' prepCTDS( data, Time.unit = "auto", rast, directions = 4, zero.idx = integer(), print.iter = FALSE, interpMethod = "ShortestPath", spatialCovs = NULL, spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, covNames = NULL, ncores = 1, hierLevels, coordLevel, ... )
prepCTDS(data, ...) ## Default S3 method: prepCTDS( data, Time.unit = "auto", rast, directions = 4, zero.idx = integer(), print.iter = FALSE, interpMethod = "ShortestPath", spatialCovs = NULL, spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, covNames = NULL, ncores = 1, ... ) ## S3 method for class 'hierarchical' prepCTDS( data, Time.unit = "auto", rast, directions = 4, zero.idx = integer(), print.iter = FALSE, interpMethod = "ShortestPath", spatialCovs = NULL, spatialCovs.grad = NULL, crw = TRUE, normalize.gradients = FALSE, grad.point.decreasing = FALSE, covNames = NULL, ncores = 1, hierLevels, coordLevel, ... )
data |
Either a |
... |
further arguments passed to or from other methods |
Time.unit |
Character string indicating units for time difference between observations (e.g. 'auto', 'secs', 'mins', 'hours', 'days', 'weeks'). Ignored unless |
rast |
A raster object or raster stack object that will define the discrete-space grid cells for the CTMC movement path. |
directions |
Integer. Either 4 (indicating a "Rook's neighborhood" of 4 neighboring grid cells) or 8 (indicating a "King's neighborhood" of 8 neighboring grid cells). |
zero.idx |
Integer vector of the indices of raster cells that are not passable and should be excluded. These are cells where movement should be impossible. Default is zero.idx=integer(). |
print.iter |
Logical. If true, then the progress stepping through each observed location in |
interpMethod |
Specifies interpolation method. Either "ShortestPath", which uses the shortest graphical path on the raster graph, or "LinearInterp", which linearly interpolates between observed locations. "ShortestPath" is slower, slightly more accurate, and allows for impassible barriers specified through "zero.idx". "LinearInterp" is faster but does not allow for impassible barriers. |
spatialCovs |
List of |
spatialCovs.grad |
List of |
crw |
Logical. If TRUE (default), an autocovariate is created for each cell visited in the CTMC movement path. The autocovariate is a unit-length directional vector pointing from the center of the most recent grid cell to the center of the current grid cell. |
normalize.gradients |
Logical. Default is FALSE. If TRUE, then all gradient covariates for |
grad.point.decreasing |
Logical. Default is FALSE. If TRUE, then the gradient covariates are positive in the direction of decreasing values of the covariate. If FALSE, then the gradient covariates are positive in the direction of increasing values of the covariate (like a true gradient). |
covNames |
Character vector indicating the names of any covariates in |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
hierLevels |
Character vector indicating the levels of the hierarchy and their order, from top (coarsest scale) to bottom (finest scale), that are included in |
coordLevel |
Character string indicating the level of the hierarchy for the location data. If specified, then |
Any moves to adjacent cells (as defined by directions
; see adjacent
) between times t and t+1 are assumed to occur at time t+1 (where tau
= t+1 - t), regardless of the length of the time step (tau
) and where in the current cell or in the neighboring cell that the corresponding (continuous-space) locations reside.
This is done to facilitate the temporal alignment of the discrete-space movement process with other data streams and spatio-temporal covariates. The appropriateness of this simplifying assumption will depend on tau
, the scale of movement, and the raster resolution.
If there are move(s) across multiple cells within a time step, data stream(s) other than z
(the categorical CTDS data stream) are set to NA
for these time step(s).
These NA
data stream values must be manually set based on the time spent in each cell if they are to be included in subsequent analysis.
Time-varying spatial covariates for these time step(s) are set based on the z-value for the initial cell and, if different, these must also be set manually based on z-value(s) at the time cell(s) were crossed.
All such instances are indicated wherever the 'cellCross' field is > 0.
If data
is a crwData
object, the momentuHMMData
object created by prepCTDS
is based on the best predicted
locations (i.e., crwData$crwPredict$mu.x
and crwData$crwPredict$mu.y
). Prior to using prepCTDS
, additional data streams or covariates unrelated to location (including z-values associated with
spatialCovs
and spatialCovs.grad
raster stacks or bricks) can be merged with the crwData
object using crawlMerge
.
A momentuHMMData
object of class ctds
that can be passed to fitCTHMM
(or as a list to MIfitCTHMM
), where:
ID |
The ID(s) of the observed animal(s) |
time |
Time stamp for each observation |
x |
Easting coordinate of the current cell |
y |
Northing coordinate of the current cell |
z |
Categorical CTDS data stream indicating which cell was moved to, where |
tau |
Time difference (in |
cellCross |
Integer indexing move(s) across multiple cells within a time step (if any) |
... |
Additional data streams (if any) |
... |
Covariates (if any). Covariates are needed for all neighboring cells and are indexed by |
Hanks E. M., Hooten M. B., and Alldredge M. W. 2015. Continuous-time Discrete-space Models for Animal Movement. The Annals of Applied Statistics 9:145-165
Preprocessing of the data streams, including calculation of step length, turning angle, and covariates from location data to be suitable for
analysis using fitHMM
.
prepData(data, ...) ## Default S3 method: prepData( data, type = c("UTM", "LL"), coordNames = c("x", "y"), covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, CT = FALSE, Time.name = "time", Time.unit = "auto", ... ) ## S3 method for class 'hierarchical' prepData( data, type = c("UTM", "LL"), coordNames = c("x", "y"), covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, CT = FALSE, Time.name = "time", Time.unit = "auto", hierLevels, coordLevel, ... )
prepData(data, ...) ## Default S3 method: prepData( data, type = c("UTM", "LL"), coordNames = c("x", "y"), covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, CT = FALSE, Time.name = "time", Time.unit = "auto", ... ) ## S3 method for class 'hierarchical' prepData( data, type = c("UTM", "LL"), coordNames = c("x", "y"), covNames = NULL, spatialCovs = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, altCoordNames = NULL, gradient = FALSE, CT = FALSE, Time.name = "time", Time.unit = "auto", hierLevels, coordLevel, ... )
data |
Either a data frame of data streams or a |
... |
further arguments passed to or from other methods |
type |
|
coordNames |
Names of the columns of coordinates in the |
covNames |
Character vector indicating the names of any covariates in |
spatialCovs |
List of |
centers |
2-column matrix providing the x-coordinates (column 1) and y-coordinates (column 2) for any activity centers (e.g., potential
centers of attraction or repulsion) from which distance and angle covariates will be calculated based on the location data. If no row names are provided, then generic names are generated
for the distance and angle covariates (e.g., 'center1.dist', 'center1.angle', 'center2.dist', 'center2.angle'); otherwise the covariate names are derived from the row names
of |
centroids |
List where each element is a data frame containing the x-coordinates ('x'), y-coordinates ('y'), and times (with user-specified name, e.g., ‘time’) for centroids (i.e., dynamic activity centers where the coordinates can change over time)
from which distance and angle covariates will be calculated based on the location data. If any centroids are specified, then |
angleCovs |
Character vector indicating the names of any circular-circular regression angular covariates in |
altCoordNames |
Character string indicating an alternative name for the returned location data. If provided, then |
gradient |
Logical indicating whether or not to calculate gradients of |
CT |
Logical indicating whether or not |
Time.name |
Character string indicating name of the time column in data (for continuous-time HMMs). Default: |
Time.unit |
Character string indicating units for time difference between observations (e.g. 'auto', 'secs', 'mins', 'hours', 'days', 'weeks'). Ignored unless |
hierLevels |
Character vector indicating the levels of the hierarchy and their order, from top (coarsest scale) to bottom (finest scale), that are included in |
coordLevel |
Character string indicating the level of the hierarchy for the location data. If specified, then |
If data
is a crwData
(or crwHierData
) object, the momentuHMMData
(or momentuHierHMMData
) object created by prepData
includes step lengths and turning angles calculated from the best predicted
locations (i.e., crwData$crwPredict$mu.x
and crwData$crwPredict$mu.y
). Prior to using prepData
, additional data streams or covariates unrelated to location (including z-values associated with
spatialCovs
raster stacks or bricks) can be merged with the crwData
(or crwHierData
) object using crawlMerge
.
For continuous-time data (i.e. CT=TRUE
), if data
is not a crwData
object returned by crawlWrap
, then data
must include a 'Time.name' field indicating the time column in data
.
For hierarchical data, data
must include a 'level' field indicating the level of the hierarchy for each observation, and, for location data identified by coordNames
, the coordLevel
argument must indicate the level of the hierarchy at which the location data are obtained.
An object momentuHMMData
or momentuHierHMMData
, i.e., a dataframe of:
ID |
The ID(s) of the observed animal(s) |
... |
Data streams (e.g., 'step', 'angle', etc.) |
x |
Either easting or longitude (if |
y |
Either norting or latitude (if |
... |
Covariates (if any) |
crawlMerge
, crawlWrap
, crwData
coord1 <- c(1,2,3,4,5,6,7,8,9,10) coord2 <- c(1,1,1,2,2,2,1,1,1,2) cov1 <- rnorm(10) data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1") # include additional data stream named 'omega' omega <- rbeta(10,1,1) data <- data.frame(coord1=coord1,coord2=coord2,omega=omega,cov1=cov1) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1") # include 'forest' example raster layer as covariate data <- data.frame(coord1=coord1*1000,coord2=coord2*1000) spatialCov <- list(forest=forest) d <- prepData(data,coordNames=c("coord1","coord2"),spatialCovs=spatialCov) # include 2 activity centers data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1", centers=matrix(c(0,10,0,10),2,2,dimnames=list(c("c1","c2"),NULL))) # include centroid data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1,time=1:10) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1", centroid=list(centroid=data.frame(x=coord1+rnorm(10), y=coord2+rnorm(10), time=1:10))) # Include angle covariate that needs conversion to # turning angle relative to previous movement direction u <- rnorm(10) # horizontal component v <- rnorm(10) # vertical component cov2 <- atan2(v,u) data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1,cov2=cov2) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1", angleCovs="cov2")
coord1 <- c(1,2,3,4,5,6,7,8,9,10) coord2 <- c(1,1,1,2,2,2,1,1,1,2) cov1 <- rnorm(10) data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1") # include additional data stream named 'omega' omega <- rbeta(10,1,1) data <- data.frame(coord1=coord1,coord2=coord2,omega=omega,cov1=cov1) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1") # include 'forest' example raster layer as covariate data <- data.frame(coord1=coord1*1000,coord2=coord2*1000) spatialCov <- list(forest=forest) d <- prepData(data,coordNames=c("coord1","coord2"),spatialCovs=spatialCov) # include 2 activity centers data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1", centers=matrix(c(0,10,0,10),2,2,dimnames=list(c("c1","c2"),NULL))) # include centroid data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1,time=1:10) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1", centroid=list(centroid=data.frame(x=coord1+rnorm(10), y=coord2+rnorm(10), time=1:10))) # Include angle covariate that needs conversion to # turning angle relative to previous movement direction u <- rnorm(10) # horizontal component v <- rnorm(10) # vertical component cov2 <- atan2(v,u) data <- data.frame(coord1=coord1,coord2=coord2,cov1=cov1,cov2=cov2) d <- prepData(data,coordNames=c("coord1","coord2"),covNames="cov1", angleCovs="cov2")
miHMM
Print miHMM
## S3 method for class 'miHMM' print(x, ...)
## S3 method for class 'miHMM' print(x, ...)
x |
A |
... |
Currently unused. For compatibility with generic method. |
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) miHMM <- momentuHMM:::miHMM(list(miSum=MIpool(HMMfits),HMMfits=HMMfits)) print(miHMM) ## End(Not run)
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) miHMM <- momentuHMM:::miHMM(list(miSum=MIpool(HMMfits),HMMfits=HMMfits)) print(miHMM) ## End(Not run)
miSum
Print miSum
## S3 method for class 'miSum' print(x, ...)
## S3 method for class 'miSum' print(x, ...)
x |
A |
... |
Currently unused. For compatibility with generic method. |
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # Pool estimates miSum <- MIpool(HMMfits) print(miSum) ## End(Not run)
## Not run: # Extract data from miExample obsData <- miExample$obsData # error ellipse model err.model <- list(x= ~ ln.sd.x - 1, y = ~ ln.sd.y - 1, rho = ~ error.corr) # Fit crawl to obsData crwOut <- crawlWrap(obsData,theta=c(4,0),fixPar=c(1,1,NA,NA), err.model=err.model) # Fit four imputations bPar <- miExample$bPar HMMfits <- MIfitHMM(crwOut,nSims=4,poolEstimates=FALSE, nbStates=2,dist=list(step="gamma",angle="vm"), Par0=bPar$Par,beta0=bPar$beta, formula=~cov1+cos(cov2), estAngleMean=list(angle=TRUE), covNames=c("cov1","cov2")) # Pool estimates miSum <- MIpool(HMMfits) print(miSum) ## End(Not run)
momentuHMM
Print momentuHMM
## S3 method for class 'momentuHMM' print(x, ...) ## S3 method for class 'momentuHierHMM' print(x, ...)
## S3 method for class 'momentuHMM' print(x, ...) ## S3 method for class 'momentuHierHMM' print(x, ...)
x |
A |
... |
Currently unused. For compatibility with generic method. |
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m print(m)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m print(m)
The pseudo-residuals of momentuHMM models, as described in Zucchini and McDonad (2009).
pseudoRes(m, ncores = 1)
pseudoRes(m, ncores = 1)
m |
A |
ncores |
number of cores to use for parallel processing |
If some turning angles in the data are equal to pi, the corresponding pseudo-residuals will not be included. Indeed, given that the turning angles are defined on (-pi,pi], an angle of pi results in a pseudo-residual of +Inf (check Section 6.2 of reference for more information on the computation of pseudo-residuals).
A continuity adjustment (adapted from Harte 2017) is made for discrete probability distributions. When the data are near the boundary (e.g. 0 for “pois”; 0 and 1 for “bern”), then the pseudo residuals can be a poor indicator of lack of fit.
For multiple imputation analyses, if m
is a miHMM
object or a list of momentuHMM
objects, then
the pseudo-residuals are individually calculated for each model fit. Note that pseudo-residuals for miSum
objects (as returned by MIpool
) are based on pooled parameter
estimates and the means of the data values across all imputations (and therefore may not be particularly meaningful).
If m
is a momentuHMM
, miHMM
, or miSum
object, a list of pseudo-residuals for each data stream (e.g., 'stepRes', 'angleRes') is returned.
If m
is a list of momentuHMM
objects, then a list of length length(m)
is returned where each element is a list of pseudo-residuals for each data stream.
Harte, D. 2017. HiddenMarkov: Hidden Markov Models. R package version 1.8-8.
Zucchini, W. and MacDonald, I.L. 2009. Hidden Markov Models for Time Series: An Introduction Using R. Chapman & Hall (London).
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m res <- pseudoRes(m) stats::qqnorm(res$stepRes) stats::qqnorm(res$angleRes)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m res <- pseudoRes(m) stats::qqnorm(res$stepRes) stats::qqnorm(res$angleRes)
Approximate individual-level random effects estimation for state transition probabilities based on Burnham & White (2002)
randomEffects( m, Xformula = ~1, alpha = 0.95, ncores = 1, nlmPar = list(), fit = TRUE, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), modelName = NULL, ... )
randomEffects( m, Xformula = ~1, alpha = 0.95, ncores = 1, nlmPar = list(), fit = TRUE, retryFits = 0, retrySD = NULL, optMethod = "nlm", control = list(), modelName = NULL, ... )
m |
A |
Xformula |
Formula for the design matrix of the random effects model. The default |
alpha |
Significance level of the confidence intervals. Default: 0.95 (i.e. 95% CIs). |
ncores |
number of cores to use for parallel processing |
nlmPar |
List of parameters to pass to the optimization function |
fit |
|
retryFits |
Non-negative integer indicating the number of times to attempt to iteratively fit the model using random perturbations of the current parameter estimates as the
initial values for likelihood optimization. See |
retrySD |
An optional list of scalars or vectors indicating the standard deviation to use for normal perturbations of each working scale parameter when |
optMethod |
The optimization method to be used. See |
control |
A list of control parameters to be passed to |
modelName |
An optional character string providing a name for the fitted model. See |
... |
further arguments passed to or from other methods. Not currently used. |
A randomEffects
model similar to a momentuHMM
object, but including the additional random effect components:
varcomp |
A list of length |
traceG |
The trace of the projection matrix for each random effect. |
Burnham, K.P. and White, G.C. 2002. Evaluation of some random effects methodology applicable to bird ringing data. Journal of Applied Statistics 29: 245-264.
McClintock, B.T. 2021. Worth the effort? A practical examination of random effects in hidden Markov models for animal telemetry data. Methods in Ecology and Evolution doi:10.1111/2041-210X.13619.
## Not run: # simulated data with normal random effects # and binary individual covariate nbAnimals <- 5 # should be larger for random effects estimation obsPerAnimal <- 110 indCov <- rbinom(nbAnimals,1,0.5) # individual covariate betaCov <- c(-0.5,0.5) # covariate effects mu <- c(-0.1,0.1) # mean for random effects sigma <- c(0.2,0.4) # sigma for random effects beta0 <- cbind(rnorm(nbAnimals,mu[1],sigma[1]), rnorm(nbAnimals,mu[2],sigma[2])) reData <- simData(nbAnimals=nbAnimals,obsPerAnimal=obsPerAnimal,nbStates=2, dist=list(step="gamma"),formula=~0+ID+indCov, Par=list(step=c(1,10,1,2)), beta=rbind(beta0,betaCov), covs=data.frame(indCov=rep(indCov,each=obsPerAnimal))) # fit null model nullFit <- fitHMM(reData,nbStates=2, dist=list(step="gamma"), Par0=list(step=c(1,10,1,2))) # fit covariate model covFit <- fitHMM(reData,nbStates=2, dist=list(step="gamma"),formula=~indCov, Par0=list(step=c(1,10,1,2)), beta0=rbind(mu,betaCov)) # fit fixed effects model fixFit <- fitHMM(reData,nbStates=2, dist=list(step="gamma"),formula=~0+ID, Par0=list(step=c(1,10,1,2)), beta0=beta0) # fit random effect model reFit <- randomEffects(fixFit) # fit random effect model with individual covariate reCovFit <- randomEffects(fixFit, Xformula=~indCov) # compare by AICc AIC(nullFit,covFit,fixFit,reFit,reCovFit, n=nrow(reData)) ## End(Not run)
## Not run: # simulated data with normal random effects # and binary individual covariate nbAnimals <- 5 # should be larger for random effects estimation obsPerAnimal <- 110 indCov <- rbinom(nbAnimals,1,0.5) # individual covariate betaCov <- c(-0.5,0.5) # covariate effects mu <- c(-0.1,0.1) # mean for random effects sigma <- c(0.2,0.4) # sigma for random effects beta0 <- cbind(rnorm(nbAnimals,mu[1],sigma[1]), rnorm(nbAnimals,mu[2],sigma[2])) reData <- simData(nbAnimals=nbAnimals,obsPerAnimal=obsPerAnimal,nbStates=2, dist=list(step="gamma"),formula=~0+ID+indCov, Par=list(step=c(1,10,1,2)), beta=rbind(beta0,betaCov), covs=data.frame(indCov=rep(indCov,each=obsPerAnimal))) # fit null model nullFit <- fitHMM(reData,nbStates=2, dist=list(step="gamma"), Par0=list(step=c(1,10,1,2))) # fit covariate model covFit <- fitHMM(reData,nbStates=2, dist=list(step="gamma"),formula=~indCov, Par0=list(step=c(1,10,1,2)), beta0=rbind(mu,betaCov)) # fit fixed effects model fixFit <- fitHMM(reData,nbStates=2, dist=list(step="gamma"),formula=~0+ID, Par0=list(step=c(1,10,1,2)), beta0=beta0) # fit random effect model reFit <- randomEffects(fixFit) # fit random effect model with individual covariate reCovFit <- randomEffects(fixFit, Xformula=~indCov) # compare by AICc AIC(nullFit,covFit,fixFit,reFit,reCovFit, n=nrow(reData)) ## End(Not run)
modelName
for a momentuHMM
, miHMM
, HMMfits
, or miSum
objectSet modelName
for a momentuHMM
, miHMM
, HMMfits
, or miSum
object
setModelName(model, modelName)
setModelName(model, modelName)
model |
|
modelName |
Character string providing a name for the model. See |
model
object with new modelName
field
m <- example$m mName <- setModelName(m, modelName="example")
m <- example$m mName <- setModelName(m, modelName="example")
stateNames
for a momentuHMM
, miHMM
, HMMfits
, or miSum
objectSet stateNames
for a momentuHMM
, miHMM
, HMMfits
, or miSum
object
setStateNames(model, stateNames)
setStateNames(model, stateNames)
model |
|
stateNames |
Character string providing state names for the model. See |
model
object with new stateNames
field
m <- example$m mName <- setStateNames(m, stateNames=c("encamped","exploratory"))
m <- example$m mName <- setStateNames(m, stateNames=c("encamped","exploratory"))
Simulates data from a (multivariate) continuous-time discrete-space hidden Markov model based on Blackwell et al. (2016), where specification of kappa
determines the potential times of state switches when formula
includes covariates (larger kappa
means potential switches are more frequent). Note that state transitions can only occur when there is a cell transition (unless moveState
is set to FALSE
), and time-varying covariates are assumed piece-wise constant between observations. Movement data are assumed to be in Cartesian coordinates (not longitude/latitude) and can be generated with or without observation error attributable to location measurement error.
simCTDS( nbAnimals = 1, nbStates = 2, dist, Par, beta = NULL, delta = NULL, formula = ~1, formulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbCovs = 0, rast, spatialCovs = NULL, spatialCovs.grad = NULL, directions = 4, normalize.gradients = FALSE, grad.point.decreasing = FALSE, zero.idx = integer(), moveState = FALSE, obsPerAnimal = c(500, 1500), initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, betaRef = NULL, stateNames = NULL, model = NULL, matchModelObs = TRUE, states = FALSE, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL, keepSwitch = FALSE, kappa = NULL )
simCTDS( nbAnimals = 1, nbStates = 2, dist, Par, beta = NULL, delta = NULL, formula = ~1, formulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbCovs = 0, rast, spatialCovs = NULL, spatialCovs.grad = NULL, directions = 4, normalize.gradients = FALSE, grad.point.decreasing = FALSE, zero.idx = integer(), moveState = FALSE, obsPerAnimal = c(500, 1500), initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, betaRef = NULL, stateNames = NULL, model = NULL, matchModelObs = TRUE, states = FALSE, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL, keepSwitch = FALSE, kappa = NULL )
nbAnimals |
Number of observed individuals to simulate. |
nbStates |
Number of behavioural states to simulate. |
dist |
A named list indicating the probability distributions of the data streams. A data stream named |
Par |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in If |
beta |
Matrix of regression parameters for the transition rates. |
delta |
Initial value for the initial distribution of the HMM. Default: |
formula |
Regression formula for the transition rate covariates. Default: |
formulaDelta |
Regression formula for the initial distribution. Default: |
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. Default: |
covs |
Covariate values to include in the simulated data, as a dataframe. The names of any covariates specified by |
nbCovs |
Number of covariates to simulate (0 by default). Does not need to be specified if
|
rast |
A raster object or raster stack object that will define the discrete-space grid cells for the CTMC movement path. |
spatialCovs |
List of |
spatialCovs.grad |
List of |
directions |
Integer. Either 4 (indicating a "Rook's neighborhood" of 4 neighboring grid cells) or 8 (indicating a "King's neighborhood" of 8 neighboring grid cells). |
normalize.gradients |
Logical. Default is FALSE. If TRUE, then all gradient covariates for |
grad.point.decreasing |
Logical. If TRUE, then the gradient covariates are positive in the direction of decreasing values of the covariate. If FALSE, then the gradient covariates are positive in the direction of increasing values of the covariate (like a true gradient). |
zero.idx |
Integer vector of the indices of raster cells that are not passable and should be excluded. These are cells where movement should be impossible. Default is zero.idx=integer(). |
moveState |
Logical indicating whether or not transitions out of the current state are forbidden when the animal does not move out of the current cell. Default: |
obsPerAnimal |
Either 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. Alternatively, |
initialPosition |
2-vector providing the x- and y-coordinates of the initial position for all animals. Alternatively, |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
betaRef |
Numeric vector of length |
stateNames |
Optional character vector of length nbStates indicating state names. |
model |
A |
matchModelObs |
If |
states |
|
lambda |
Observation rate. |
errorEllipse |
List providing the upper bound for the semi-major axis ( |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
export |
Character vector of the names of any additional objects or functions in the global environment that are used in |
keepSwitch |
Logical indicating whether or not to return the (typically unobserved) data at the times when potential state switches could have occurred. Default: |
kappa |
List of the form Alternatively, |
simCTDS
assumes the snapshot property applies to all data stream distributions (i.e. observations are "instantaneous") except for the continuous-time discrete-space (ctds
), (multivariate) normal random walk (rw_norm
, rw_mvnorm2
, rw_mvnorm3
), and Poisson (pois
) distributions. For these particular distributions, the observed data are not "instantaneous"; they depend on the time interval between observations and, hence, the state sequence during the entire interval.
If fitting with
fitCTHMM
(or MIfitCTHMM
), it is critical that the frequency of observations (specified by lambda
) is high relative to the serial correlation in the hidden state process (specified by beta
) in order for the discrete-time approximation of fitCTHMM
to be reasonably accurate for these distributions.
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).
When covariates are not included in formulaDelta
(i.e. formulaDelta=NULL
), then delta
is specified as a vector of length nbStates
that
sums to 1. When covariates are included in formulaDelta
, then delta
must be specified
as a k x (nbStates
-1) matrix of working parameters, where k is the number of regression coefficients and the columns correspond to states 2:nbStates
. For example, in a 3-state
HMM with formulaDelta=~cov1+cov2
, the matrix delta
has three rows (intercept + two covariates)
and 2 columns (corresponding to states 2 and 3). The initial distribution working parameters are transformed to the real scale as exp(covsDelta*Delta)/rowSums(exp(covsDelta*Delta))
, where covsDelta
is the N x k design matrix, Delta=cbind(rep(0,k),delta)
is a k x nbStates
matrix of working parameters,
and N=length(unique(data$ID))
.
If the simulated data have no measurement error (i.e., errorEllipse=NULL
), a momentuHMMData
object,
i.e., a dataframe of:
ID |
The ID(s) of the observed animal(s) |
time |
Numeric time of each observed observation |
... |
Data streams as specified by |
x |
Easting location cell centroid |
y |
Northing location cell centroid |
... |
Covariates (if any). If |
If simulated location data include measurement error (i.e., errorEllipse!=NULL
), a dataframe of:
ID |
The ID(s) of the observed animal(s) |
time |
Numeric time of each observed (and missing) observation |
z |
Categorical index indicating observed cell movement (where |
x |
Easting observed location cell centroid |
y |
Norting observed location cell centroid |
tau |
Time difference between consecutive observations |
... |
Data streams that are not derived from location (if applicable) |
... |
Covariates at true ( |
mux |
Easting true location cell centroid |
muy |
Northingtrue location cell centroid |
error_semimajor_axis |
error ellipse semi-major axis (if applicable) |
error_semiminor_axis |
error ellipse semi-minor axis (if applicable) |
error_ellipse_orientation |
error ellipse orientation (if applicable) |
ln.sd.x |
log of the square root of the x-variance of bivariate normal error (if applicable; required for error ellipse models in |
ln.sd.y |
log of the square root of the y-variance of bivariate normal error (if applicable; required for error ellipse models in |
error.corr |
correlation term of bivariate normal error (if applicable; required for error ellipse models in |
Simulates data from a (multivariate) continuous-time hidden Markov model based on Blackwell et al. (2016), where specification of kappa
determines the potential times of state switches when formula
includes covariates (larger kappa
means potential switches are more frequent). Note that any time-varying covariates are assumed piece-wise constant between observations. Movement data are assumed to be in Cartesian coordinates (not longitude/latitude) and can be generated with or without observation error attributable to location measurement error.
simCTHMM( nbAnimals = 1, nbStates = 2, dist, Par, beta = NULL, delta = NULL, formula = ~1, formulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbCovs = 0, spatialCovs = NULL, obsPerAnimal = c(500, 1500), initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, betaRef = NULL, mvnCoords = NULL, stateNames = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, keepSwitch = FALSE, kappa = NULL, TMB = FALSE )
simCTHMM( nbAnimals = 1, nbStates = 2, dist, Par, beta = NULL, delta = NULL, formula = ~1, formulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbCovs = 0, spatialCovs = NULL, obsPerAnimal = c(500, 1500), initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, betaRef = NULL, mvnCoords = NULL, stateNames = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, keepSwitch = FALSE, kappa = NULL, TMB = FALSE )
nbAnimals |
Number of observed individuals to simulate. |
nbStates |
Number of behavioural states to simulate. |
dist |
A named list indicating the probability distributions of the data streams. |
Par |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in If |
beta |
Matrix of regression parameters for the transition rates. |
delta |
Initial distribution of the HMM. Default: |
formula |
Regression formula for the transition rate covariates. Default: |
formulaDelta |
Regression formula for the initial distribution. Default: |
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. Default: |
covs |
Covariate values to include in the simulated data, as a dataframe. The names of any covariates specified by |
nbCovs |
Number of covariates to simulate (0 by default). Does not need to be specified if
|
spatialCovs |
List of |
obsPerAnimal |
Either 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. Alternatively, |
initialPosition |
2-vector providing the x- and y-coordinates of the initial position for all animals. Alternatively, |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
betaRef |
Numeric vector of length |
mvnCoords |
Character string indicating the name of location data that are to be simulated using a multivariate normal distribution. For example, if |
stateNames |
Optional character vector of length nbStates indicating state names. |
model |
A |
matchModelObs |
If |
states |
|
retrySims |
Number of times to attempt to simulate data within the spatial extent of |
lambda |
Observation rate. |
errorEllipse |
List providing the upper bound for the semi-major axis ( |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
export |
Character vector of the names of any additional objects or functions in the global environment that are used in |
gradient |
Logical indicating whether or not to calculate gradients of |
keepSwitch |
Logical indicating whether or not to return the (typically unobserved) data at the times when potential state switches could have occurred. Default: |
kappa |
List of the form Alternatively, If |
TMB |
Logical indicating whether or not data should be simulated based on TMB model specifications. Default: |
simCTHMM
assumes the snapshot property applies to all data stream distributions (i.e. observations are "instantaneous") except for the (multivariate) normal random walk (rw_norm
, rw_mvnorm2
, rw_mvnorm3
) and Poisson (pois
) distributions. For these particular distributions, the observed data are not "instantaneous"; they depend on the time interval between observations and, hence, the state sequence during the entire interval.
If fitting with
fitCTHMM
(or MIfitCTHMM
), it is critical that the frequency of observations (specified by lambda
) is high relative to the serial correlation in the hidden state process (specified by beta
) in order for the discrete-time approximation of fitCTHMM
to be reasonably accurate for these distributions.
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).
When covariates are not included in formulaDelta
(i.e. formulaDelta=NULL
), then delta
is specified as a vector of length nbStates
that
sums to 1. When covariates are included in formulaDelta
, then delta
must be specified
as a k x (nbStates
-1) matrix of working parameters, where k is the number of regression coefficients and the columns correspond to states 2:nbStates
. For example, in a 3-state
HMM with formulaDelta=~cov1+cov2
, the matrix delta
has three rows (intercept + two covariates)
and 2 columns (corresponding to states 2 and 3). The initial distribution working parameters are transformed to the real scale as exp(covsDelta*Delta)/rowSums(exp(covsDelta*Delta))
, where covsDelta
is the N x k design matrix, Delta=cbind(rep(0,k),delta)
is a k x nbStates
matrix of working parameters,
and N=length(unique(data$ID))
.
If the simulated data have no measurement error (i.e., errorEllipse=NULL
), a momentuHMMData
object,
i.e., a dataframe of:
ID |
The ID(s) of the observed animal(s) |
time |
Numeric time of each observed observation |
... |
Data streams as specified by |
x |
Either easting or longitude (if data streams include valid non-negative distribution for 'step') |
y |
Either norting or latitude (if data streams include valid non-negative distribution for 'step') |
... |
Covariates (if any) |
If simulated location data include measurement error (i.e., errorEllipse!=NULL
), a dataframe of:
ID |
The ID(s) of the observed animal(s) |
time |
Numeric time of each observed (and missing) observation |
x |
Either easting or longitude observed location |
y |
Either norting or latitude observed location |
... |
Data streams that are not derived from location (if applicable) |
... |
Covariates at true ( |
mux |
Either easting or longitude true location |
muy |
Either norting or latitude true location |
error_semimajor_axis |
error ellipse semi-major axis (if applicable) |
error_semiminor_axis |
error ellipse semi-minor axis (if applicable) |
error_ellipse_orientation |
error ellipse orientation (if applicable) |
ln.sd.x |
log of the square root of the x-variance of bivariate normal error (if applicable; required for error ellipse models in |
ln.sd.y |
log of the square root of the y-variance of bivariate normal error (if applicable; required for error ellipse models in |
error.corr |
correlation term of bivariate normal error (if applicable; required for error ellipse models in |
Blackwell, P.G., Niu, M., Lambert, M.S. and LaPoint, S.D., 2016. Exact Bayesian inference for animal movement in continuous time. Methods in Ecology and Evolution, 7 (2), 184-195.
Simulates data from a (multivariate) hidden Markov model. Movement data are assumed to be in Cartesian coordinates (not longitude/latitude) and can be generated with or without observation error attributable to temporal irregularity or location measurement error.
simData( nbAnimals = 1, nbStates = 2, dist, Par, beta = NULL, delta = NULL, formula = ~1, formulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbCovs = 0, spatialCovs = NULL, zeroInflation = NULL, oneInflation = NULL, circularAngleMean = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, obsPerAnimal = c(500, 1500), initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, betaRef = NULL, mvnCoords = NULL, stateNames = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = NULL, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, TMB = FALSE, ... ) simHierData( nbAnimals = 1, hierStates, hierDist, Par, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbHierCovs = NULL, spatialCovs = NULL, zeroInflation = NULL, oneInflation = NULL, circularAngleMean = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, obsPerLevel, initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, mvnCoords = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = NULL, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, TMB = FALSE, ... )
simData( nbAnimals = 1, nbStates = 2, dist, Par, beta = NULL, delta = NULL, formula = ~1, formulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbCovs = 0, spatialCovs = NULL, zeroInflation = NULL, oneInflation = NULL, circularAngleMean = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, obsPerAnimal = c(500, 1500), initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, betaRef = NULL, mvnCoords = NULL, stateNames = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = NULL, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, TMB = FALSE, ... ) simHierData( nbAnimals = 1, hierStates, hierDist, Par, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbHierCovs = NULL, spatialCovs = NULL, zeroInflation = NULL, oneInflation = NULL, circularAngleMean = NULL, centers = NULL, centroids = NULL, angleCovs = NULL, obsPerLevel, initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, mvnCoords = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = NULL, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, TMB = FALSE, ... )
nbAnimals |
Number of observed individuals to simulate. |
nbStates |
Number of behavioural states to simulate. |
dist |
A named list indicating the probability distributions of the data streams. Currently
supported distributions are 'bern', 'beta', 'cat', 'exp', 'gamma', 'lnorm', 'logis', 'negbinom', 'norm', 'mvnorm2' (bivariate normal distribution), 'mvnorm3' (trivariate normal distribution),
'pois', 'rw_norm' (normal random walk), 'rw_mvnorm2' (bivariate normal random walk), 'rw_mvnorm3' (trivariate normal random walk), 'vm', 'vmConsensus', 'weibull', and 'wrpcauchy'. For example,
|
Par |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in If |
beta |
Matrix of regression parameters for the transition probabilities (more information in "Details"). |
delta |
Initial distribution of the HMM. Default: |
formula |
Regression formula for the transition probability covariates. Default: |
formulaDelta |
Regression formula for the initial distribution. Default: |
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. Default: |
covs |
Covariate values to include in the simulated data, as a dataframe. The names of any covariates specified by |
nbCovs |
Number of covariates to simulate (0 by default). Does not need to be specified if
|
spatialCovs |
List of |
zeroInflation |
A named list of logicals indicating whether the probability distributions of the data streams should be zero-inflated. If |
oneInflation |
A named list of logicals indicating whether the probability distributions of the data streams should be one-inflated. If |
circularAngleMean |
An optional named list indicating whether to use circular-linear (FALSE) or circular-circular (TRUE)
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. For example,
Alternatively, |
centers |
2-column matrix providing the x-coordinates (column 1) and y-coordinates (column 2) for any activity centers (e.g., potential
centers of attraction or repulsion) from which distance and angle covariates will be calculated based on the simulated location data. These distance and angle
covariates can be included in |
centroids |
List where each element is a data frame consisting of at least |
angleCovs |
Character vector indicating the names of any circular-circular regression angular covariates in |
obsPerAnimal |
Either 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. Alternatively, |
initialPosition |
2-vector providing the x- and y-coordinates of the initial position for all animals. Alternatively, |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of Design matrices specified using formulas allow standard functions in R formulas
(e.g., |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability
distribution parameters for each data stream. For example, for a 2-state model using the wrapped Cauchy ('wrpcauchy') distribution for
a data stream named 'angle' with |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
betaRef |
Numeric vector of length |
mvnCoords |
Character string indicating the name of location data that are to be simulated using a multivariate normal distribution. For example, if |
stateNames |
Optional character vector of length nbStates indicating state names. |
model |
A |
matchModelObs |
If |
states |
|
retrySims |
Number of times to attempt to simulate data within the spatial extent of |
lambda |
Observation rate for location data. If |
errorEllipse |
List providing the upper bound for the semi-major axis ( |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
export |
Character vector of the names of any additional objects or functions in the global environment that are used in |
gradient |
Logical indicating whether or not to calculate gradients of |
TMB |
Logical indicating whether or not data should be simulated based on TMB model specifications. Default: |
... |
further arguments passed to or from other methods |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). Default: |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
nbHierCovs |
A hierarchical data structure |
obsPerLevel |
A hierarchical data structure |
simHierData
is very similar to simData
except that instead of simply specifying the number of states (nbStates
), distributions (dist
), observations (obsPerAnimal
), covariates (nbCovs
), and a single t.p.m. formula (formula
), the hierStates
argument specifies the hierarchical nature of the states,
the hierDist
argument specifies the hierarchical nature of the data streams, the obsPerLevel
argument specifies the number of observations for each level of the hierarchy, the nbHierCovs
argument specifies the number of covariates for each level of the hierarchy, and the hierFormula
argument specifies a t.p.m. formula for each level of the hierarchy.
All of the hierarchical arguments in simHierData
are specified as Node
objects from the data.tree
package.
x- and y-coordinate location data are generated only if valid 'step' and 'angle' data streams are specified. Vaild distributions for 'step' include 'gamma', 'weibull', 'exp', and 'lnorm'. Valid distributions for 'angle' include 'vm' and 'wrpcauchy'. If only a valid 'step' data stream is specified, then only x-coordinates are generated.
If DM
is specified for a particular data stream, then the initial values are specified on
the working (i.e., beta) scale of the parameters. The working scale of each parameter is determined by the link function used.
The function getParDM
is intended to help with obtaining initial values on the working scale when specifying a design matrix and other
parameter constraints.
Simulated data that are temporally regular (i.e., lambda=NULL
) and without location measurement error (i.e., errorEllipse=NULL
) are returned
as a momentuHMMData
(or momentuHierHMMData
) object suitable for analysis using fitHMM
.
Simulated location data that are temporally-irregular (i.e., lambda>0
) and/or with location measurement error (i.e., errorEllipse!=NULL
) are returned
as a data frame suitable for analysis using crawlWrap
.
The matrix 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 formula
covariates, the matrix 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), beta
has one row, for the intercept.
State-specific formulas can be specified in DM
using special formula functions. These special functions can take
the names paste0("state",1:nbStates)
(where the integer indicates the state-specific formula). For example,
DM=list(step=list(mean=~cov1+state1(cov2),sd=~cov2+state2(cov1)))
includes cov1
on the mean parameter for all states, cov2
on the mean parameter for state 1, cov2
on the sd parameter for all states, and cov1
on the sd parameter for state 2.
State- and parameter-specific formulas can be specified for transition probabilities in formula
using special formula functions.
These special functions can take the names paste0("state",1:nbStates)
(where the integer indicates the current state from which transitions occur),
paste0("toState",1:nbStates)
(where the integer indicates the state to which transitions occur),
or paste0("betaCol",nbStates*(nbStates-1))
(where the integer indicates the column of the beta
matrix). For example with nbStates=3
,
formula=~cov1+betaCol1(cov2)+state3(cov3)+toState1(cov4)
includes cov1
on all transition probability parameters, cov2
on the beta
column corresponding
to the transition from state 1->2, cov3
on transition probabilities from state 3 (i.e., beta
columns corresponding to state transitions 3->1 and 3->2),
and cov4
on transition probabilities to state 1 (i.e., beta
columns corresponding to state transitions 2->1 and 3->1).
Cyclical relationships (e.g., hourly, monthly) may be simulated using the consinor(x,period)
special formula function for covariate x
and sine curve period of time length period
. For example, if
the data are hourly, a 24-hour cycle can be simulated using ~cosinor(cov1,24)
, where the covariate cov1
is a repeating series
of integers 0,1,...,23,0,1,...,23,0,1,...
(note that simData
will not do this for you, the appropriate covariate must be specified using the covs
argument; see example below).
The cosinor(x,period)
function converts x
to 2 covariates
cosinorCos(x)=cos(2*pi*x/period)
and consinorSin(x)=sin(2*pi*x/period
for inclusion in the model (i.e., 2 additional parameters per state). The amplitude of the sine wave
is thus sqrt(B_cos^2 + B_sin^2)
, where B_cos
and B_sin
are the working parameters correponding to cosinorCos(x)
and cosinorSin(x)
, respectively (e.g., see Cornelissen 2014).
When the circular-circular regression model is used, the special function angleFormula(cov,strength,by)
can be used in DM
for the mean of angular distributions (i.e. 'vm', 'vmConsensus', and 'wrpcauchy'),
where cov
is an angle covariate (e.g. wind direction), strength
is a positive real covariate (e.g. wind speed), and by
is an optional factor variable for individual- or group-level effects (e.g. ID, sex). This allows angle covariates to be weighted based on their strength or importance at time step t as in
Rivest et al. (2016).
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' (or 'obsPerLevel' for simHierData
), 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).
For simData
, when covariates are not included in formulaDelta
(i.e. formulaDelta=NULL
), then delta
is specified as a vector of length nbStates
that
sums to 1. When covariates are included in formulaDelta
, then delta
must be specified
as a k x (nbStates
-1) matrix of working parameters, where k is the number of regression coefficients and the columns correspond to states 2:nbStates
. For example, in a 3-state
HMM with formulaDelta=~cov1+cov2
, the matrix delta
has three rows (intercept + two covariates)
and 2 columns (corresponding to states 2 and 3). The initial distribution working parameters are transformed to the real scale as exp(covsDelta*Delta)/rowSums(exp(covsDelta*Delta))
, where covsDelta
is the N x k design matrix, Delta=cbind(rep(0,k),delta)
is a k x nbStates
matrix of working parameters,
and N=length(unique(data$ID))
.
For simHierData
, 'delta' in hierDelta
must be specified
as a k x (nbStates
-1) matrix of working parameters, where k is the number of regression coefficients and the columns correspond to states 2:nbStates
.
If the simulated data are temporally regular (i.e., lambda=NULL
) with no measurement error (i.e., errorEllipse=NULL
), an object momentuHMMData
(or momentuHierHMMData
),
i.e., a dataframe of:
ID |
The ID(s) of the observed animal(s) |
... |
Data streams as specified by |
x |
Either easting or longitude (if data streams include valid non-negative distribution for 'step') |
y |
Either northing or latitude (if data streams include valid non-negative distribution for 'step') |
... |
Covariates (if any) |
If simulated location data are temporally irregular (i.e., lambda>0
) and/or include measurement error (i.e., errorEllipse!=NULL
), a dataframe of:
time |
Numeric time of each observed (and missing) observation |
ID |
The ID(s) of the observed animal(s) |
x |
Either easting or longitude observed location |
y |
Either northing or latitude observed location |
... |
Data streams that are not derived from location (if applicable) |
... |
Covariates at temporally-regular true ( |
mux |
Either easting or longitude true location |
muy |
Either northing or latitude true location |
error_semimajor_axis |
error ellipse semi-major axis (if applicable) |
error_semiminor_axis |
error ellipse semi-minor axis (if applicable) |
error_ellipse_orientation |
error ellipse orientation (if applicable) |
ln.sd.x |
log of the square root of the x-variance of bivariate normal error (if applicable; required for error ellipse models in |
ln.sd.y |
log of the square root of the y-variance of bivariate normal error (if applicable; required for error ellipse models in |
error.corr |
correlation term of bivariate normal error (if applicable; required for error ellipse models in |
Cornelissen, G. 2014. Cosinor-based rhythmometry. Theoretical Biology and Medical Modelling 11:16.
McClintock BT, London JM, Cameron MF, Boveng PL. 2015. Modelling animal movement using the Argos satellite telemetry location error ellipse. Methods in Ecology and Evolution 6(3):266-277.
Rivest, LP, Duchesne, T, Nicosia, A, Fortin, D, 2016. A general angular regression model for the analysis of data on animal movement in ecology. Journal of the Royal Statistical Society: Series C (Applied Statistics), 65(3):445-463.
Leos-Barajas, V., Gangloff, E.J., Adam, T., Langrock, R., van Beest, F.M., Nabe-Nielsen, J. and Morales, J.M. 2017. Multi-scale modeling of animal movement and general behavior data using hidden Markov models with hierarchical structures. Journal of Agricultural, Biological and Environmental Statistics, 22 (3), 232-248.
# 1. Pass a fitted model to simulate from # (m is a momentuHMM 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) ## Not run: # 2. Pass the parameters of the model to simulate from stepPar <- c(1,10,1,5,0.2,0.3) # mean_1, mean_2, sd_1, sd_2, zeromass_1, zeromass_2 anglePar <- c(pi,0,0.5,2) # mean_1, mean_2, concentration_1, concentration_2 omegaPar <- c(1,10,10,1) # shape1_1, shape1_2, shape2_1, shape2_2 stepDist <- "gamma" angleDist <- "vm" omegaDist <- "beta" data <- simData(nbAnimals=4,nbStates=2,dist=list(step=stepDist,angle=angleDist,omega=omegaDist), Par=list(step=stepPar,angle=anglePar,omega=omegaPar),nbCovs=2, zeroInflation=list(step=TRUE), obsPerAnimal=obsPerAnimal) # 3. Include covariates # (note that it is useless to specify "nbCovs", which are overruled # by the number of columns of "cov") cov <- data.frame(temp=log(rnorm(500,20,5))) stepPar <- c(log(10),0.1,log(100),-0.1,log(5),log(25)) # working scale parameters for step DM anglePar <- c(pi,0,0.5,2) # mean_1, mean_2, concentration_1, concentration_2 stepDist <- "gamma" angleDist <- "vm" data <- simData(nbAnimals=2,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=list(mean=~temp,sd=~1)), covs=cov, obsPerAnimal=obsPerAnimal) # 4. Include example 'forest' spatial covariate raster layer # nbAnimals and obsPerAnimal kept small to reduce example run time spatialCov<-list(forest=forest) data <- simData(nbAnimals=1,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par=list(step=c(100,1000,50,100),angle=c(0,0,0.1,5)), beta=matrix(c(5,-10,-25,50),nrow=2,ncol=2,byrow=TRUE), formula=~forest,spatialCovs=spatialCov, obsPerAnimal=250,states=TRUE, retrySims=100) # 5. Specify design matrix for 'omega' data stream # natural scale parameters for step and angle stepPar <- c(1,10,1,5) # shape_1, shape_2, scale_1, scale_2 anglePar <- c(pi,0,0.5,0.7) # mean_1, mean_2, concentration_1, concentration_2 # working scale parameters for omega DM omegaPar <- c(log(1),0.1,log(10),-0.1,log(10),-0.1,log(1),0.1) stepDist <- "weibull" angleDist <- "wrpcauchy" omegaDist <- "beta" data <- simData(nbStates=2,dist=list(step=stepDist,angle=angleDist,omega=omegaDist), Par=list(step=stepPar,angle=anglePar,omega=omegaPar),nbCovs=2, DM=list(omega=list(shape1=~cov1,shape2=~cov2)), obsPerAnimal=obsPerAnimal,states=TRUE) # 6. Include temporal irregularity and location measurement error lambda <- 2 # expect 2 observations per time step errorEllipse <- list(M=50,m=25,r=180) obsData <- simData(model=m,obsPerAnimal=obsPerAnimal, lambda=lambda, errorEllipse=errorEllipse) # 7. Cosinor and state-dependent formulas nbStates<-2 dist<-list(step="gamma") Par<-list(step=c(100,1000,50,100)) # include 24-hour cycle on all transition probabilities # include 12-hour cycle on transitions from state 2 formula=~cosinor(hour24,24)+state2(cosinor(hour12,12)) # specify appropriate covariates covs<-data.frame(hour24=0:23,hour12=0:11) beta<-matrix(c(-1.5,1,1,NA,NA,-1.5,-1,-1,1,1),5,2) # row names for beta not required but can be helpful rownames(beta)<-c("(Intercept)", "cosinorCos(hour24, 24)", "cosinorSin(hour24, 24)", "cosinorCos(hour12, 12)", "cosinorSin(hour12, 12)") data.cos<-simData(nbStates=nbStates,dist=dist,Par=Par, beta=beta,formula=formula,covs=covs) # 8. Piecewise constant B-spline on step length mean and angle concentration nObs <- 1000 # length of simulated track cov <- data.frame(time=1:nObs) # time covariate for splines dist <- list(step="gamma",angle="vm") stepDM <- list(mean=~splines2::bSpline(time,df=2,degree=0),sd=~1) angleDM <- list(mean=~1,concentration=~splines2::bSpline(time,df=2,degree=0)) DM <- list(step=stepDM,angle=angleDM) Par <- list(step=c(log(1000),1,-1,log(100)),angle=c(0,log(10),2,-5)) data.spline<-simData(obsPerAnimal=nObs,nbStates=1,dist=dist,Par=Par,DM=DM,covs=cov) # 9. Initial state (delta) based on covariate nObs <- 100 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.01,0.75)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs))) # sex covariate formulaDelta <- ~ sex + 0 # Female begins in state 1, male begins in state 2 delta <- matrix(c(-100,100),2,1,dimnames=list(c("sexF","sexM"),"state 2")) data.delta<-simData(nbAnimals=2,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, delta=delta,formulaDelta=formulaDelta,covs=cov, beta=matrix(-1.5,1,2),states=TRUE) ## End(Not run)
# 1. Pass a fitted model to simulate from # (m is a momentuHMM 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) ## Not run: # 2. Pass the parameters of the model to simulate from stepPar <- c(1,10,1,5,0.2,0.3) # mean_1, mean_2, sd_1, sd_2, zeromass_1, zeromass_2 anglePar <- c(pi,0,0.5,2) # mean_1, mean_2, concentration_1, concentration_2 omegaPar <- c(1,10,10,1) # shape1_1, shape1_2, shape2_1, shape2_2 stepDist <- "gamma" angleDist <- "vm" omegaDist <- "beta" data <- simData(nbAnimals=4,nbStates=2,dist=list(step=stepDist,angle=angleDist,omega=omegaDist), Par=list(step=stepPar,angle=anglePar,omega=omegaPar),nbCovs=2, zeroInflation=list(step=TRUE), obsPerAnimal=obsPerAnimal) # 3. Include covariates # (note that it is useless to specify "nbCovs", which are overruled # by the number of columns of "cov") cov <- data.frame(temp=log(rnorm(500,20,5))) stepPar <- c(log(10),0.1,log(100),-0.1,log(5),log(25)) # working scale parameters for step DM anglePar <- c(pi,0,0.5,2) # mean_1, mean_2, concentration_1, concentration_2 stepDist <- "gamma" angleDist <- "vm" data <- simData(nbAnimals=2,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par=list(step=stepPar,angle=anglePar), DM=list(step=list(mean=~temp,sd=~1)), covs=cov, obsPerAnimal=obsPerAnimal) # 4. Include example 'forest' spatial covariate raster layer # nbAnimals and obsPerAnimal kept small to reduce example run time spatialCov<-list(forest=forest) data <- simData(nbAnimals=1,nbStates=2,dist=list(step=stepDist,angle=angleDist), Par=list(step=c(100,1000,50,100),angle=c(0,0,0.1,5)), beta=matrix(c(5,-10,-25,50),nrow=2,ncol=2,byrow=TRUE), formula=~forest,spatialCovs=spatialCov, obsPerAnimal=250,states=TRUE, retrySims=100) # 5. Specify design matrix for 'omega' data stream # natural scale parameters for step and angle stepPar <- c(1,10,1,5) # shape_1, shape_2, scale_1, scale_2 anglePar <- c(pi,0,0.5,0.7) # mean_1, mean_2, concentration_1, concentration_2 # working scale parameters for omega DM omegaPar <- c(log(1),0.1,log(10),-0.1,log(10),-0.1,log(1),0.1) stepDist <- "weibull" angleDist <- "wrpcauchy" omegaDist <- "beta" data <- simData(nbStates=2,dist=list(step=stepDist,angle=angleDist,omega=omegaDist), Par=list(step=stepPar,angle=anglePar,omega=omegaPar),nbCovs=2, DM=list(omega=list(shape1=~cov1,shape2=~cov2)), obsPerAnimal=obsPerAnimal,states=TRUE) # 6. Include temporal irregularity and location measurement error lambda <- 2 # expect 2 observations per time step errorEllipse <- list(M=50,m=25,r=180) obsData <- simData(model=m,obsPerAnimal=obsPerAnimal, lambda=lambda, errorEllipse=errorEllipse) # 7. Cosinor and state-dependent formulas nbStates<-2 dist<-list(step="gamma") Par<-list(step=c(100,1000,50,100)) # include 24-hour cycle on all transition probabilities # include 12-hour cycle on transitions from state 2 formula=~cosinor(hour24,24)+state2(cosinor(hour12,12)) # specify appropriate covariates covs<-data.frame(hour24=0:23,hour12=0:11) beta<-matrix(c(-1.5,1,1,NA,NA,-1.5,-1,-1,1,1),5,2) # row names for beta not required but can be helpful rownames(beta)<-c("(Intercept)", "cosinorCos(hour24, 24)", "cosinorSin(hour24, 24)", "cosinorCos(hour12, 12)", "cosinorSin(hour12, 12)") data.cos<-simData(nbStates=nbStates,dist=dist,Par=Par, beta=beta,formula=formula,covs=covs) # 8. Piecewise constant B-spline on step length mean and angle concentration nObs <- 1000 # length of simulated track cov <- data.frame(time=1:nObs) # time covariate for splines dist <- list(step="gamma",angle="vm") stepDM <- list(mean=~splines2::bSpline(time,df=2,degree=0),sd=~1) angleDM <- list(mean=~1,concentration=~splines2::bSpline(time,df=2,degree=0)) DM <- list(step=stepDM,angle=angleDM) Par <- list(step=c(log(1000),1,-1,log(100)),angle=c(0,log(10),2,-5)) data.spline<-simData(obsPerAnimal=nObs,nbStates=1,dist=dist,Par=Par,DM=DM,covs=cov) # 9. Initial state (delta) based on covariate nObs <- 100 dist <- list(step="gamma",angle="vm") Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.01,0.75)) # create sex covariate cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs))) # sex covariate formulaDelta <- ~ sex + 0 # Female begins in state 1, male begins in state 2 delta <- matrix(c(-100,100),2,1,dimnames=list(c("sexF","sexM"),"state 2")) data.delta<-simData(nbAnimals=2,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par, delta=delta,formulaDelta=formulaDelta,covs=cov, beta=matrix(-1.5,1,2),states=TRUE) ## End(Not run)
Simulates data from an approximate (multivariate) hierarchical continuous-time discrete-space hidden Markov model. Note that the state active at observation time t determines the (state-dependent) observation distribution from observation time t to time t+1 (see Details) and time-varying covariates are assumed piece-wise constant between observations. Movement data are assumed to be in Cartesian coordinates (not longitude/latitude) and can be generated with or without observation error attributable to location measurement error.
simHierCTDS( nbAnimals = 1, hierStates, hierDist, Par, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbHierCovs = NULL, rast, spatialCovs = NULL, spatialCovs.grad = NULL, directions = 4, normalize.gradients = FALSE, grad.point.decreasing = FALSE, zero.idx = integer(), moveState = FALSE, obsPerLevel, initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, model = NULL, matchModelObs = TRUE, states = FALSE, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL )
simHierCTDS( nbAnimals = 1, hierStates, hierDist, Par, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbHierCovs = NULL, rast, spatialCovs = NULL, spatialCovs.grad = NULL, directions = 4, normalize.gradients = FALSE, grad.point.decreasing = FALSE, zero.idx = integer(), moveState = FALSE, obsPerLevel, initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, model = NULL, matchModelObs = TRUE, states = FALSE, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL )
nbAnimals |
Number of observed individuals to simulate. |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
Par |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in If |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). Default: |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. Default: |
covs |
Covariate values to include in the simulated data, as a dataframe. The names of any covariates specified by |
nbHierCovs |
A hierarchical data structure |
rast |
A raster object or raster stack object that will define the discrete-space grid cells for the CTMC movement path. |
spatialCovs |
List of |
spatialCovs.grad |
List of |
directions |
Integer. Either 4 (indicating a "Rook's neighborhood" of 4 neighboring grid cells) or 8 (indicating a "King's neighborhood" of 8 neighboring grid cells). |
normalize.gradients |
Logical. Default is FALSE. If TRUE, then all gradient covariates for |
grad.point.decreasing |
Logical. If TRUE, then the gradient covariates are positive in the direction of decreasing values of the covariate. If FALSE, then the gradient covariates are positive in the direction of increasing values of the covariate (like a true gradient). |
zero.idx |
Integer vector of the indices of raster cells that are not passable and should be excluded. These are cells where movement should be impossible. Default is zero.idx=integer(). |
moveState |
Logical indicating whether or not transitions out of the current state are forbidden when the animal does not move out of the current cell. Default: |
obsPerLevel |
A hierarchical data structure |
initialPosition |
2-vector providing the x- and y-coordinates of the initial position for all animals. Alternatively, |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
model |
A |
matchModelObs |
If |
states |
|
lambda |
Observation rate. |
errorEllipse |
List providing the upper bound for the semi-major axis ( |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
export |
Character vector of the names of any additional objects or functions in the global environment that are used in |
Importantly, simHierCTDS
assumes the state active at observation time t determines the (state-dependent) observation distribution from observation time t to time t+1 (i.e. state switches only occur at the times of observations).
The snapshot property applies to all data stream distributions (i.e. observations are "instantaneous") except for the continuous-time discrete space (ctds
), (multivariate) normal random walk (rw_norm
, rw_mvnorm2
, rw_mvnorm3
), and Poisson (pois
) distributions. For these particular distributions, the observed data are not "instantaneous"; they depend on the time interval between observations .
It is critical that the frequency of observations (specified by
lambda
) is high relative to the serial correlation in the hidden state process (specified by hierBeta
) in order for this discrete-time approximation to be reasonably accurate.
Model specification in simHierCTDS
is similar to simCTDS
except that instead of simply specifying the number of states (nbStates
), distributions (dist
), observations (obsPerAnimal
), covariates (nbCovs
), and a single t.p.m. formula (formula
), the hierStates
argument specifies the hierarchical nature of the states,
the hierDist
argument specifies the hierarchical nature of the data streams, the obsPerLevel
argument specifies the number of observations for each level of the hierarchy, the nbHierCovs
argument specifies the number of covariates for each level of the hierarchy, and the hierFormula
argument specifies a t.p.m. formula for each level of the hierarchy.
All of the hierarchical arguments in simHierCTDS
are specified as Node
objects from the data.tree
package.
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 'obsPerLevel', 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).
In hierDelta
, 'delta' must be specified
as a k x (nbStates
-1) matrix of working parameters, where k is the number of regression coefficients and the columns correspond to states 2:nbStates
.
Simulates data from an approximate (multivariate) hierarchical continuous-time hidden Markov model. Note that the state active at observation time t determines the (state-dependent) observation distribution from observation time t to time t+1 (see Details) and time-varying covariates are assumed piece-wise constant between observations. Movement data are assumed to be in Cartesian coordinates (not longitude/latitude) and can be generated with or without observation error attributable to location measurement error.
simHierCTHMM( nbAnimals = 1, hierStates, hierDist, Par, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbHierCovs = NULL, spatialCovs = NULL, obsPerLevel, initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, mvnCoords = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, TMB = FALSE )
simHierCTHMM( nbAnimals = 1, hierStates, hierDist, Par, hierBeta = NULL, hierDelta = NULL, hierFormula = NULL, hierFormulaDelta = NULL, mixtures = 1, formulaPi = NULL, covs = NULL, nbHierCovs = NULL, spatialCovs = NULL, obsPerLevel, initialPosition = c(0, 0), DM = NULL, userBounds = NULL, workBounds = NULL, mvnCoords = NULL, model = NULL, matchModelObs = FALSE, states = FALSE, retrySims = 0, lambda = 1, errorEllipse = NULL, ncores = 1, export = NULL, gradient = FALSE, TMB = FALSE )
nbAnimals |
Number of observed individuals to simulate. |
hierStates |
A hierarchical model structure |
hierDist |
A hierarchical data structure |
Par |
A named list containing vectors of initial state-dependent probability distribution parameters for
each data stream specified in If |
hierBeta |
A hierarchical data structure |
hierDelta |
A hierarchical data structure |
hierFormula |
A hierarchical formula structure for the transition probability covariates for each level of the hierarchy ('formula'). Default: |
hierFormulaDelta |
A hierarchical formula structure for the initial distribution covariates for each level of the hierarchy ('formulaDelta'). Default: |
mixtures |
Number of mixtures for the state transition probabilities (i.e. discrete random effects *sensu* DeRuiter et al. 2017). Default: |
formulaPi |
Regression formula for the mixture distribution probabilities. Default: |
covs |
Covariate values to include in the simulated data, as a dataframe. The names of any covariates specified by |
nbHierCovs |
A hierarchical data structure |
spatialCovs |
List of |
obsPerLevel |
A hierarchical data structure |
initialPosition |
2-vector providing the x- and y-coordinates of the initial position for all animals. Alternatively, |
DM |
An optional named list indicating the design matrices to be used for the probability distribution parameters of each data
stream. Each element of |
userBounds |
An optional named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
workBounds |
An optional named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters. For each matrix, the first column pertains to the lower bound and the second column the upper bound.
For data streams, each element of |
mvnCoords |
Character string indicating the name of location data that are to be simulated using a multivariate normal distribution. For example, if |
model |
A |
matchModelObs |
If |
states |
|
retrySims |
Number of times to attempt to simulate data within the spatial extent of |
lambda |
Observation rate. |
errorEllipse |
List providing the upper bound for the semi-major axis ( |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). |
export |
Character vector of the names of any additional objects or functions in the global environment that are used in |
gradient |
Logical indicating whether or not to calculate gradients of |
TMB |
Logical indicating whether or not data should be simulated based on TMB model specifications. Default: |
Importantly, simHierCTHMM
assumes the state active at observation time t determines the (state-dependent) observation distribution from observation time t to time t+1 (i.e. state switches only occur at the times of observations).
The snapshot property applies to all data stream distributions (i.e. observations are "instantaneous") except for the (multivariate) normal random walk (rw_norm
, rw_mvnorm2
, rw_mvnorm3
) and Poisson (pois
) distributions. For these particular distributions, the observed data are not "instantaneous"; they depend on the time interval between observations .
It is critical that the frequency of observations (specified by
lambda
) is high relative to the serial correlation in the hidden state process (specified by hierBeta
) in order for this discrete-time approximation to be reasonably accurate.
Model specification in simHierCTHMM
is similar to simCTHMM
except that instead of simply specifying the number of states (nbStates
), distributions (dist
), observations (obsPerAnimal
), covariates (nbCovs
), and a single t.p.m. formula (formula
), the hierStates
argument specifies the hierarchical nature of the states,
the hierDist
argument specifies the hierarchical nature of the data streams, the obsPerLevel
argument specifies the number of observations for each level of the hierarchy, the nbHierCovs
argument specifies the number of covariates for each level of the hierarchy, and the hierFormula
argument specifies a t.p.m. formula for each level of the hierarchy.
All of the hierarchical arguments in simHierCTHMM
are specified as Node
objects from the data.tree
package.
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 'obsPerLevel', 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).
In hierDelta
, 'delta' must be specified
as a k x (nbStates
-1) matrix of working parameters, where k is the number of regression coefficients and the columns correspond to states 2:nbStates
.
Leos-Barajas, V., Gangloff, E.J., Adam, T., Langrock, R., van Beest, F.M., Nabe-Nielsen, J. and Morales, J.M. 2017. Multi-scale modeling of animal movement and general behavior data using hidden Markov models with hierarchical structures. Journal of Agricultural, Biological and Environmental Statistics, 22 (3), 232-248.
Simulates observed location data subject to temporal irregularity and/or location measurement error
simObsData(data, lambda, errorEllipse, ...) ## S3 method for class 'momentuHMMData' simObsData(data, lambda, errorEllipse, ...) ## S3 method for class 'momentuHierHMMData' simObsData(data, lambda, errorEllipse, coordLevel, ...)
simObsData(data, lambda, errorEllipse, ...) ## S3 method for class 'momentuHMMData' simObsData(data, lambda, errorEllipse, ...) ## S3 method for class 'momentuHierHMMData' simObsData(data, lambda, errorEllipse, coordLevel, ...)
data |
A |
lambda |
Observation rate for location data. If |
errorEllipse |
List providing the bounds for the semi-major axis ( |
... |
further arguments passed to or from other methods |
coordLevel |
Level of the hierarchy in which the location data are obtained |
Simulated location data that are temporally-irregular (i.e., lambda>0
) and/or with location measurement error (i.e., errorEllipse!=NULL
) are returned
as a data frame suitable for analysis using crawlWrap
.
A dataframe of:
time |
Numeric time of each observed (and missing) observation |
ID |
The ID(s) of the observed animal(s) |
x |
Either easting or longitude observed location |
y |
Either norting or latitude observed location |
... |
Data streams that are not derived from location (if applicable) |
... |
Covariates at temporally-regular true ( |
mux |
Either easting or longitude true location |
muy |
Either norting or latitude true location |
error_semimajor_axis |
error ellipse semi-major axis (if applicable) |
error_semiminor_axis |
error ellipse semi-minor axis (if applicable) |
error_ellipse_orientation |
error ellipse orientation (if applicable) |
ln.sd.x |
log of the square root of the x-variance of bivariate normal error (if applicable; required for error ellipse models in |
ln.sd.y |
log of the square root of the y-variance of bivariate normal error (if applicable; required for error ellipse models in |
error.corr |
correlation term of bivariate normal error (if applicable; required for error ellipse models in |
McClintock BT, London JM, Cameron MF, Boveng PL. 2015. Modelling animal movement using the Argos satellite telemetry location error ellipse. Methods in Ecology and Evolution 6(3):266-277.
# extract momentuHMMData example data <- example$m$data lambda <- 2 # expect 2 observations per time step errorEllipse <- list(M=c(0,50),m=c(0,50),r=c(0,180)) obsData1 <- simObsData(data,lambda=lambda,errorEllipse=errorEllipse) errorEllipse <- list(M=50,m=50,r=180) obsData2 <- simObsData(data,lambda=lambda,errorEllipse=errorEllipse)
# extract momentuHMMData example data <- example$m$data lambda <- 2 # expect 2 observations per time step errorEllipse <- list(M=c(0,50),m=c(0,50),r=c(0,180)) obsData1 <- simObsData(data,lambda=lambda,errorEllipse=errorEllipse) errorEllipse <- list(M=50,m=50,r=180) obsData2 <- simObsData(data,lambda=lambda,errorEllipse=errorEllipse)
For a given model, computes the probability of the process being in the different states at each time point.
stateProbs(m, hierarchical = FALSE)
stateProbs(m, hierarchical = FALSE)
m |
A |
hierarchical |
Logical indicating whether or not to return a list of state probabilities for each level of a hierarchical HMM. Ignored unless |
The matrix of state probabilities, with element [i,j] the probability of being in state j in observation i.
Zucchini, W. and MacDonald, I.L. 2009. Hidden Markov Models for Time Series: An Introduction Using R. Chapman & Hall (London).
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m sp <- stateProbs(m)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m sp <- stateProbs(m)
Calculates the stationary probabilities of each state based on covariate values.
stationary(model, covs, covIndex)
stationary(model, covs, covIndex)
model |
|
covs |
Either a data frame or a design matrix of covariates. If |
covIndex |
Integer vector indicating specific rows of the data to be used in the calculations. This can be useful for reducing unnecessarily long computation times, e.g., when |
A list of length model$conditions$mixtures
where each element is a matrix of stationary state probabilities for each mixture. For each matrix, each row corresponds to
a row of covs, and each column corresponds to a state.
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m # data frame of covariates stationary(m, covs = data.frame(cov1 = 0, cov2 = 0)) # design matrix (each column corresponds to row of m$mle$beta) stationary(m, covs = matrix(c(1,0,cos(0)),1,3)) # get stationary distribution for first 3 observations stationary(m, covIndex = c(1,2,3))
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m # data frame of covariates stationary(m, covs = data.frame(cov1 = 0, cov2 = 0)) # design matrix (each column corresponds to row of m$mle$beta) stationary(m, covs = matrix(c(1,0,cos(0)),1,3)) # get stationary distribution for first 3 observations stationary(m, covIndex = c(1,2,3))
Written in C++.
stationary_rcpp(A)
stationary_rcpp(A)
A |
transition rate matrix (of dimension nbStates x nbStates) |
row vector of stationary distribution probabilities
momentuHMMData
Summary momentuHMMData
## S3 method for class 'momentuHMMData' summary(object, dataNames = c("step", "angle"), animals = NULL, ...) ## S3 method for class 'momentuHierHMMData' summary(object, dataNames = c("step", "angle", "level"), animals = NULL, ...)
## S3 method for class 'momentuHMMData' summary(object, dataNames = c("step", "angle"), animals = NULL, ...) ## S3 method for class 'momentuHierHMMData' summary(object, dataNames = c("step", "angle", "level"), animals = NULL, ...)
object |
A |
dataNames |
Names of the variables to summarize. Default is |
animals |
Vector of indices or IDs of animals for which data will be summarized.
Default: |
... |
Currently unused. For compatibility with generic method. |
# data is a momentuHMMData object (as returned by prepData), automatically loaded with the package data <- example$m$data summary(data,dataNames=c("step","angle","cov1","cov2"))
# data is a momentuHMMData object (as returned by prepData), automatically loaded with the package data <- example$m$data summary(data,dataNames=c("step","angle","cov1","cov2"))
Calculate proportion of time steps assigned to each state (i.e. “activity budgets”)
timeInStates(m, by = NULL, alpha = 0.95, ncores = 1) ## S3 method for class 'momentuHMM' timeInStates(m, by = NULL, alpha = 0.95, ncores = 1) ## S3 method for class 'HMMfits' timeInStates(m, by = NULL, alpha = 0.95, ncores = 1) ## S3 method for class 'miHMM' timeInStates(m, by = NULL, alpha = 0.95, ncores = 1)
timeInStates(m, by = NULL, alpha = 0.95, ncores = 1) ## S3 method for class 'momentuHMM' timeInStates(m, by = NULL, alpha = 0.95, ncores = 1) ## S3 method for class 'HMMfits' timeInStates(m, by = NULL, alpha = 0.95, ncores = 1) ## S3 method for class 'miHMM' timeInStates(m, by = NULL, alpha = 0.95, ncores = 1)
m |
A |
by |
A character vector indicating any groupings by which to calculate the proportions, such as individual (“ID”) or group-level (e.g. sex or age class) covariates. Default is |
alpha |
Significance level for calculating confidence intervals of pooled estimates. Default: 0.95. Ignored unless |
ncores |
Number of cores to use for parallel processing. Default: 1 (no parallel processing). Ignored unless |
If m
is a continuous-time HMM, then the time in states is calculated based on the actual amount of time that was spent in each state (based on m$data$dt
) instead of the proportion of time steps assigned to each state.
If m
is a momentuHMM
object, a data frame containing the estimated activity budgets for each state (grouped according to by
). If m
is a miHMM
or HMMfits
object, a list containing the activity budget
estimates, standard errors, lower bounds, and upper bounds across all imputations.
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m timeInStates(m) timeInStates(m, by = "ID")
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m timeInStates(m) timeInStates(m, by = "ID")
Computation of the transition probability matrix, as a function of the covariates and the regression
parameters. Written in C++. Used in viterbi
.
trMatrix_rcpp( nbStates, beta, covs, betaRef, CT = FALSE, dt = as.numeric(c()), aInd = as.integer(c()), rateMatrix = FALSE, kappa = NA_real_, check = TRUE )
trMatrix_rcpp( nbStates, beta, covs, betaRef, CT = FALSE, dt = as.numeric(c()), aInd = as.integer(c()), rateMatrix = FALSE, kappa = NA_real_, check = TRUE )
nbStates |
Number of states |
beta |
Matrix of regression parameters |
covs |
Matrix of covariate values |
betaRef |
Indices of reference elements for t.p.m. multinomial logit link. |
CT |
logical indicating discrete-time approximation of a continuous-time model |
dt |
numeric vector of length |
aInd |
Vector of indices of the rows at which the data (i.e. |
rateMatrix |
logical indicating whether or not to return the transition rate matrix. Ignored unless |
kappa |
maximum allowed value for the row sums of the off-diagonal elements in the state transition rate matrix, such that the minimum value for the diagonal elements is |
check |
logical indicating whether or not to check transition probability matrix for issues. Ignored unless |
Three dimensional array trMat
, such that trMat[,,t]
is the transition matrix at
time t.
turnAngle(x, y, z, type = "UTM", angleCov = FALSE)
turnAngle(x, y, z, type = "UTM", angleCov = FALSE)
x |
First point |
y |
Second point |
z |
Third point |
type |
|
angleCov |
logical indicating to not return NA when x=y or y=z. Default: FALSE (i.e. NA is returned if x=y or y=z). |
The angle between vectors (x,y) and (y,z).
If type='LL'
then turning angle is calculated based on initial bearings using bearing
.
## Not run: x <- c(0,0) y <- c(4,6) z <- c(10,7) momentuHMM:::turnAngle(x,y,z) ## End(Not run)
## Not run: x <- c(0,0) y <- c(4,6) z <- c(10,7) momentuHMM:::turnAngle(x,y,z) ## End(Not run)
For a given model, reconstructs the most probable states sequence, using the Viterbi algorithm.
viterbi(m, hierarchical = FALSE)
viterbi(m, hierarchical = FALSE)
m |
An object |
hierarchical |
Logical indicating whether or not to return a list of Viterbi-decoded states for each level of a hierarchical HMM. Ignored unless |
The sequence of most probable states. If hierarchical
is TRUE
, then a list of the most probable states for each level of the hierarchy is returned.
Zucchini, W. and MacDonald, I.L. 2009. Hidden Markov Models for Time Series: An Introduction Using R. Chapman & Hall (London).
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m # reconstruction of states sequence states <- viterbi(m)
# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package m <- example$m # reconstruction of states sequence states <- viterbi(m)
Scales each parameter from the set of real numbers, back to its natural interval. Used during the optimization of the log-likelihood.
w2n( wpar, bounds, parSize, nbStates, nbCovs, estAngleMean, circularAngleMean, consensus, stationary, fullDM, DMind, nbObs, dist, Bndind, nc, meanind, covsDelta, workBounds, covsPi, TMB = FALSE )
w2n( wpar, bounds, parSize, nbStates, nbCovs, estAngleMean, circularAngleMean, consensus, stationary, fullDM, DMind, nbObs, dist, Bndind, nc, meanind, covsDelta, workBounds, covsPi, TMB = FALSE )
wpar |
Vector of working parameters. |
bounds |
Named list of 2-column matrices specifying bounds on the natural (i.e, real) scale of the probability distribution parameters for each data stream. |
parSize |
Named list indicating the number of natural parameters of the data stream probability distributions |
nbStates |
The number of states of the HMM. |
nbCovs |
The number of beta covariates. |
estAngleMean |
Named list indicating whether or not to estimate the angle mean for data streams with angular distributions ('vm' and 'wrpcauchy'). |
circularAngleMean |
Named list indicating whether to use circular-linear or circular-circular
regression on the mean of circular distributions ('vm' and 'wrpcauchy') for turning angles. See |
consensus |
Named list indicating whether to use the circular-circular regression consensus model |
stationary |
|
fullDM |
Named list containing the full (i.e. not shorthand) design matrix for each data stream. |
DMind |
Named list indicating whether |
nbObs |
Number of observations in the data. |
dist |
Named list indicating the probability distributions of the data streams. |
Bndind |
Named list indicating whether |
nc |
indicator for zeros in fullDM |
meanind |
index for circular-circular regression mean angles with at least one non-zero entry in fullDM |
covsDelta |
data frame containing the delta model covariates |
workBounds |
named list of 2-column matrices specifying bounds on the working scale of the probability distribution, transition probability, and initial distribution parameters |
covsPi |
data frame containing the pi model covariates |
TMB |
logical indicating whether or not |
A list of:
... |
Matrices containing the natural parameters for each data stream (e.g., 'step', 'angle', etc.) |
beta |
Matrix of regression coefficients of the transition probabilities |
delta |
Initial distribution |
## Not run: m<-example$m nbStates <- 2 nbCovs <- 2 parSize <- list(step=2,angle=2) par <- list(step=c(t(m$mle$step)),angle=c(t(m$mle$angle))) bounds <- m$conditions$bounds beta <- matrix(rnorm(6),ncol=2,nrow=3) delta <- c(0.6,0.4) #working parameters wpar <- momentuHMM:::n2w(par,bounds,list(beta=beta),log(delta[-1]/delta[1]),nbStates, m$conditions$estAngleMean,NULL,m$conditions$Bndind, m$conditions$dist) #natural parameter p <- momentuHMM:::w2n(wpar,bounds,parSize,nbStates,nbCovs,m$conditions$estAngleMean, m$conditions$circularAngleMean,lapply(m$conditions$dist,function(x) x=="vmConsensus"), m$conditions$stationary,m$conditions$fullDM, m$conditions$DMind,1,m$conditions$dist,m$conditions$Bndind, matrix(1,nrow=length(unique(m$data$ID)),ncol=1),covsDelta=m$covsDelta, workBounds=m$conditions$workBounds) ## End(Not run)
## Not run: m<-example$m nbStates <- 2 nbCovs <- 2 parSize <- list(step=2,angle=2) par <- list(step=c(t(m$mle$step)),angle=c(t(m$mle$angle))) bounds <- m$conditions$bounds beta <- matrix(rnorm(6),ncol=2,nrow=3) delta <- c(0.6,0.4) #working parameters wpar <- momentuHMM:::n2w(par,bounds,list(beta=beta),log(delta[-1]/delta[1]),nbStates, m$conditions$estAngleMean,NULL,m$conditions$Bndind, m$conditions$dist) #natural parameter p <- momentuHMM:::w2n(wpar,bounds,parSize,nbStates,nbCovs,m$conditions$estAngleMean, m$conditions$circularAngleMean,lapply(m$conditions$dist,function(x) x=="vmConsensus"), m$conditions$stationary,m$conditions$fullDM, m$conditions$DMind,1,m$conditions$dist,m$conditions$Bndind, matrix(1,nrow=length(unique(m$data$ID)),ncol=1),covsDelta=m$covsDelta, workBounds=m$conditions$workBounds) ## End(Not run)
Loop for computation of design matrix (X) times the working scale parameters (B). Written in C++. Used in w2n
.
XBloop_rcpp( DM, Xvec, nbObs, nr, nc, circularAngleMean, consensus, rindex, cindex, nbStates, refCoeff = 1 )
XBloop_rcpp( DM, Xvec, nbObs, nr, nc, circularAngleMean, consensus, rindex, cindex, nbStates, refCoeff = 1 )
DM |
design matrix |
Xvec |
working parameters |
nbObs |
number of observations |
nr |
number of rows in design matrix |
nc |
number of column in design matrix |
circularAngleMean |
indicator for whether or not circular-circular regression model |
consensus |
indicator for whether or not circular-circular regression consensus model |
rindex |
row index for design matrix |
cindex |
column index for design matrix |
nbStates |
number of states |
refCoeff |
intercept coefficient for circular-circular regression model |
XB matrix