AIC

Description

Akaike information criterion of momentuHMM model(s).

Usage

## S3 method for class 'momentuHMM'
AIC(object, ..., k = 2, n = NULL)

Arguments

Value

The AIC of the model(s) provided. If several models are provided, the AICs are output in ascending order.

Examples

# 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

Description

Calculate Akaike information criterion model weights

Usage

AICweights(..., k = 2, n = NULL)

Arguments

Details

• 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.

Value

The AIC weights of the models. If multiple imputation objects are provided, then the mean model weights (and standard deviations) are provided.

Examples

## 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)

Confidence intervals for working (i.e., beta) parameters

Description

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:

Usage

CIbeta(m, alpha = 0.95)

Arguments

Details

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.

Value

A list of the following objects:

Examples

# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
m <- example$m

CIbeta(m)

Confidence intervals for the natural (i.e., real) parameters

Description

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.

Usage

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)

Arguments

Details

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).

Value

A list of the following objects:

Examples

# 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)

Constructor of HMMfits objects

Description

Constructor of HMMfits objects

Usage

HMMfits(m)

Arguments

Value

An object HMMfits.

Fit HMMs to multiple imputation data

Description

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.

Usage

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,
  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,
  ...
)

## 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,
  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,
  ...
)

Arguments

Details

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.

Value

If nSims>1, poolEstimates=TRUE, and fit=TRUE, a miHMM object, i.e., a list consisting of:

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.

References

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.

See Also

crawlWrap, crwPostIS, crwSimulator, fitHMM, getParDM, MIpool, prepData

Examples

# 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)

Calculate pooled parameter estimates and states across multiple imputations

Description

Calculate pooled parameter estimates and states across multiple imputations

Usage

MIpool(im, alpha = 0.95, ncores = 1, covs = NULL, na.rm = FALSE)

Arguments

Details

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.

Value

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:

Examples

## 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)

Get XB

Description

Loop for computation of design matrix (X) times the working scale parameters (B). Written in C++. Used in w2n.

Usage

XBloop_rcpp(
  DM,
  Xvec,
  nbObs,
  nr,
  nc,
  circularAngleMean,
  consensus,
  rindex,
  cindex,
  nbStates,
  refCoeff = 1
)

Arguments

Value

XB matrix

Matrix of all probabilities

Description

Used in functions viterbi, logAlpha, logBeta.

Usage

allProbs(m)

Arguments

Value

Matrix of all probabilities.

Examples

## Not run: 
P <- momentuHMM:::allProbs(m=example$m)

## End(Not run)

Check parameter length and order for a fitHMM (or MIfitHMM) model

Description

Prints parameters with labels based on DM, formula, and/or formulaDelta. See fitHMM for further argument details.

Usage

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,
  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,
  fixPar = NULL,
  prior = NULL,
  ...
)

Arguments

See Also

fitHMM, MIfitHMM

Examples

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)

Convert standard direction angles (in radians relative to the x-axis) to turning angle covariates suitable for circular-circular regression on the angle mean

Description

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.

Usage

circAngles(refAngle, data, coordNames = c("x", "y"))

Arguments

Value

A vector of turning angles between the movement direction at time step t-1 and refAngle at time t

Examples

# 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)

Merge crwData or crwHierData object with additional data streams and/or covariates

Description

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.

Usage

crawlMerge(crwData, data, Time.name)

Arguments

Details

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")

Value

A crwData object

Examples

## 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)

Fit and predict tracks for using crawl

Description

Wrapper function for fitting crawl::crwMLE models and predicting locations with crawl::crwPredict for multiple individuals.

Usage

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,
  ...
)

Arguments

Details

• 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.

Value

A crwData or crwHierData object, i.e. a list of:

The crwData object is used in MIfitHMM analyses that account for temporal irregularity or location measurement error.

See Also

MIfitHMM, simData

Examples

## 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)

Constructor of crwData objects

Description

Constructor of crwData objects

Usage

crwData(m)

Arguments

Value

An object crwData.

See Also

crawlWrap, MIfitHMM

Constructor of crwHierData objects

Description

Constructor of crwHierData objects

Usage

crwHierData(m)

Arguments

Value

An object crwHierData.

See Also

crawlWrap, MIfitHMM

Constructor of crwHierSim objects

Description

Constructor of crwHierSim objects

Usage

crwHierSim(m)

Arguments

Value

An object crwHierSim.

Constructor of crwSim objects

Description

Constructor of crwSim objects

Usage

crwSim(m)

Arguments

Value

An object crwSim.

Bernoulli density function

Description

Probability density function of the Bernoulli distribution (written in C++)

Usage

dbern_rcpp(x, prob, foo)

Arguments

Value

Vector of densities

Probability density function of the beta distribution (written in C++)

Description

Probability density function of the beta distribution (written in C++)

Usage

dbeta_rcpp(x, shape1, shape2)

Arguments

Value

Vector of densities

Categorical density function

Description

Probability density function of the categorical distribution (written in C++)

Usage

dcat_rcpp(x, prob, foo)

Arguments

Value

Vector of densities

Exponential density function

Description

Probability density function of the exponential distribution (written in C++)

Usage

dexp_rcpp(x, rate, foo)

Arguments

Value

Vector of densities

Gamma density function

Description

Probability density function of the gamma distribution (written in C++)

Usage

dgamma_rcpp(x, mu, sigma)

Arguments

Value

Vector of densities

Calculate distance between points y and z and turning angle between points x, y, and z

Description

Calculate distance between points y and z and turning angle between points x, y, and z

Usage

distAngle(x, y, z, type = "UTM", angleCov = TRUE)

Arguments

Details

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.

Value

2-vector with first element the distance between y and z and second element the turning angle between (x,y) and (y,z).

Log-normal density function

Description

Probability density function of the log-normal distribution (written in C++)

Usage

dlnorm_rcpp(x, meanlog, sdlog)

Arguments

Value

Vector of densities

logistic density function

Description

Probability density function of the logistic distribution (written in C++)

Usage

dlogis_rcpp(x, location, scale)

Arguments

Value

Vector of densities

C++ implementation of multivariate Normal probability density function for multiple inputs

Description

C++ implementation of multivariate Normal probability density function for multiple inputs

Usage

dmvnorm_rcpp(x, mean, varcovM)

Arguments

Value

matrix of densities of dimension K x n.

negative binomial density function

Description

Probability density function of the negative binomial distribution (written in C++)

Usage

dnbinom_rcpp(x, mu, size)

Arguments

Value

Vector of densities

Normal density function

Description

Probability density function of the normal distribution (written in C++)

Usage

dnorm_rcpp(x, mean, sd)

Arguments

Value

Vector of densities

Poisson density function

Description

Probability density function of the Poisson distribution (written in C++)

Usage

dpois_rcpp(x, rate, foo)

Arguments

Value

Vector of densities

student t density function

Description

Probability density function of non-central student t (written in C++)

Usage

dt_rcpp(x, df, ncp)

Arguments

Value

Vector of densities

Von Mises density function

Description

Probability density function of the Von Mises distribution, defined as a function of the modified Bessel function of order 0 (written in C++)

Usage

dvm_rcpp(x, mu, kappa)

Arguments

Value

Vector of densities

Weibull density function

Description

Probability density function of the Weibull distribution (written in C++)

Usage

dweibull_rcpp(x, shape, scale)

Arguments

Value

Vector of densities

Wrapped Cauchy density function

Description

Probability density function of the wrapped Cauchy distribution (written in C++)

Usage

dwrpcauchy_rcpp(x, mu, rho)

Arguments

Value

Vector of densities

Example dataset

Description

These data are used in the examples and tests of functions to keep them as short as possible.

Usage

example

miExample

forest

Details

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)

Description

Expand vector of free working parameters to vector of all working parameters including any fixed parameters (used in fitHMM.R and nLogLike.R)

Usage

expandPar(
  optPar,
  optInd,
  fixPar,
  wparIndex,
  betaCons,
  deltaCons,
  nbStates,
  nbCovsDelta,
  stationary,
  nbCovs,
  nbRecovs = 0,
  mixtures = 1,
  nbCovsPi = 0
)

Arguments

Value

A vector of all working parameters including any fixed parameters

Examples

## 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)

Fit a multivariate HMM to the data

Description

Fit a (multivariate) hidden Markov model to the data provided, using numerical optimization of the log-likelihood function.

Usage

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,
  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,
  optMethod = "nlm",
  control = list(),
  prior = NULL,
  modelName = NULL,
  ...
)

Arguments

Details

• 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.

Value

A momentuHMM or momentuHierHMM object, i.e. a list of:

References

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.

See Also

getParDM, prepData, simData

simHierData

Examples

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

Description

Convert hierarchical HMM structure to a conventional HMM

Usage

formatHierHMM(
  data,
  hierStates,
  hierDist,
  hierBeta = NULL,
  hierDelta = NULL,
  hierFormula = NULL,
  hierFormulaDelta = NULL,
  mixtures = 1,
  workBounds = NULL,
  betaCons = NULL,
  deltaCons = NULL,
  fixPar = NULL,
  checkData = TRUE
)

Arguments

Value

A list of arguments needed for specifying a hierarchical HMM as a conventional HMM in fitHMM or MIfitHMM, including:

Get names of any covariates used in probability distribution parameters

Description

Get names of any covariates used in probability distribution parameters

Usage

getCovNames(m, p, distname)

Arguments

Value

A list of:

Get design matrix

Description

Loop for creating full design matrix (X) from pseudo-design matrix (DM). Written in C++. Used in getDM.

Usage

getDM_rcpp(X, covs, DM, nr, nc, cov, nbObs)

Arguments

Value

full design matrix (X)

Get starting values from momentuHMM, miHMM, or miSum object returned by fitHMM, MIfitHMM, or MIpool

Description

Get starting values from momentuHMM, miHMM, or miSum object returned by fitHMM, MIfitHMM, or MIpool

Usage

getPar(m)

Arguments

Value

A list of parameter values (Par, beta, delta) that can be used as starting values in fitHMM or MIfitHMM

See Also

getPar0, getParDM

Examples

# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
m <- example$m
Par <- getPar(m)

Get starting values for new model from existing momentuHMM or momentuHierHMM model fit

Description

For 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).

Usage

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,
  ...
)

Arguments

Details

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.

Value

A named list containing starting values suitable for Par0 and beta0 arguments in fitHMM or MIfitHMM:

See Also

getPar, getParDM, fitHMM, MIfitHMM

Examples

# 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)

Get starting values on working scale based on design matrix and other parameter constraints

Description

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.

Usage

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,
  ...
)

## 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,
  ...
)

Arguments

Details

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.

Value

A list of parameter values that can be used as starting values (Par0) in fitHMM or MIfitHMM

See Also

getPar, getPar0, fitHMM, MIfitHMM

Examples

# 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)

Transition probability matrix

Description

Computation of the transition probability matrix for each time step as a function of the covariates and the regression parameters.

Usage

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,
  ...
)

## S3 method for class 'hierarchical'
getTrProbs(
  data,
  hierStates,
  hierBeta,
  workBounds = NULL,
  hierFormula = NULL,
  mixtures = 1,
  hierDist,
  getCI = FALSE,
  covIndex = NULL,
  alpha = 0.95,
  ...
)

Arguments

Value

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.

Examples

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)

Is HMMfits

Description

Check that an object is of class HMMfits.

Usage

is.HMMfits(x)

Arguments

Value

TRUE if x is of class HMMfits, FALSE otherwise.

Is crwData

Description

Check that an object is of class crwData. Used in MIfitHMM.

Usage

is.crwData(x)

Arguments

Value

TRUE if x is of class crwData, FALSE otherwise.

Is crwHierData

Description

Check that an object is of class crwHierData. Used in MIfitHMM.

Usage

is.crwHierData(x)

Arguments

Value

TRUE if x is of class crwHierData, FALSE otherwise.

Is crwHierSim

Description

Check that an object is of class crwHierSim.

Usage

is.crwHierSim(x)

Arguments

Value

TRUE if x is of class crwHierSim, FALSE otherwise.

Is crwSim

Description

Check that an object is of class crwSim.

Usage

is.crwSim(x)

Arguments

Value

TRUE if x is of class crwSim, FALSE otherwise.

Is miHMM

Description

Check that an object is of class miHMM.

Usage

is.miHMM(x)

Arguments

Value

TRUE if x is of class miHMM, FALSE otherwise.

Is miSum

Description

Check that an object is of class miSum.

Usage

is.miSum(x)

Arguments

Value

TRUE if x is of class miSum, FALSE otherwise.

Is momentuHMM

Description

Check that an object is of class momentuHMM. Used in CIreal, CIbeta, plotPR, plotStates, pseudoRes, stateProbs, and viterbi.

Usage

is.momentuHMM(x)

Arguments

Value

TRUE if x is of class momentuHMM, FALSE otherwise.

Is momentuHMMData

Description

Check that an object is of class momentuHMMData. Used in fitHMM.

Usage

is.momentuHMMData(x)

Arguments

Value

TRUE if x is of class momentuHMMData, FALSE otherwise.

Is momentuHierHMM

Description

Check that an object is of class momentuHierHMM. Used in CIreal, CIbeta, plotPR, plotStates, pseudoRes, stateProbs, and viterbi.

Usage

is.momentuHierHMM(x)

Arguments

Value

TRUE if x is of class momentuHierHMM, FALSE otherwise.

Is momentuHierHMMData

Description

Check that an object is of class momentuHierHMMData. Used in fitHMM.

Usage

is.momentuHierHMMData(x)

Arguments

Value

TRUE if x is of class momentuHierHMMData, FALSE otherwise.

Forward log-probabilities

Description

Used in stateProbs and pseudoRes.

Usage

logAlpha(m)

Arguments

Value

A list of length model$conditions$mixtures where each element is a matrix of forward log-probabilities for each mixture.

Examples

## 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)

Backward log-probabilities

Description

Used in stateProbs.

Usage

logBeta(m)

Arguments

Value

A list of length model$conditions$mixtures where each element is a matrix of backward log-probabilities for each mixture.

Examples

## 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)

Constructor of miHMM objects

Description

Constructor of miHMM objects

Usage

miHMM(m)

Arguments

Value

An object miHMM.

Constructor of miSum objects

Description

Constructor of miSum objects

Usage

miSum(m)

Arguments

Value

An object miSum.

Mixture probabilities

Description

For a fitted model, this function computes the probability of each individual being in a particular mixture

Usage

mixtureProbs(m, getCI = FALSE, alpha = 0.95)

Arguments

Details

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.

Value

The matrix of individual mixture probabilities, with element [i,j] the probability of individual i being in mixture j

References

Maruotti, A., and T. Ryden. 2009. A semiparametric approach to hidden Markov models under longitudinal observations. Statistics and Computing 19: 381-393.

Examples

## 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)

Constructor of momentuHMM objects

Description

Constructor of momentuHMM objects

Usage

momentuHMM(m)

Arguments

Value

An object momentuHMM.

Constructor of momentuHMMData objects

Description

Constructor of momentuHMMData objects

Usage

momentuHMMData(data)

Arguments

Value

An object momentuHMMData.

Constructor of momentuHierHMM objects

Description

Constructor of momentuHierHMM objects

Usage

momentuHierHMM(m)

Arguments

Value

An object momentuHierHMM.

Constructor of momentuHierHMMData objects

Description

Constructor of momentuHierHMMData objects

Usage

momentuHierHMMData(data)

Arguments

Value

An object momentuHierHMMData.

Scaling function: natural to working parameters.

Description

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.

Usage

n2w(par, bounds, beta, delta = NULL, nbStates, estAngleMean, DM, Bndind, dist)

Arguments

Value

A vector of unconstrained parameters.

Examples

## 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

Description

Negative log-likelihood function

Usage

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
)

Arguments

Value

The negative log-likelihood of the parameters given the data.

Examples

## 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)

Negative log-likelihood

Description

Computation of the negative log-likelihood (forward algorithm - written in C++)

Usage

nLogLike_rcpp(
  nbStates,
  covs,
  data,
  dataNames,
  dist,
  Par,
  aInd,
  zeroInflation,
  oneInflation,
  stationary,
  knownStates,
  betaRef,
  mixtures
)

Arguments

Value

Negative log-likelihood

Parameters definition

Description

Parameters definition

Usage

parDef(
  dist,
  nbStates,
  estAngleMean,
  zeroInflation,
  oneInflation,
  DM,
  userBounds = NULL
)

Arguments

Value

A list of:

Examples

## 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)

Plot crwData

Description

Plot observed locations, error ellipses (if applicable), predicted locations, and prediction intervals from crwData or crwHierData object.

Usage

## S3 method for class 'crwData'
plot(
  x,
  animals = NULL,
  compact = FALSE,
  ask = TRUE,
  plotEllipse = TRUE,
  crawlPlot = FALSE,
  ...
)

Arguments

Details

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.

See Also

crwPredictPlot

Examples

## 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)

Plot miHMM

Description

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.

Usage

## 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,
  ...
)

Arguments

Details

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.

Examples

## 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)

Plot miSum

Description

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.

Usage

## 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,
  ...
)

Arguments

Details

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.

Examples

## 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)

Plot momentuHMM

Description

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.

Usage

## 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,
  ...
)

Arguments

Details

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).

Examples

# 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)

Plot momentuHMMData or momentuHierHMMData

Description

Plot momentuHMMData or momentuHierHMMData

Usage

## S3 method for class 'momentuHMMData'
plot(
  x,
  dataNames = c("step", "angle"),
  animals = NULL,
  compact = FALSE,
  ask = TRUE,
  breaks = "Sturges",
  ...
)

Arguments

Examples

# 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)

Plot pseudo-residuals

Description

Plots time series, qq-plots (against the standard normal distribution) using qqPlot, and sample ACF functions of the pseudo-residuals for each data stream

Usage

plotPR(m, lag.max = NULL, ncores = 1)

Arguments

Details

• 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).

References

Zucchini, W. and MacDonald, I.L. 2009. Hidden Markov Models for Time Series: An Introduction Using R. Chapman & Hall (London).

Examples

# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
m <- example$m

plotPR(m)

Plot observations on satellite image

Description

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.

Usage

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
)

Arguments

Details

If the plot displays the message "Sorry, we have no imagery here", try a lower level of zoom.

References

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 observations on raster image

Description

Plot tracking data over a raster layer.

Usage

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
)

Arguments

Examples

## 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 states

Description

Plot the states and states probabilities.

Usage

plotStates(m, animals = NULL, ask = TRUE)

Arguments

Examples

# 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

Description

Plot stationary state probabilities

Usage

plotStationary(
  model,
  covs = NULL,
  col = NULL,
  plotCI = FALSE,
  alpha = 0.95,
  return = FALSE,
  ...
)

Arguments

Examples

# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
m <- example$m

plotStationary(m)

Preprocessing of the data streams and covariates

Description

Preprocessing of the data streams, including calculation of step length, turning angle, and covariates from location data to be suitable for analysis using fitHMM.

Usage

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,
  ...
)

## 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,
  hierLevels,
  coordLevel,
  ...
)

Arguments

Details

• 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 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.

Value

An object momentuHMMData or momentuHierHMMData, i.e., a dataframe of:

See Also

crawlMerge, crawlWrap, crwData

crwHierData

Examples

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")

Print miHMM

Description

Print miHMM

Usage

## S3 method for class 'miHMM'
print(x, ...)

Arguments

Examples

## 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)

Print miSum

Description

Print miSum

Usage

## S3 method for class 'miSum'
print(x, ...)

Arguments

Examples

## 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)

Print momentuHMM

Description

Print momentuHMM

Usage

## S3 method for class 'momentuHMM'
print(x, ...)

## S3 method for class 'momentuHierHMM'
print(x, ...)

Arguments

Examples

# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
m <- example$m

print(m)

Pseudo-residuals

Description

The pseudo-residuals of momentuHMM models, as described in Zucchini and McDonad (2009).

Usage

pseudoRes(m, ncores = 1)

Arguments

Details

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).

Value

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.

References

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).

Examples

# 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)

Random effects estimation

Description

Approximate individual-level random effects estimation for state transition probabilities based on Burnham & White (2002)

Usage

randomEffects(
  m,
  Xformula = ~1,
  alpha = 0.95,
  ncores = 1,
  nlmPar = list(),
  fit = TRUE,
  retryFits = 0,
  retrySD = NULL,
  optMethod = "nlm",
  control = list(),
  modelName = NULL,
  ...
)

Arguments

Value

A randomEffects model similar to a momentuHMM object, but including the additional random effect components:

References

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.

Examples

## 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)

Set modelName for a momentuHMM, miHMM, HMMfits, or miSum object

Description

Set modelName for a momentuHMM, miHMM, HMMfits, or miSum object

Usage

setModelName(model, modelName)

Arguments

Value

model object with new modelName field

Examples

m <- example$m
mName <- setModelName(m, modelName="example")

Set stateNames for a momentuHMM, miHMM, HMMfits, or miSum object

Description

Set stateNames for a momentuHMM, miHMM, HMMfits, or miSum object

Usage

setStateNames(model, stateNames)

Arguments

Value

model object with new stateNames field

Examples

m <- example$m
mName <- setStateNames(m, stateNames=c("encamped","exploratory"))

Simulation tool

Description

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.

Usage

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,
  states = FALSE,
  retrySims = 0,
  lambda = NULL,
  errorEllipse = NULL,
  ncores = 1
)

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,
  states = FALSE,
  retrySims = 0,
  lambda = NULL,
  errorEllipse = NULL,
  ncores = 1
)

Arguments

Details

• 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 hierarhcial 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 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.

Value

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:

If simulated location data are temporally irregular (i.e., lambda>0) and/or include measurement error (i.e., errorEllipse!=NULL), a dataframe of:

References

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.

See Also

prepData, simObsData

Examples

# 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)                

Observation error simulation tool

Description

Simulates observed location data subject to temporal irregularity and/or location measurement error

Usage

simObsData(data, lambda, errorEllipse, ...)

## S3 method for class 'momentuHMMData'
simObsData(data, lambda, errorEllipse, ...)

## S3 method for class 'momentuHierHMMData'
simObsData(data, lambda, errorEllipse, coordLevel, ...)

Arguments

Details

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.

Value

A dataframe of:

References

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.

See Also

crawlWrap, prepData, simData

simHierData

Examples

# 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)

State probabilities

Description

For a given model, computes the probability of the process being in the different states at each time point.

Usage

stateProbs(m, hierarchical = FALSE)

Arguments

Value

The matrix of state probabilities, with element [i,j] the probability of being in state j in observation i.

References

Zucchini, W. and MacDonald, I.L. 2009. Hidden Markov Models for Time Series: An Introduction Using R. Chapman & Hall (London).

Examples

# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
m <- example$m

sp <- stateProbs(m)

Stationary state probabilities

Description

Calculates the stationary probabilities of each state based on covariate values.

Usage

stationary(model, covs, covIndex)

Arguments

Value

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.

Examples

# 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))

Summary momentuHMMData

Description

Summary momentuHMMData

Usage

## 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, ...)

Arguments

Examples

# 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”)

Description

Calculate proportion of time steps assigned to each state (i.e. “activity budgets”)

Usage

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)

Arguments

Value

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.

Examples

# m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
m <- example$m
timeInStates(m)
timeInStates(m, by = "ID")

Transition probability matrix

Description

Computation of the transition probability matrix, as a function of the covariates and the regression parameters. Written in C++. Used in viterbi.

Usage

trMatrix_rcpp(nbStates, beta, covs, betaRef)

Arguments

Value

Three dimensional array trMat, such that trMat[,,t] is the transition matrix at time t.

Turning angle

Description

Used in prepData and simData.

Usage

turnAngle(x, y, z, type = "UTM", angleCov = FALSE)

Arguments

Value

The angle between vectors (x,y) and (y,z).

If type='LL' then turning angle is calculated based on initial bearings using bearing.

Examples