Hello,

I have successfully fit a 3 state HMM to relocation data from 35 individual deer with two covariates. Recently I received data from 6 additional individuals and have been unable to complete my re-analyses. The new data set is structurally the same as the old, however when I run the fitHMM function this error is returned:

#Error in nlm(nLogLike, optPar, nbStates, newformula, p$bounds, p$parSize,:non-finite value supplied by 'nlm'

In reading through the old moveHMM vingette it seems that this may be the result of covariate values that are too large. I standardized my covariate values and re-ran the function with no improvement. It seems unlikely to me that covariate values would be the problem in this situation because the same range of covariate values worked successfully on my original dataset.

I have also tried playing around with my starting parameter values, but I've only been able to produce a new error:

#Error: Last global step failed to locate a point lower than x. Either x is an approximate local minimum of the function,

the function is too non-linear for this algorithm, or steptol is too large.

Below is a summary of the code I'm currently using:

deerPrep <- prepData(tr, type = "UTM", coordNames = c("x","y"), covNames = c("dist.meters","tsh.hours"))

summary(deerPrep)

HMM data for 41 individuals:

100 -- 12736 observations
13 -- 12736 observations
17 -- 12736 observations
170 -- 11236 observations
20 -- 12736 observations
253 -- 12736 observations
255 -- 8932 observations
256 -- 10376 observations
259 -- 1590 observations
261 -- 12736 observations
262 -- 12736 observations
263 -- 12736 observations
265 -- 12736 observations
269 -- 12736 observations
27 -- 12736 observations
272 -- 5614 observations
273 -- 6724 observations
277 -- 12736 observations
278 -- 1621 observations
281 -- 12736 observations
284 -- 12736 observations
289 -- 12736 observations
293 -- 12736 observations
294 -- 12736 observations
295 -- 12736 observations
297 -- 12736 observations
299 -- 12737 observations
321 -- 12736 observations
329 -- 12736 observations
331 -- 11620 observations
337 -- 12736 observations
339 -- 12736 observations
341 -- 12736 observations
342 -- 12737 observations
344 -- 12736 observations
347 -- 12736 observations
3470 -- 12736 observations
348 -- 6724 observations
350 -- 12736 observations
97 -- 12736 observations
99 -- 12736 observations

Data summaries:

  step           angle       

Min. : 0 Min. :-3.142
1st Qu.: 7 1st Qu.:-1.601
Median : 17 Median : 0.000
Mean : 48 Mean : 0.001
3rd Qu.: 54 3rd Qu.: 1.608
Max. :2838 Max. : 3.142
NA's :10598 NA's :14200

Initial parameters

mu0 <- c(7, 50, 150) # Mean step length
sigma0 <- c(10, 32, 100) # SD of the step length
zm0 <- c(0.001, 0, 0) # Basically split the difference in the zero-mass parameter for the proportion of step lengths that are expected to be zero in each state
am0 <- c(pi, 0, 0) # Turning angle mean
kappa0 <- c(0.1, 0.1, 0.1) # Turning angle concentration parameter (kappa between 0 and 1 for wrpcauchy)

#deerHMM3s <- fitHMM(deerPrep,
nbStates = 3,
dist = list(step = "gamma",angle = "wrpcauchy"),
Par0 = list(step = c(mu0, sigma0,zm0), angle = c(am0, kappa0)),
formula = ~tsh.hours+dist.meters,
estAngleMean = list(angle=TRUE),
stationary = FALSE,
stateNames=stateNames,
nlmPar = list(print.level=2))

I'm looking for a way to troubleshoot these errors, specifically the first one mentioned:

#Error in nlm(nLogLike, optPar, nbStates, newformula, p$bounds, p$parSize, : non-finite value supplied by 'nlm'

Thank you!

Hi Brett,
I think when I tried to reinstall momentuHMM after updating R, RStudio, and Rtools, when I tried to install either from CRAN or GitHub, momentuHMM was trying to use the old version of Rtools, which I deleted from my computer. I imagine I can get around it by re-installing the old Rtools and using that, but I wanted to put this issue on your radar in case it wasn't yet. Thanks!
-Robbie Emmet

Hello, I am working with data from GPS tags deployed on harbour seals. I am trying to fit a model with environmental covariates and mixtures to capture inter-individual heterogeneity in the responses. Following the example with pilot whales in the vignettes.

As I am focussing my analysis on foraging trips at sea, I have removed sections of the data during which seals are hauled-out or doing short coastal trips. Therefore, I had to use my trip_ID as my ID column when preparing the data (prepData()), to avoid the error of not continuous data for each individual. However, now that I am running the models with mixtures, I would like to account for the seal individual (seal_ID) variability and not of each individual trip (trip_ID). So, I have tried to switch the columns and set as my ID the seal_ID. The model is actually running and estimating stationary probability distributions, but when it comes to estimate the most probable state sequence using the viterbi() function this error is generated:

> viterbi(ms)
Error in tm[[mix]][, , i] : subscript out of bounds

Is there any way to obtain the state sequence even with the gaps in the time series? Do you have any suggestions on how to overcome this error?

On another matter, I am trying to use as covariates temperature data that were collected by the tags. As there are various gaps in the data, so far, I have discarded trips with missing temperature readings, but that is decreasing a lot my sample size. I was wondering if there is any way to fit the model with missing values in the covariate columns.
I hope my questions were clear enough, please do not hesitate to ask for clarifications.

Hi Brett,

I'm trying to prescribe some expected relationships between behavioral states as in lines 31 - 60 of the fur seal example

I'm having a hard time understanding the structure of the design matrices and how to change them to indicate the relationships I wish to. I'm thinking that the suffix numbers on the means and standard deviations refer to the stateNames list, but I'm confused as to the meaning of "(intercept)" and the semantic implication of 1's and 0's within the matrix. I'm presuming this is also a strict binary (i.e. no fractional values), is this correct?

Could you help me understand how to read / write matrices such as angleDM and stepDM in the example?

Thank you much,

Travis

I am attempting to fit a hierarchical model with 4 coarse states (daily movement average) and 4 fine states (individuals steps). I am using step length and turning angle for the fine state data streams and for the coarse I am using a daily average of Brownian motion variance (BMV_mean) estimated from a dynamic Brownian Bridge movement model fit to each movement track. The goal of the model is to identify when changes in seasonal movement patterns occur by looking at both the individuals steps and the average of the steps within a day.

The issue I am having appears to occur when I include workBounds, userBounds, and/or fixPar information to the model. Prior to their inclusion, the model runs, albeit not well. I include the workBounds to restrict the direction of the coefficients for BMV_mean, step length and turning angle concentration. The userBounds was included to prevent the parameters from running tooi high, and to specifically restrict step length mean and zero-mass parameters for identifying fine scale states.

I include the fixPar objects for the beta and delta to restrict the transition between coarse scale states and to restrict the starting states to being within the first coarse state. We have 4 coarse states; Wintering, Dispersal, PreNesting, Nesting. I wish to start all individuals in the Wintering state, and restrict transition paths to 2 options: W -> D -> P -> N or W -> P -> N.

The issue is, when including any combination of these arguments in the fitHMM function, the model runs for a short period of time (<1 hour), returns no errors, but the regression coefficients for the transition matrix and initial distribution are the same as those returned by getTRProbs(). When I attempt to plot the model, an error is returned

plot(turk.hhmm.4.4,ask=FALSE)
Decoding state sequence... Error in stSeq[nbObs] <- which.max(tmxi[nbObs, ]) : 
  replacement has length zero

I am unsure whether the issue lies in how I am introducing the values within the arguments of fitHMM or something more complicated. I am a little lost on how to diagnose the issue beyond trial and error of excluding arguments, which is how I have (hopefull) narrowed the issue down to these arguments. I was unable to find many examples of how to implement these arguments for hierarchical models so I attempted to follow the instructions of the vignette and package documentation for fitHMM. I have included the full code for my model below.

#############################################
### Define the hierarchical HMM structure ###
#############################################
### define hierarchical HMM: 
#### states 1-4 = coarse state 1 (Wintering)
#### states 5-8 = coarse state 2 (Dispersal)
#### states 9-12 = coarse state 3 (PreNesting)
#### states 13-16 = coarse state 4 (Nesting)
hierStates <- data.tree::Node$new("turkey HHMM states")
hierStates$AddChild("Winter")   # resident/foraging
hierStates$Winter$AddChild("W1", state=1)
hierStates$Winter$AddChild("W2", state=2)
hierStates$Winter$AddChild("W3", state=3)
hierStates$Winter$AddChild("W4", state=4)
hierStates$AddChild("Dispersal")  # mobile/foraging
hierStates$Dispersal$AddChild("D1", state=5)
hierStates$Dispersal$AddChild("D2", state=6)
hierStates$Dispersal$AddChild("D3", state=7)
hierStates$Dispersal$AddChild("D4", state=8)
hierStates$AddChild("PreNesting")    # travelling/migrating
hierStates$PreNesting$AddChild("P1", state=9)
hierStates$PreNesting$AddChild("P2", state=10)
hierStates$PreNesting$AddChild("P3", state=11)
hierStates$PreNesting$AddChild("P4", state=12)
hierStates$AddChild("Nesting")    # travelling/migrating
hierStates$Nesting$AddChild("N1", state=13)
hierStates$Nesting$AddChild("N2", state=14)
hierStates$Nesting$AddChild("N3", state=15)
hierStates$Nesting$AddChild("N4", state=16)

print(hierStates,"state")

nbStates <- length(hierStates$Get("state",filterFun=data.tree::isLeaf))

###################################################
### Parameter Distributions and Starting Values ###
###################################################
# data stream distributions: 
## level 1 = coarse level (BMV_mean="gamma")
## level 2 = fine level (step="gamma", angle = "wrpcauchy")
hierDist <- data.tree::Node$new("turkey HHMM dist")
hierDist$AddChild("level1")
hierDist$level1$AddChild("BMV_mean", dist="gamma")
hierDist$AddChild("level2")
hierDist$level2$AddChild("step", dist="gamma")
hierDist$level2$AddChild("angle", dist="wrpcauchy")

print(hierDist,"dist")

### defining start values based on those reported by Adam et al
BMV.mu0 <- c(75, 500, 300,1)
BMV.sigma0 <- c(50, 200, 100,50)
step.mu0 <- step.sigma0 <- step.zero <- angle.con <- list()
step.mu0[[1]] <- step.mu0[[2]] <- step.mu0[[3]] <- step.mu0[[4]] <-c(5, 20, 150, 300)
# step.mu0[[2]] <- c(5, 20, 150) #For if you want to specify different start values
# step.mu0[[3]] <- c(5, 20, 150)
# step.mu0[[4]] <- c(5, 20, 150)
step.sigma0[[1]] <- step.sigma0[[2]] <- step.sigma0[[3]] <- step.sigma0[[4]] <- c(20, 40, 100, 200)
# step.sigma0[[2]] <- c(0.043, 0.047, 0.342)
# step.sigma0[[3]] <- c(0.109, 0.462, 1.878)
# step.sigma0[[4]] <- c(0.109, 0.462, 1.878)
step.zero[[1]] <- step.zero[[2]] <- step.zero[[3]] <- step.zero[[4]] <- c(0.985, 0.006, 0.005, 0.004)
# step.zero[[2]] <- c(0.043, 0.047, 0.342)
# step.zero[[3]] <- c(0.109, 0.462, 1.878)
# step.zero[[4]] <- c(0.109, 0.462, 1.878)
angle.con[[1]] <- angle.con[[2]] <- angle.con[[3]] <- angle.con[[4]] <- c(0.1, 0.2, 0.3, 0.4)
# angle.con[[2]] <- c(1.791e-08, 0.035, 0.002)
# angle.con[[3]] <- c(0.012, 2.933e-04, 3.462e-09)
# angle.con[[4]] <- c(0.012, 2.933e-04, 3.462e-09)

Par0 <- list(BMV_mean=c(rep(BMV.mu0,each=4),rep(BMV.sigma0,each=4)),
             step=c(unlist(step.mu0),unlist(step.sigma0),unlist(step.zero)),
             angle=c(unlist(angle.con)))

####################################
### Design Matrix Specifications ###
####################################
BMV.DM <- matrix(c(rep(c(1,1,0,0,0,0,0,0),4),
                   rep(c(1,1,1,1,0,0,0,0),4),
                   rep(c(1,1,1,0,0,0,0,0),4),
                   rep(c(1,0,0,0,0,0,0,0),4),
                   rep(c(0,0,0,0,1,0,0,0),4),
                   rep(c(0,0,0,0,0,1,0,0),4),
                   rep(c(0,0,0,0,0,0,1,0),4),
                   rep(c(0,0,0,0,0,0,0,1),4)),
                 nrow = 2*nbStates, byrow = T,
                 ncol = 8,
                 dimnames = list(paste0(rep(c("mean_","sd_"),each=nbStates),1:nbStates),
                                 c(paste0(rep(c("mean_","sd_"),each=4),1:length(hierStates$children)))))

step.DM <- matrix(c(rep(c(1,0,0,0,0,0,0,0,0,0,0,0,
                          1,1,0,0,0,0,0,0,0,0,0,0,
                          1,1,1,0,0,0,0,0,0,0,0,0,
                          1,1,1,1,0,0,0,0,0,0,0,0),4),
                    rep(c(0,0,0,0,1,0,0,0,0,0,0,0,
                          0,0,0,0,0,1,0,0,0,0,0,0,
                          0,0,0,0,0,0,1,0,0,0,0,0,
                          0,0,0,0,0,0,0,1,0,0,0,0),4),
                    rep(c(0,0,0,0,0,0,0,0,1,0,0,0,
                          0,0,0,0,0,0,0,0,0,1,0,0,
                          0,0,0,0,0,0,0,0,0,0,1,0,
                          0,0,0,0,0,0,0,0,0,0,0,1),4)),
                 nrow = 3*nbStates, byrow = T,
                 ncol = 12,
                 dimnames = list(paste0(rep(c("mean_","sd_","zero_"),each=nbStates),1:nbStates),
                                 c(paste0(rep(c("mean_","sd_","zero_"),each=4),1:4))))

angle.DM <- matrix(c(rep(c(1,0,0,0,
                           1,1,0,0,
                           1,1,1,0,
                           1,1,1,1),4)),
                  nrow = 1*nbStates, byrow = T,
                  ncol = 4,
                  dimnames = list(paste0(rep(c("con_"),each=nbStates),1:nbStates),
                                  c(paste0(rep(c("con_"),each=4),1:4))))

  
DM <- list(BMV_mean = BMV.DM,
           step = step.DM,
           angle = angle.DM)


#############################################
### Specify Work Bounds for distributions ###
#############################################
#These correspond to the columns of the DM
#define the directions of the differences
BMVworkBounds <- matrix(c(-Inf,Inf, 
                           0,Inf,
                           0,Inf,
                           0,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf), 
                         nrow = ncol(BMV.DM), byrow = T,
                         dimnames = list(colnames(BMV.DM), c("lower", "upper")))

stepworkBounds <- matrix(c(-Inf,Inf,
                           0,Inf,
                           0,Inf,
                           0,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf,
                           -Inf,Inf), 
                     nrow = ncol(step.DM), byrow = T,
                     dimnames = list(colnames(step.DM), c("lower", "upper")))

angleworkBounds <- matrix(c(-Inf,Inf,
                            -Inf,0,
                            -Inf,0,
                            -Inf,0),
                      nrow = ncol(angle.DM),
                      byrow=TRUE, dimnames=list(colnames(angle.DM),
                                                c("lower","upper"))) 

workBounds<-list(BMV_mean = BMVworkBounds,
                 step=stepworkBounds,
                 angle=angleworkBounds)

#####################################################
### Specify User Defined Bounds for distributions ###
#####################################################
#These correspond to the actual parameter values for each row
BMVBounds <- matrix(c(rep(c(0,1000),4),
                      rep(c(0,1000),4),
                      rep(c(0,1000),4),
                      rep(c(0,1000),4),
                      rep(c(0,Inf),4),
                      rep(c(0,Inf),4),
                      rep(c(0,Inf),4),
                      rep(c(0,Inf),4)), 
                     nrow = 4*ncol(BMV.DM), byrow = T,
                     dimnames = list(rep(colnames(BMV.DM),each=4), c("lower", "upper")))

stepBounds <- matrix(c(rep(c(0,5,
                             0,100,
                             0,200,
                             0,400),4),
                       rep(c(0,Inf,
                             0,Inf,
                             0,Inf,
                             0,Inf),4),
                       rep(c(.97,1,
                             0,.03,
                             0,.03,
                             0,.03),4)), 
                     nrow = 4*ncol(step.DM), byrow = T,
                     dimnames = list(rep(colnames(step.DM),each=4), c("lower", "upper")))

angleBounds <- matrix(c(rep(c(0,1,
                              0,1,
                              0,1,
                              0,1),4)),
                      nrow = 4*ncol(angle.DM), byrow=TRUE, 
                      dimnames=list(rep(colnames(angle.DM),each=4), c("lower","upper"))) 

userBounds <- list(BMV_mean = BMVBounds,
                   step = stepBounds,
                   angle = angleBounds)

################################
### Initial parameter Values ###
################################
# get initial parameter values for data stream probability distributions on the working scale
Par <- getParDM(turkData,
                hierStates=hierStates,
                hierDist=hierDist,
                Par=Par0,
                DM=DM)

##################################################
### Define Hierarchical Transition Probability ###
##################################################
hierFormula <- data.tree::Node$new("turkey HHMM formula")
hierFormula$AddChild("level1", formula=~1)
hierFormula$AddChild("level2", formula=~1)


#number of repeats of the betas will depend on the number of covariates in the above formulas, repeats will be level specific
hierBeta <- data.tree::Node$new("turkey beta")
hierBeta$AddChild("level1",beta=matrix(c(-.1,-.1,-100,
                                         -100,-.1,-100,
                                         -100,-100,-.1,
                                         -100,-100,-100),
                                       1,length(hierStates$children)*(length(hierStates$children)-1), byrow = T))
hierBeta$AddChild("level2")
hierBeta$level2$AddChild("Winter",beta=matrix(rep(c(2,3,-3,
                                                    2,3,-3,
                                                    2,3,-3,
                                                    2,3,-3),
                                                  1),1,length(hierStates$Winter$children)*(length(hierStates$Winter$children)-1),byrow=TRUE))
hierBeta$level2$AddChild("Dispersal",beta=matrix(rep(c(0,2,3,
                                                       0,2,3,
                                                       0,2,3,
                                                       0,2,3),
                                                     1),1,length(hierStates$Dispersal$children)*(length(hierStates$Dispersal$children)-1),byrow=TRUE))
hierBeta$level2$AddChild("PreNesting",  beta=matrix(rep(c(1,3,1,
                                                          1,3,1,
                                                          1,3,1,
                                                          1,3,1),
                                                        1),1,length(hierStates$PreNesting$children)*(length(hierStates$PreNesting$children)-1),byrow=TRUE))
hierBeta$level2$AddChild("Nesting",  beta=matrix(rep(c(1,-1,-3,
                                                       1,-1,-3,
                                                       1,-1,-3,
                                                       1,-1,-3),
                                                     1),1,length(hierStates$Nesting$children)*(length(hierStates$Nesting$children)-1),byrow=TRUE))
print(hierBeta,"beta")

#Fix the betas for the coarse scale to prevent switching to previous states
hierFixBeta <- data.tree::Node$new("turkey beta")
hierFixBeta$AddChild("level1",fixPar=matrix(c(NA,NA,-100,
                                         -100,NA,-100,
                                         -100,-100,NA,
                                         -100,-100,-100),
                                         1, length(hierStates$children)*(length(hierStates$children)-1), byrow = T))
hierFixBeta$AddChild("level2")
hierFixBeta$level2$AddChild("Winter",fixPar=matrix(rep(c(NA,NA,NA,
                                                    NA,NA,NA,
                                                    NA,NA,NA,
                                                    NA,NA,NA),
                                                    1),1,length(hierStates$Winter$children)*(length(hierStates$Winter$children)-1),byrow=TRUE))
hierFixBeta$level2$AddChild("Dispersal",fixPar=matrix(rep(c(NA,NA,NA,
                                                       NA,NA,NA,
                                                       NA,NA,NA,
                                                       NA,NA,NA),
                                                       1),1,length(hierStates$Dispersal$children)*(length(hierStates$Dispersal$children)-1),byrow=TRUE))
hierFixBeta$level2$AddChild("PreNesting",  fixPar=matrix(rep(c(NA,NA,NA,
                                                          NA,NA,NA,
                                                          NA,NA,NA,
                                                          NA,NA,NA),
                                                          1),1,length(hierStates$PreNesting$children)*(length(hierStates$PreNesting$children)-1),byrow=TRUE))
hierFixBeta$level2$AddChild("Nesting",  fixPar=matrix(rep(c(NA,NA,NA,
                                                       NA,NA,NA,
                                                       NA,NA,NA,
                                                       NA,NA,NA),
                                                       1),1,length(hierStates$Nesting$children)*(length(hierStates$Nesting$children)-1),byrow=TRUE))
print(hierFixBeta,"fixPar")

###################################
### Define Initial Distribution ###
###################################
hierDelta <- data.tree::Node$new("turkey HHMM delta")
hierDelta$AddChild("level1",delta=matrix(c(-100, -100, -100),1))
hierDelta$AddChild("level2")
hierDelta$level2$AddChild("Winter",delta=matrix(c(0,0,0),1))
hierDelta$level2$AddChild("Dispersal",delta=matrix(c(0,0,0),1))
hierDelta$level2$AddChild("PreNesting",delta=matrix(c(0,0,0),1))
hierDelta$level2$AddChild("Nesting",delta=matrix(c(0,0,0),1))

print(hierDelta,"delta")

#Fix the initial distribution to start in Winter
hierFixDelta <- data.tree::Node$new("turkey HHMM delta")
hierFixDelta$AddChild("level1",fixPar=matrix(c(-100, -100, -100),1))
hierFixDelta$AddChild("level2")
hierFixDelta$level2$AddChild("Winter",fixPar=matrix(c(NA,NA,NA),1))
hierFixDelta$level2$AddChild("Dispersal",fixPar=matrix(c(NA,NA,NA),1))
hierFixDelta$level2$AddChild("PreNesting",fixPar=matrix(c(NA,NA,NA),1))
hierFixDelta$level2$AddChild("Nesting",fixPar=matrix(c(NA,NA,NA),1))

print(hierFixDelta,"fixPar")

##################################
### Check Model Specifications ###
##################################
# check hierarchical model specification and parameters
checkPar0(turkData,
          hierStates=hierStates,
          hierDist=hierDist,
          Par0=Par,
          hierFormula=hierFormula,
          DM=DM, 
          workBounds = workBounds,
          userBounds = userBounds,
          hierBeta=hierBeta,
          hierDelta=hierDelta,
          fixPar=list(beta=hierFixBeta, delta=hierFixDelta)
          )


#################################
### Run the Hierarchical HHMM ###
#################################
#prevents working parameters from straying along boundary
#https://github.com/bmcclintock/momentuHMM/issues/41
prior <- function(par){sum(dnorm(par,0,10,log=TRUE))}

# compute hierarchical state transition probabilities based on initial values
iTrProbs <- getTrProbs(turkData,
                       hierStates=hierStates,
                       hierDist=hierDist,
                       Par0=Par,
                       hierFormula=hierFormula,
                       DM=DM, 
                       workBounds = workBounds, 
                       userBounds = userBounds,
                       hierBeta=hierBeta,
                       hierDelta=hierDelta,
                       fixPar=list(beta=hierFixBeta, delta=hierFixDelta)
                       )
iTrProbs$level1$gamma[,,1] # tpm at first time step for level1
lapply(iTrProbs$level2,function(x) x$gamma[,,1]) # tpm at first time step for level2

print(Sys.time())
turk.hhmm.4.4 <- fitHMM(turkData,
                        hierStates=hierStates,
                        hierDist=hierDist,
                        hierFormula=hierFormula,
                        Par0=Par,
                        hierBeta=hierBeta,
                        hierDelta=hierDelta,
                        DM=DM, 
                        workBounds = workBounds,
                        userBounds = userBounds,
                        fixPar=list(beta=hierFixBeta, delta=hierFixDelta),
                        # prior = prior,
                        ncores = detectCores()/2
                        )
turk.hhmm.4.4
print(Sys.time())

Dear Brett (and others),

I am trying to expand a step length and turning angle model to also include proportion data (values between 0 and 1) - that is the proportion of time recorded wet by a saltwater immersion logger (deployed alongside a GPS). I am aiming to refine models of seabird behaviour by including wet events, i.e. contact with the ocean surface.

From looking at the vignette - https://cran.r-project.org/web/packages/momentuHMM/vignettes/momentuHMM.pdf - there doesn't appear to be a distribution that incorporates both 0 and 1 values and anything in between, unless I have missed something here? Unfortunately, the data are rather skewed towards 0 and 1 values, so using a normal distribution or something similar, is probably not possible.

I would appreciated anything that you could suggest.

Many thanks,

Tommy

Hi,

I am currently preparing data for an analysis using momentuHMM. My data is temporally-irregular and has some measurement error and I would like to use the multiple imputation approach (as I have had difficulties with fitting state-switching SSMs to this data). I've read the various other issues (e.g. #42) that relate to the use of pseduo-design matrices to avoid label switching and I think I'm starting to get it. I note that it is trivial and I'm fairly certain I've had the "a-ha!" moment you describe in #42, but I just want to be sure.

It is my understanding that for a 3-state HMM which requires zero-inflation the following would work:

nbStates <- 3
dist <- list(step = "gamma", angle = "wrpcauchy")
stepDM <- matrix(c(1,0,0,0,0,0,0,0,0,
                   0,1,0,0,0,0,0,0,0,
                   0,0,1,0,0,0,0,0,0,
                   0,0,0,1,0,0,0,0,0,
                   0,0,0,0,1,0,0,0,0,
                   0,0,0,0,0,1,0,0,0,
                   0,0,0,0,0,0,1,0,0,
                   0,0,0,0,0,0,0,1,0,
                   0,0,0,0,0,0,0,0,1),
                 nrow = 3*nbStates,
                 byrow = T,
                 dimnames = list(c(paste0("mu_", 1:nbStates),
                                   paste0("sd_", 1:nbStates),
                                   paste0("zeromass_", 1:nbStates)),
                                 c(paste0("mu_", 1:nbStates),
                                   paste0("sd_", 1:nbStates),
                                   paste0("zeromass_", 1:nbStates))))

Then, if I wish to impose a constraint such that mu_3 > mu_1 I would change the 0 to a 1 in the corresponding cell (i.e. the third row in the first column) and specify the bounds using an upper bound for mu_1 that is less than the upper bound for mu_3, e.g.:

stepuserBounds <- matrix(c(0,0,0,0,0,0,0,0,0, 
                           2,2,4,2,2,4,0.1,0.2,0.3), # numbers chosen arbitrarily for illustration
                         9,2,
                         dimnames=list(colnames(stepDM),c("lower","upper")))

All of which can then be passed to getParDM() to convert the initial parameters from the natural scale to the working scale.

Have I had the "a-ha!" moment or do I have more reading to do? Any advice you can provide would be extremely appreciated!

Thanks,
Dan

I'm able to make non-spatial HMMs by just passing NA values to my coordinate columns with prepData(). I initially tried to do this using the momentuHMMData() constructor, but I coulden't as this is not exported to the namespace for the user to access, even though its listed in the manual and the documentation implies that it can be used to make a non-spatial momentuHMMData object.

I was wondering if there is a reason that momentuHMMData has no @export flag. I'm not aware of any way to make a non-spatial momentuHMMData object, even through the resulting models seem to work fine without spatial data.

Interested to hear your thoughts!

In the previous version of momentuHMM, I used "initial.state = ..." in my crawlWrap code.

firsts <- aggregate(as.numeric(row.names(daily_locs)), list(daily_locs$ID), min)$x

get_init <- function(x){
     out <- c(daily_locs$x[x],0, daily_locs$y[x],0)
     return(out)
}
inits <- lapply(firsts, get_init)

crw_mod <- crawlWrap(obsData = daily_locs, timestep = "15 min", initial.state = inits,
                                        attempts = 15, retryFits = 20)

However, in the recent release, initial.state has been removed and fixPar and theta arguments are required - apparently only on Windows (runs fine on Mac and Linux).

This code also runs correctly on my Windows OS:
crawl::crwMLE(data=subset(daily_locs, ID==1), attempts = 15)

Fix: The code runs fine with fixPar = c(NA, NA) and theta = c(0,0)

crw_mod <- crawlWrap(obsData = daily_locs, timestep = "15 min", fixPar = c(NA, NA), theta = c(0,0), attempts = 15, retryFits = 20)

Discussed in #98

For MIfitHMM, there is a way to specify activity centers that may act as attractive/repulsive agents for the movement model. Is there a way to specify a more complex underwater structure, such as a pipe or large artificial reef that may be used in a similar way? Transforming a shapeline to a series of points significantly increases processing time for model fitting, as the model is then required to fit angle and distance covariates to each individual point, not just the one closest to the animal at that point in space and time. However, to ignore such a structure would bias the model output.

Hi,
I am trying to run a model with four covariates (formula ~ a * b * c * d) but I'm running into trouble with the computational time needed. Do you have any suggestions for fitting the HMM in parallel or other ways to speed up the process?
Thanks!

I hope it is okay to post here @TheoMichelot @bmcclintock, if there is a different forum please let me know.

I have been trying to run the HMM on my data but have come across a few things I couldn't find answers to online. Please see info and plots below.

My data: Vultures, GPS tags collecting hourly points from 5am to 7pm, daily for 8 months.

Here is a histogram of the time diff after I have thinned to one hour using the move package:

Here is the same histogram after adding in the NAs:

- Should I pad with NAs from 7pm to 5am? Or is everyday a new track?

When selecting initial parameter values, the turning angles histogram is Uniform (see plot) - how do I determine initial parameter values from this?

I struggled to find initial step parameters as well from the data given the step hist is from 0 - 30km

I've uploaded the CSV in case it helps but if anyone has any guidance or help I would be greatly appreciated!

data_na.csv

I tried using random starting points but the code wouldn't work and I couldn't get it to work.

I'm using the crawlWrap function to do a non-linear interpolation to fill in missing locations in GPS tracks of seabirds (because it's GPS data, I'm not concerned about location error). I'm not sure how to determine what values to specify for "theta", or what exactly theta is.

If I'm not specifying an error, activity, or drift model, is theta the vector of 2 starting parameters for the movement model (sigma and beta)? As explained in the crwMLE help: "The movement parameters (sigma for velocity variation and beta for velocity autocorrelation)". Do you have any suggestions for estimating theta values from the actual data?

I ran through the Elephant example in the momentuHMM vignette, and tried a variety of similar theta values with my data. It did work sometimes, but I often get an error: "Error in optim(init$par, crwN2ll, method = method, hessian = need.hess, : L-BFGS-B needs finite values of 'fn' ". Is this an error due to inappropriate values of theta?

My initial analyses were in the moveHMM environment, which was suggested to standardize covariates prior to fitting the HMM.

However since I have transitioned to working with momentuHMM, I have noticed that the momentuHMM tutorials do not standardize covariates prior to fitting the model, except for one occasion I found in the tutorial on page 59 when the boat distance covariate was standardized:
fulmar_data$boat.dist <- scale(fulmar_data$boat.dist)

Am I missing something here or is it not necessary to standardize the covariates prior to training in momentuHMM?

great code thanks
may you clarify :
will it work for multivariate time series prediction both regression and classification

each row means the time series, and columns represent different continues or/and discrete/category observation.

1
where all values are continues values
weight height age target
time 1| 56 160 34 1.2
time 2| 77 170 54 3.5
time 3| 87 167 43 0.7
time 4| 55 198 72 0.5
time 5| 88 176 32 2.3
etc

2
or even will it work for multivariate time series where values are mixture of continues and categorical values
for example 2 dimensions have continues values and 3 dimensions are categorical values

		color        weight     gender  height  age      target 

time 1| black 56 m 160 34 yes
time 2| white 77 f 170 54 no
time 3| yellow 87 m 167 43 yes
time 4| white 55 m 198 72 no
time 5| white 88 f 176 32 yes

etc

Hi Brett,

Thanks for creating the package momentuHMM and all the time you spend replying users' doubts.

I am assessing the effect of different environmental co-variables in the activity behaviour of a marine fish species using acoustic telemetry. I am running a HMM by individual using fitHMM() since I have only a few individuals and I think that the patterns can be different among them. So my data stream is the average activity of my fish at 15-min time bins (meanAcc). I am just defining two states (Active and Resting) after observing the histograms for meanAcc for all the individuals (two clear frequency peaks in meanAcc).

My doubt comes when I plot the histogram of the activity values and the estimated state-dependent probability distributions using plot() in the same way you do it in the provided examples. If I do not use fixPar(), I get estimated state-dependent probability distributions that do not seem proper to me (at least visually). I have tried different initial parameters but the state-dependent probability distributions end always in the same fit independently of the initial values. However, when I use fixPar I really get good state-dependent probability distributions.

This is the fit when I do not use fixPar (independently of the initial parameters):

This is the fit when I do use fixPar:

Thus, my question is: is it a problem to do not leave to the algorithm to estimate the parameters for the state-dependent probability distributions? I wonder if in this case, for example, it would be advisable to use fixPar or not.

Thanks for your time Brett.

I’m working through the momentuHMM caracara tutorial and am wondering how the starting values for optimization (gamma means & SD) were determined? It is my understanding that each value is associated with one of the 4 (assumed) states, but I’ve not been able to determine how to calculate or estimate these values.

I realize this may not be a simple answer, but any clarification is welcome!

## starting values for optimization
# gamma means for states 1,...,N
mu <- c(0.005,0.01,0.05,0.1)
# gamma standard deviations for states 1,...,N
sig <- c(0.01,0.02,0.03,0.06)
Par0 <- list(VDBA = c(mu,sig))

Is there a possibility in momentuHMM to use temporary centers of attraction? For example, I have a dataset of killer whale telemetry data, and fishery catch locations. The fishery lasts for a few hours (remains relatively stationary) and then stops. The boat leaves, but I have no AIS data from the vessels, only the locations of where they were fishing. I would like to investigate if this attracts the killer whales by using a 3 state HMM (foraging, transit, attracted to fishery). The fishing events are only temporary, but several may occur at the same time, in which case the whales may be attracted to one of them. Is there a way to model this in momentuHMM? I think I need to use the centroids option, but in predata, I get the error message "centroids time does not span data; each observation time must have a corresponding entry in centroids".

Dear Dr Brett,

I am studying the behaviour of a marine fish species using acoustic telemetry. I wish to use as data streams the number of detections (Detect) and mean activity (meanAcc) at intervals of 30 minutes, and then check the effect of different environmental variables on the behaviours defined by these two data streams. Specifically, I would like to identify three behaviours: "Hidden", "Resting" and "Foraging". The behaviour "Hidden" is determined just using the data stream "number of detections" so that, when it is 0, I consider that my fish is hidden. The other two behaviours ("Resting" and "Foraging") are determined by the data stream "meanAcc".

When I try to run a model with those premises, I find a problem: for the behaviour "hidden" I have NA for the data stream "meanAcc", since for these intervals I do not have activity data.

Below I show an easy example to work with using Photoperiod as covariate:

1 myData <- prepData(data=Data[,c("meanAcc","Detect","Photoperiod")], covNames=c("Photoperiod"),coordNames=NULL)
2
3 stateNames <- c("Hidden","Active","Resting")                                                         
4 dist = list(meanAcc = "gamma", Detect = "pois")          
5                                     
6 Par0_m1 <- list(meanAcc=c(NA,1.1,0.1,NA,1.1,0.1),Detect =c(0,20,20))
7 m1 <- fitHMM(data = myData, nbStates = 3, dist = dist, Par0 = Par0_m1, stateNames = stateNames)
8 states.m1 <- viterbi(m1)
9
10 formula <- ~ Photoperiod
11 Par0_m2 <- getPar0(model=m1, formula=formula)
12 m2 <- fitHMM(data = myData, nbStates = 3, dist = dist, Par0 = Par0_m2$Par, beta0=Par0_m2$beta, estAngleMean = list(angle=FALSE), stateNames = stateNames, formula=formula)
13 states.m2 <- viterbi(m2)

Line 7 produces two errors if I fix my initial parameters as I did:

1. Error in n2w(Par0[i], p$bounds, NULL, NULL, nbStates, inputs$estAngleMean, : Check the parameter bounds for Detect (the initial parameters should be strictly between the bounds of their parameter space).

2. Error in fitHMM.momentuHMMData(data = myData, nbStates = 3,: Scaling error. Check initial parameter values and bounds.

I guess the Error 1 it is due to set as 0 the number of detections for the behaviour Hidden. I assume that this can be solved just setting the number of detections as 0.001 instead of as 0, right? (I did it and Error 1 disappeared). Regarding Error 2, I do not know what I should do considering that when my animal is Hidden (Detect=0) there is no activity data (meanAcc= NA).

There appears to be a compatibility issue between older mixture models (trainied under pre v1.5.0, develop version of momentuHMM) and the recent updated version (v1.5.0) of momentuHMM develop.

When using the getPar0() function from newer version of momentuHMM library with older momentuHMM mixture models, I got the following error:

Error in !m$condition$fit : invalid argument type

And I’ve never received this error when using older version of momentuHMM develop.

I tried re-installing the older version of momentuHMM (v1.4.3) the getPar() function works, but now the fitHMM() does not work as it does not identify the "mixture" parameter:

ERROR : unused argument (mixture = nmix)

Where nmix is specified number of mixture. Note that I reinstalled older "develop" version of momentuHMM this way:

devtools::install_github("bmcclintock/[email protected]", ref = "develop")

Is there any way I can make my already trained mixture models work with the current version of momentuHMM?

On the example data from the northern fur seal script the following line executes without issue:

crwOut_nfs <- crawlWrap(nfsData, timeStep = "4 hours")

I've processed another data source to be in a format that seems identical to the above. Both objects consist of 7 columns named "ID","time","x","y","ln.sd.x","ln.sd.y","error.corr" (in that order) with the same datatypes in each column.

On execution of this line:

crwOut_fox <- crawlWrap(crawl_data, timeStep = "4 hours")

I receive the following error:

Error in crwData(list(crwFits = model_fits, crwPredict = predData)) :
Can't construct crwData object: fields are missing

which seems to refer to line 626 in crawlWrap.R

I thought the issue might be the definition of predData within the HeirInd() control block at line 584 but this wouldn't explain why the NFS code executes. What could be causing this error?

The command I used to retrieve the seal data
load(url("https://raw.github.com/bmcclintock/momentuHMM/master/vignettes/nfsData.RData"))
also returns foragedives and the "Par" list. Should either of these be relevant?

fox_crawl.xlsx

I am working with GPS data at a fine scale. The datasets have some big gaps and the location error is not negligible.
When I interpolate missing gaps using crawlWrap I get some locations at unrealistic places, for example, land for marine animals. This can have a significant impact on my analysis because the resistance of the landscape is the main covariate. I do have a "barrier" raster file and I'm able to constrain movement (e.g. from land) when using simData but not when interpolating.
Is there any way to deal with this issue using momentuHMM without introducing NAs to the dataset?

I know it is sort of a basic question but I could not find any resource to help me with this.

Thanks.
Cheers

I am trying to use the crawlWrap function to impute some missing locations in my data prior to fitting a HMM and I'm having issues with the optimization. In the elephant example of the vignette I see an argument (initial.state) where the initial values of a and P can be specified, but the current version of crawlWrap does not appear to use that argument. Is there a new way to input these initial values?

Thanks!

-JF

Hi Brett,
I'm trying to run a HSMM by expanding the number of states in my HMM. I'm hoping to go from a 4 state model to an 8 state model (3 states with geometric distributions and 5 states that expand the last state with dwell times pulled from a negative binomial distribution.

I made the expanded transition matrix following Zucchini's advice but I'm not sure how to incorporate it into the model. I believe I can fix transitions probabilities to 0 with fixPar$beta0 but I'm not sure how to fix the transition probabilities based on the dwell-time distribution.

Any help would be appreciated. Thanks!

Dear Brett,

I used the stationary function to calculate the stationary probabilities based on different covariate values in a 3-states HMM. I am now wondering if there is any way to obtain a CI and/or SE for the stationary probabilities obtained? I see that these are calculated and plotted within 'plotStationary', but it would be great to extract the values (similarly to what CIreal does for state transition probabilities). Is there a function that I somehow missed?

Thank you for any insight,
All the best,
Federico

I am using CRAWL to regularize my polar bear tracks for a location every 4 hr. I imputed 20 tracks to account for uncertainty and ran the MIfitHMM function for the 20 tracks (MI_mod). I would like to assign the most likely sequence of states (n=2) for each location using the Viterbi algorithm function (viterbi(MI_mod)) but that function only takes momentuHMM or momentHierHMM objects. How do I accomplish with a miHMM object, given by the multiple imputation method?

I have provided data for 1 bear and my code for fitting the MI_mod below:
pb_21678_MIdata_5.18.22.csv

folder where multiple imputation (from CRAWL) data is located

indir<-'C:\Users\....'

csvlist<-c(list.files(indir,pattern = ".csv$",ignore.case=T))

for(i in csvlist){
raw<-read.csv(paste(indir,i,sep='/'))

Convert times to POSIX

raw$date <- as.POSIXct(raw$datetime, format = "%Y-%m-%d %H:%M:%S")

Create list of imputed tracks

ind_imp <- c(0, which(raw$datetime[-1] < raw$datetime[-nrow(raw)]), nrow(raw))
data_list <- list()
raw$imp <- NA
for(i in 1:(length(ind_imp) - 1)) {
subdata <- prepData(raw[(ind_imp[i] + 1):ind_imp[i + 1],],
coordNames = c("X", "Y"))
data_list[[i]] <- subdata
raw$imp[(ind_imp[i] + 1):ind_imp[i + 1]] <- i
}

Fit model to multiple imputation data

MI_mod <- MIfitHMM(miData = data_list, nbStates = 2,
                 Par0 = list(step = c(100, 300, 100, 300),
                             angle = c(pi, pi, 0.5, 0.1)),
                 dist = list(step = "gamma", angle = 'wrpcauchy'),
                 estAngleMean = list(angle = TRUE))

save the model output summary

sink("MI_mod_21678_5.22.txt")
print(MI_mod)
sink()

}

MI_mod
plot(MI_mod)

Assign the probable sequence of states with the Viterbi algorithm

viterbi(MI_mod)
**Error here because MI_mod is a miHMM object but viterbi() requires a momentHMM or momentHierHMM object

Hi Brett,

I have been experimenting with momentuHMM and am very impressed with the extensive functionality you have built into this package.

At the moment, I am fitting an HMM formulation to animal tracking data including influence from a single raster covariate, with the intention to simulate from this model. The raster is spatially referenced as simple lon/lat (coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0). Passing the tracking dataframe to prepData() with x/y locations specified in lon/lat seems to work very sensibly. The step lengths calculated in are all sensible values in km, and the covariate values are all extracted properly from the raster with the cellFromXY call etc. The output from fitHMM() is all good.

However, when I then use this model to simulate tracks, even with very few (e.g. N=3) timesteps - which in the dataset were originally 6h intervals – immediately the error is thrown: "Movement is beyond the spatial extent of the … raster. Try expanding the extent of the raster." Further investigation with/without the SpatialCov shows that the simulation step lengths (which ought to be km) are transiting lon/lat grid coordinates as if they were easting/northing units (i.e. also specified in km). So, for example a single step length of 8km might travel from 30S, 150E straight to 157E. I guess a workaround could be to respecify the raster for simulation purposes using easting/northing, but it seems a peculiar thing to have to do at this point? So, I thought I would post it here. Issue was encountered using both CRAN version 1.4.3 (2018-11-21) and the development version 1.5.0 (2019-04-16). Apologies if I have missed something in the simData help file …

Thank you for any insight,
Sophie Bestley

Hi Brett,
I'm wondering why you have !is.null(model$mod$hessian) as a condition to put Sigma <- model$mod$Sigma in the plotStationary function (line 34 in plotStationary.momentuHMM (at least how I see it in RStudio). I guess it has to do with the way the covariance matrix of all parameters was estimated by numerical optimization.
Thanks in advance,
Rocío

I am having an issue with momentuHMM_1.5.5 on R version 4.2.2

fitHMM incorrectly implements fixPar for angle when using estAngleMean. In the reproducible code below, the model correctly fixes the mean step length for state1 to 0.5 in the fit_2s_fx HMM. However, the angle for state1 is not fixed to pi. Interestingly, the estimated mean turn angles are different in fit_2s vs fit_2s_fx.

library(momentuHMM)
set.seed(123)

# Generate 100 x y locations
track <- data.frame(x = cumsum(rnorm(100, mean = 0, sd = 1)),
                    y = cumsum(rnorm(100, mean = 0, sd = 1)))
prep <- prepData(track)

# define model
stateNames <- c("state1", "state2")
nbStates <- length(stateNames)
dist <- list(step = "gamma", 
             angle = "vm")
Par0 <- list(step=c(0.5,1,0.5,1),angle=c(pi,0,0.1,0.2))
fixPar <- list(step = c(0.5,NA,NA,NA),
               angle= c(pi, NA, rep(NA, 2)))
#fit HMM
fit_2s <- fitHMM(prep, stateNames = stateNames, 
                 nbStates = nbStates, dist = dist, 
                 estAngleMean = list(angle = T),
                 Par0 = Par0)
fit_2s_fx <- fitHMM(prep, stateNames = stateNames, 
                    nbStates = nbStates, dist = dist, 
                    estAngleMean = list(angle = T),
                    Par0 = Par0, fixPar = fixPar)
# view results
fit_2s$mle$step
fit_2s_fx$mle$step
fit_2s$mle$angle
fit_2s_fx$mle$angle
fit_2s_fx$conditions$fixPar$angle

Here is my output of the 2 models:

> # view results
> fit_2s$mle$step
        state1    state2
mean 1.1489887 1.4300370
sd   0.6730796 0.6194929
> fit_2s_fx$mle$step
        state1  state2
mean 0.5000000 1.17991
sd   0.1762279 0.67677
> fit_2s$mle$angle
                state1     state2
mean          2.694855 -0.9625877
concentration 0.257626  5.0528374
> fit_2s_fx$mle$angle
                 state1     state2
mean          -1.610949 -2.8646546
concentration  2.491160  0.1213415
> fit_2s_fx$conditions$fixPar$angle
[1] 3.141593       NA       NA       NA

Hi Brett,

I am trying to compute pseudo residuals for HMM assuming bivariate Gaussian as emission distribution. I noticed that you made use of the Mahalanobis distance to compute the cumulative density function for bivariate Gaussian (mvnorm2). So that the cumulative probability and pseudo residual can be calculated for each time t.
There are two kinds of pseudo-residual, one is uniform pseudo-residual and the other is normal pseudo-residual according to Zucchini et al book.

However, in your code for pseudoRes function:
if(dist[[j]] %in% mvndists){
genRes[[paste0(j,"Res")]][i] <- stats::qchisq(genRes[[paste0(j,"Res")]][i],df=ndim)
}
it seems that you assumed the pseudo-residuals are distributed Chi-squared with ndim degree of freedom instead of Normal.
May I know the reason why you choose qchisq?
Thank you for your useful package and look forward to your reply.

Bests,
Jiaxing

Hi Brett, Théo,

I am looking at the turtle example and I am confused by one of the starting values that you specify for the state-dependent probability distribution parameters.

In your code you specify the following:

nbStates<-2
dist<-list(step="gamma",angle="wrpcauchy")
estAngleMean<-list(angle=TRUE)
circularAngleMean<-list(angle=TRUE)

DM<-list(step=list(mean=~state2(w:angle_osc),sd=~1),
         angle=list(mean=~state2(d),concentration=~1))

Par0<-list(step=c(9.132130, 9.200282, 0, 8.579123, 8.640819),
           angle=c(0, -16.0417715, -0.4749668))

turtleFits<-MIfitHMM(miTurtleData$miData, 
                     ncores=ncores,       
                     nbStates=nbStates,
                     dist=dist,
                     Par0=Par0,
                     DM=DM,
                     estAngleMean=estAngleMean,
                     circularAngleMean=circularAngleMean,
                     retryFits=retryFits)

The value that I am confused by is the third value in the Par0 object (the zero; in bold above, and below). What is this zero for? Is it a zero-mass parameter?

I thought that in a 2-state HMM, using gamma distribution for step length, you would specify a total of 4 step length parameters (2 means, and 2 SDs). Here however, you have specified 5?

> Par0$step
[1] 9.132130 9.200282 0.000000 8.579123 8.640819

I am trying to repeat your analysis using my own data, and I am using three states: Resting, Foraging, and Transiting. I would like, as you have done, to model the transit state as being potentially influenced by ocean surface currents, and the other two states to be unaffected. I have hence tried to mirror your code using three states, as follows:

nbStates<-3
dist<-list(step="gamma",angle="wrpcauchy") 
estAngleMean<-list(angle=TRUE)
circularAngleMean<-list(angle=TRUE)

Par0<-list(step=c(0.89,0.64,0.30,0.13,0.15,0.19),
           angle=c(0,0,0,20,2.31,0.93))

DM<-list(step=list(mean=~state1(tideVel:angle_osc),sd=~1),
           angle=list(mean=~state1(tideDir2),concentration=~1))

Fits <- fitHMM(mHMM,   
               nbStates = nbStates,
               dist = dist,
               Par0 = Par0,
               DM = DM,
               estAngleMean = estAngleMean,
               circularAngleMean = circularAngleMean,
               retryFits = retryFits)

Note that in my example the first set of parameter estimates correspond to the transit state, hence I am using ~state1() as opposed to ~state2(). Also, my data do not require multiple imputations so I am using fitHMM() as opposed to MIfitHMM().

When I try to fit this model however, I get the following error:

Error in getDM(data, inputs$DM, inputs$dist, nbStates, p$parNames, p$bounds,  : 
  Dimension mismatch between Par and DM for: step, angle

If I have understood the error correctly, it arises because I have not provided enough parameters in the Par0 argument. Am missing a parameter like the additional 0 that you use above?

The problem is, I do not understand what the 0 is for and hence I do not know if I should be using zero also, or a different value. I also do not know where I should specify the additional parameter in my Par0 object because i) my model has 3 states whereas yours has 2 and ii) my starting values are in the sequence transit:forage:rest, whereas yours are forage:transit.

Thanks!

William Kay

Hi Brett,

I'm trying to use simData function to produce realizations from a fitted HMM model. I want the state sequence to start from a certain state (say, state 1) and I thought I could do it by fixing the parameter delta. However, it seems that the first state is not sampled from the initial state distribution (delta) as I though. Here is my toy example with 2-state HMM, where I wanted the sequence to start from state 1:

dist <- list(step="gamma",angle="vm")
Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.01,0.75))
set.seed(123)
data.sim<-simData(nbAnimals=2,obsPerAnimal=10,nbStates=2,dist=dist,Par=Par,
                   delta=c(1,0), states=TRUE)    
head(data.sim)
 ID       step      angle        x         y states
1  1  833.60385         NA   0.0000    0.0000      2
2  1 1041.50736 -0.1790091 341.4132  760.4817      2
3  1  951.04189  1.6199491 930.3381 1619.4957      2
4  1  115.06607 -1.8569732 120.4632 2118.0770      1
5  1 1128.34326 -2.0536619 205.9932 2195.0499      2
6  1   54.65923 -0.4443021 485.0652 1101.7625      1

So the first state seems not to honor the given initial distribution, because the first sequence starts from state 2 instead of state 1.

A minor note: delta=c(0,1) in the previous example gives and error
Error in sample.int(length(x), size, replace, prob) : NA in probability **vector.**

The example given in the documentation of simData function seems to have the same issue, since one can get a realization where female begins in state 2, see below. (Note that I added "states=TRUE" to the original example.)

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

set.seed(123)
data.delta<-simData(nbAnimals=2,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par,
                    delta=delta,formulaDelta=formulaDelta,covs=cov, states=TRUE)    
head(data.delta)
  ID       step      angle        x         y sex states
1  1  833.60385         NA   0.0000    0.0000   F      2
2  1 1041.50736 -0.1790091 341.4132  760.4817   F      2
3  1  951.04189  1.6199491 930.3381 1619.4957   F      2
4  1  115.06607 -1.8569732 120.4632 2118.0770   F      1
5  1 1128.34326 -2.0536619 205.9932 2195.0499   F      2
6  1   54.65923 -0.4443021 485.0652 1101.7625   F      1

Sampling of the first state seems to happen on line 1043:
Z[1] <- sample(1:nbStates,size=1,prob=delta0%*%gamma)

What is confusing me is that why this initial distribution delta is multiplied with gamma when sampling the first state? What am I missing here? And is there any way to guarantee that the simulated sequence starts from a certain state?

Thanks in advance,
Anna-Kaisa

In applying momentuHMM on a dataset comprised of subtracks from several collared individuals, I have 6 random effect terms that I would like to include/test for: sex, age, body condition, antibody titre, group size, and presence of juveniles in group.

In attempting to test one of the above random effects, I am struggling with the actual application of training the model with fitHMM because I do not understand how to utilize the mixture feature or how to obtain initial parameter values for the regression coefficients of the transition probabilities (beta0 values) per mixture.

I would certainly appreciate your assistance with this and thank you in advance.

Hi Brett, I'd like to forbid the state sequence to start at a certain state, so fix it's delta at 0, while still allowing for the estimation of the initial probabilities of the other states. I tried fixPar = list(delta=c(NA,0,NA)), but I get the error of: fixPar$delta must sum to 1.
I also tried, without fixing anything, giving initial values like delta_0 <- c(0.70,0.001,0.299). It still gives me a probability close to 1 of starting at the state it shouldn't. Any thoughts?

• Rocío

I am working with telemetry data with GPS tags, which have very little location error, however for certain tracks I have very big time gaps (a few hours) between locations. During these big gaps dive start and end locations are interpolated between GPS locations.
It was suggested to me that to remove this uncertain locations from the model I could have introduced NAs for the step length and turning angle for these locations before fitting the model. For example, if an interpolated location didn’t have a real GPS location within a time window the step length and turning angle of that location would become NAs. This has helped a lot my models to classify different states, especially removing sections of the trips that resulted in straight lines due to dive location interpolation. I was wondering how exactly does the momentuHMM treats and takes into account these NAs?

Hi Brett,

Exploring CIreal(example$m, alpha = 0.95, parms = "gamma"), I observe that there are 4 possible transition-state probabilities:

1. To remain in the state 1 (1 -> 1).
2. To change from the state 1 to the state 2 (1 -> 2).
3. To change from the state 2 to the state 1 (2 -> 1).
4. To reman in the state 2 (2 -> 2).

Why, when I use CIbeta(example$m)$beta$est, Do I see only two beta estimates? One for the transition state 1 -> 2 and other for 2 -> 1. Why there are no beta estimates for the transition states 1 -> 1 and 2 -> 2?

I am interested in knowing if several covariates affect the transition-state probability of my target species, and I was wondering why there is no beta for the transition states 1 -> 1 and 2 -> 2.

I am attempting to follow Theo Michelot's tutorial from his seminar with the Ecological Forecasting Initiative on 2022 Feb 7. I am working with GPS collar data (only GPS, no Argos) for polar bears and I need to regularize the track to evenly spaced time intervals (60 min) and am doing so with the crawlWrap function. I played around with the initial parameters and think I have found some that make sense. The model output looks believable for a 2-state single track (but maybe I'm wrong, see below for workflow and model output):

Load the data

raw <- read.csv("PB21615.csv") #attached
PB21615.csv

head(raw)

keep relevant columns: animal, datetime, longitude, latitude, x, y (LAEA projection)

data_all <- raw[,c(1,4,7,8,22, 23)]
head(data_all)
colnames(data_all) <- c("ID", "time", "lat", "lon", "x", "y")
data_all$time <- as.POSIXct(data_all$time, format= "%m/%d/%Y %H:%M", tz="GMT")

data <- data_all

Remove duplicate datetimes

dups <- which(duplicated(data$time)) #none in this bear's data

Plot Northing vs Easting (in LAEA centered on Wrangel Island, Russia)

ggplot(data, aes(x, y, col = ID)) +
geom_point(size = 0.5) + geom_path() +
coord_equal()

Table of time intervals in data - regularize for one location every 60 minutes

plot(table(diff(data$time)), xlim = c(0, 60),
xlab = "time interval (min)", ylab = "count")

If have very long gaps then split the track into separate tracks, but for PB we have found that crawl predicts fairly well (with posterior simulation) for anything within 14 days

Use function from utility_function.R to split data at gaps > 14 days

data_split <- split_at_gap(data = data, max_gap = 33660, shortest_track = 2460)

Predict locations on 60-min grid using crawl (through momentuHMM wrapper)

ln.prior <- function(theta) dnorm(theta[2], -4,2,log=T)

crw_out <- crawlWrap(obsData = data_split,
timeStep = "60 min",
Time.name = "time",
coord = c("x", "y"),
theta = c(log(10),3),
fixPar=c(NA,NA),
prior=ln.prior,
initialSANN = list(maxit = 5000),
retryFits=10,
attempts=10)

crw_out1 <- crawlWrap(obsData = data_split, timeStep = "60 min",
Time.name = "time", coord = c("lon", "lat"), initialSANN = list(maxit = 5000), retryFits=10)

The output can then directly be passed to prepData to get step lengths and turning angles.

data_hmm <- prepData(data = crw_out)

MomentuHMM Workflow CRAWL

data preparation

data_hmm <- prepData(data = crw_out)

Figure out best initial parameters

Plot histogram of step lengths and turn angles

hist(data_hmm$step, xlim=c(0,300), breaks=100) #most step lengths are 20m or less
hist(data_hmm$angle, breaks=seq(-pi, pi, length=15))

Observation distribution (step lengths and turning angles)

dist <- list(step = "gamma", angle = "vm")

Par0_2s <- list(step = c(10, 100, 50, 30), angle = c(.03, 2))

Fit a 2-state HMM

states <- c("ARS/rest", "Traveling")
hmm1 <- fitHMM(data_hmm, nbStates = 2, dist = dist, Par0 = Par0_2s, stateNames=states)

print parameter estimates

hmm1

Plot estimated distributions and state-coloured tracks

plot(hmm1, breaks = 10, ask = T)

###############################################################
Then I attempt to use the multiple imputation function to account for the uncertainty in predicted locations. But, when I run the function MIfitHMM I get a convergence error (Hessian is singular - I have copied the errors below the code). I tried to fix this by adding additional arguments to the crawlWrap function (retryFits, attempts, and added theta, prior, fixPar) but I still get the same errors. In this example I am only running nSims=2 to decrease processing time, but even when I run with more simulations I still get the same errors (more NaNs produced). I'm am guessing I need to adjust something within the crawlWrap function but have yet to be successful. Thanks for the help!

Fit HMM using multiple imputation

library(mitools)
library(doFuture)
hmm6 <- MIfitHMM(miData = crw_out, nSims = 2, ncores = 1, nbStates = 2,
dist = dist, Par0 = Par0_2s, retryFits = 50)
hmm6

Errors

Warning messages:
1: In MIpool(fits, alpha = alpha, ncores = ncores, na.rm = na.rm) :
Hessian is singular for HMM fit(s): 1, 2
2: In sqrt(diag(x)) : NaNs produced
3: In MIpool(fits, alpha = alpha, ncores = ncores, na.rm = na.rm) :
working parameter standard errors are not finite for HMM fits 2
4: In sqrt(diag(x)) : NaNs produced

Hi,

I have a dataset from acoustically-tagged female Giant Mud Crab (Scylla serrata) in two rivers in mid-north New South Wales (Australia). This project aims to investigate the environmental triggers of their migration to sea to spawn. The rivers are very winding (think S-shaped bends) and some of our listening stations are placed on opposite sides of a bend.

My question is related to how suitable momentuHMM is to analyse this data given this. I have calculated the distance-to-sea of each listening station using a raster grid (size = 1m^2) and a 'least cost' path analysis. An initial idea I had was to use the difference in these distances (which account for the windiness of the river) to calculate acceleration of an individual and also the increase/decrease in their distance to sea as the data streams (instead of step length and turning angle) to infer their states.

Any advice you can provide is greatly appreciated!

Kind regards,

Dan

I am having trouble with crawlWrap function.
My data seems be in identical format to example data.
However, on execution of this line:

cr<-crawlWrap(sub[,c("ID","x","y","datetime")],attempts=10,
Time.name="datetime",coord=c("x","y"),timeStep="30 min",theta=c(10, .8),fixPar=c(NA,NA),ncores=3)

I receive the following error:

Error in crwData(list(crwFits = model_fits, crwPredict = predData)) :
Can't construct crwData object: fields are missing

This is my data
test.xlsx

This is my code:
#LOAD PACKAGE
library(momentuHMM)

#LOAD DATA
test<-read.table("test.txt",h=T)
head(test)

#CORRECT DATETIME
datetime<-as.POSIXct(paste(test$date,test$time),format="%d/%m/%Y %H:%M:%S")
datetime<-datetime + test$GMT6060
test<-cbind(test,datetime)
head(test)

#ORDER DATA
sub<-test[order(test$ID,test$datetime),]
head(sub)

#FIT CONTINUOUS MOVEMENT MODEL
cr<-crawlWrap(sub[,c("ID","x","y","datetime")],attempts=10,
Time.name="datetime",coord=c("x","y"),timeStep="30 min",theta=c(10, .8),fixPar=c(NA,NA),ncores=3)
plot(cr)

What could be causing this error?
Thank you.

ps: test.txt file is equal to test.xlsx, but in a different format. Please convert it to the txt format before run the code or, if you prefer, I can send you the test.txt file by email. Thank you very much.

Hello,

I am working on GPS data collected on sheep. Looking at the distribution of turning angles, I can tell that, at least for some behavioural states, the distribution is centred on PI and not 0.
To account for this effect, I need to specify to fitHMM that the mean for one state is PI (in my case, the first state).
fitHMM seems to only accept user-specified means if estAngleMean is set to TRUE for the relevant data stream ("angle" in my case).
I thus had to set the angle means as fixed parameters for both states.
Here is my code:

`dist <- list(step = "gamma", angle = "vm")

Par0_2s <- list(step = c(10, 130, 20, 50), angle = c(pi, 0, 0.1, 3))
fixPar_2s <- list(angle = c(pi, 0, NA, NA))

hmm1 <- fitHMM(data_hmm1, nbStates = 2, dist = dist, Par0 = Par0_2s, estAngleMean=list(angle=TRUE), fixPar = fixPar_2s)`

The optimisation runs without issue.
However, the computed parameter angle means are slightly different from PI and 0 (2.9345068 for state 1 and 0.06353729 for state 2).
More surprising, when I run the same code without the fixPar parameter, I get different values (3.0594606 and 0.0522979).
Thus, it seems that the specified fixed parameters are taken into account, but not to the point that the parameters are actually fixed.

How can I keep the angular means for both states actually fixed during the optimisation and still have a mean different from 0 for some states?

I’m using momentuHMM version 1.5.4 in R version 4.1.3.

Thank you in advance for your time.

I am using momentuHMM package to work with coyote movement data that is irregularly sampled. I have had no problem using crawlWrap, prepData, and fitHMM to generate HMM using single imputation but I am having trouble implementing multiple imputation using MIfitHMM. Below is my code, and the two types of errors I have received.

crwOut <- crawlWrap(data, timeStep = "2 hours", theta=c(6.855, -0.007), fixPar=c(NA,NA))
CBcoyData <- prepData(data=crwOut)

dist <- list(step = "gamma", angle = "wrpcauchy")
Par0_m2 <- list(step=c(50,700,50,350),angle=c(3,0.3,0.3,0.5))

m2 <- fitHMM(data = CBcoyData, nbStates = 2, dist = dist, Par0 = Par0_m2, estAngleMean = list(angle=TRUE)) ### single imputation

nfsFits2 <- MIfitHMM(crwOut, nSims = 2, nbStates = 2, dist = dist, Par0 = Par0_m2, estAngleMean = list(angle = TRUE), retryFits = 30)) ### two imputations for testing

In MIfitHMM(crwOut, nSims = 2, nbStates = 2, dist = dist, Par0 = Par0_m2, :
prepData failed for imputation 1; Error in eval(.doSnowGlobals$expr, envir = .doSnowGlobals$exportenv):
crawl::crwPostIS error for individual 26.01; Error in sample.int(length(x), size, replace, prob): NA in probability vector
Check crwPostIS arguments, crawl::crwMLE model fits, and/or consult crawl documentation.

prepData failed for imputation 2; Error in eval(.doSnowGlobals$expr, envir = .doSnowGlobals$exportenv):
crawl::crwPostIS error for individual 26.01; Error in sample.int(length(x), size, replace, prob): NA in probability vector
Check crwPostIS arguments, crawl::crwMLE model fits, and/or consult crawl documentation.

My first thought was these errors might be linked to high missing rate or situations where a single fix is surrounded by a series of missing fixes. So I removed the individual listed above from data, then reran. I received the error above twice more, each time with a new offending individual (that I then removed). Then I got this error:

Drawing 2 realizations from the position process using crawl::crwPostIS...
Error in { : task 20 failed - "missing value where TRUE/FALSE needed"

I then looked at all individuals and removed 13 more individuals that had single fixes surrounded by series of missing fixes. Then I started getting this type error again:

prepData failed for imputation X; Error in eval(.doSnowGlobals$expr, envir = .doSnowGlobals$exportenv):
crawl::crwPostIS error for individual XXXX; Error in sample.int(length(x), size, replace, prob): NA in probability vector

I removed a few more individuals based on this error, but I now I am down to half my initial sample size and beginning to wonder if there is some other issue causing the errors.

Note I have also been successful with single imputation using 3 and 4 movement states but again was unsuccessful implementing multiple imputation (with the same error codes).

Is there something wrong with my data or code?

When working through the grey seal example, running the crawlWrap() function

crwOut<-crawlWrap(greySealData,err.model=list(x=~ln.sd.x-1,y=~ln.sd.y-1,rho=~error.corr),
                   fixPar=c(1,1,NA,NA),attempts=100,retryFits=20,
                   timeStep="2 hours")

Fitting 1 track(s) using crawl::crwMLE... 

returns the error:

Error in get(s, env, inherits = FALSE) : 
  argument "theta" is missing, with no default

Hello,
All the examples of using momentuHMM I have seen online have used GPS data. Is it possible to use this package in the same way but on geolocator (GLS) data (which can have high error)? I have seabird GLS data, so two location estimates every day and I'd like to try and identify migration and breeding periods. Any help would be really appreciated.
Thank you!

Hi Brett,

I have a 2-state ("1. Active/2. Resting") HMM in which, among other things, I want to include the effect of the water temperature on the parameters of the state-dependent distributions of meanActivity, which is my data stream. Similarly to what you did in the Northern fur seal example, I want to constrain the Gamma parameters such that the "Active" state tends to have higher activity values than "Resting". I guess I have to do this using DM and workBounds arguments. However, I do not know how to call DM such that I include the effect of the water temperature and force meanActivity[state=="Active"] > meanActivity[state=="Resting"].

Similarly to my case, in the Northern fur seal example it would be like including water temperature as a covariate in the HMMs to see the effect of water temperature on the parameters of the state-dependent distributions of "step", "angle" and "dive", while keeping the constraints you mentioned. If I am not wrong, in the vignettes there is no example of doing both things at once. Is there a reason?

I show the code I was able to run so far adaptating what you did in the Northern fur seal :

df <- df[,c("ID","meanActivity","WaterT")]
df <- prepData(data=df,coordNames=NULL,covNames = c("WaterT"))
nbStates <- 2
stateNames <- c("Active","Resting")
dist = list(meanActivity = "gamma")

Par0_m1 <- list(meanActivity=c(1.2,0.15,1.2,0.15))

meanActivityDM<-matrix(c(1,0,0,0,
                 1,1,0,0,
                 0,0,1,0,
                 0,0,0,1),2*length(stateNames),4,byrow=TRUE,
               dimnames=list(c(paste0("mean_",1:length(stateNames)),paste0("sd_",1:length(stateNames))),
                             c("mean_12:(Intercept)","mean_2",
                               paste0("sd_",1:length(stateNames),":(Intercept)"))))

meanActivityworkBounds <- matrix(c(-Inf,-Inf,-Inf,-Inf,Inf,0,rep(Inf,2)),4,2,
                         dimnames=list(colnames(meanActivityDM),c("lower","upper")))

DM<- list(meanActivity=meanActivityDM)

workBounds<- list(meanActivity=meanActivityworkBounds)

Par0 <- getParDM(nbStates = nbStates, dist = dist,                
                 Par = Par0_m1, DM = DM, workBounds = workBounds,    
                 estAngleMean = list(angle = FALSE))

formula <- ~ MWTB

m3 <- fitHMM(data = HJose.InviernoData, nbStates = 2, dist = dist, Par0 = Par0, 
            DM = DM, workBounds = workBounds, stateNames = stateNames, formula = formula)

This works. However, if I do not know how to include WaterT effect on DM. I tried different things:

DM<- list(meanActivity=meanActivityDM~WaterT)
DM<- list(meanActivity=list(meanActivityDM=~WaterT))
Etc.

What am I doing wrong? If my question is inappropriate I apologize in advance.

Hi Brett,

I've come back to some old code where I'm modelling step length, pitch and yaw and my previous code doesn't seem to be working anymore. It works fine if I model just step and yaw, but it really doesn't like pitch, and won't model it even on its own. The data for pitch seem fine, they are well within -pi and pi.

The error I get is this:

> mNULL_exp1_2st <- momentuHMM:::fitHMM(data = exp1prep, nbStates=2, 
                                   dist = list(step = 'gamma', pitch = 'wrpcauchy', yaw = 'wrpcauchy'),
                                   Par0 = list(step = stepPar0, pitch = PanglePar0, yaw = YanglePar0),
                                   stateNames = stateNames)
Error in get_fixParIndex(Par0, beta0, delta0, fixPar, distnames, inputs,  : fixPar$pi must be of length 1

Any idea what this error message means?

The reason I came back to this code is because I have some new data to feed through the analysis, but I can't get the code to run with the old data either. I updated R and downloaded the latest version of momentuHMM but I can't figure out what has changed and what I need to do differently now.

Any help would be hugely appreciated!

Many thanks,
Theoni

Hi Brett,

I am using your adaption of Enrico Pirotta's 6-state HMM (Northern Fulmar example) on seabird data and have since adapted the code to add an additional state. We are now looking at the effects of sex, tagging year etc on the transition probabilities. Is it possible to extract the transition probabilities and their confidence intervals from the models to look at the effect of covariates, rather than just the means that are shown when you print the model? I see that these are calculated and plotted within 'plotTPM', but with so many states and the covariates it would be great to extract the values. Is there a function that I have overlooked for extracting these values?

All the best,
Gavin

Hi Brett,

Is there an in-built function in the momentuHMM package (or if not, a workaround) to plot the estimated covariate effects on the transition probabilities? In other words, a plot of transition probability (y-axis) versus covariate effect (x-axis). Ideally I would like to generate predictions of the estimate +- 95% confidence intervals. And if possible, generate separate plots for different mixtures.

Thanks

Hi Brett and Theo! I love your package, and we are just doing a workshop on how to use it with marine data, see: https://github.com/MarieAugerMethe/DFO_HMM_2022

I have two little wishes that if you have time one day would be great to implement.

• Would it be possible to have the preData function have a option to place NAs in the response variable when we are missing a data point (e.g. we have regular time series with locations at 8:00, 8:15, 8:30, 9:00; could we get a NA place in the step length and turning angle or dive at 8:45).
• Would it be possible to have an option for parallel computing for fitHMM. It would be particularly useful for the retryFit option.

Thanks a million for the wonderful package!

Hi Brett,

I wish to include knownStates in my HMM. My model has three states: transiting, foraging and resting. I have determined all of the time periods when my animals are assumed to be "resting" (- the data are from seals and following McClintock et al., 2013 and Russell et al., 2015 I have classified proportion of time spent diving using <0.444 threshold).

I hence wish to fit a 3-state model where all of the time intervals representing "resting" states are "known", but allow the model to determine the remaining time intervals as either travelling and foraging states.

I have tried fitting the model using the knownStates argument, but when I do this the model not only includes the time intervals that I have specified as known states, but it also fits many more time intervals as additional resting states.

Is there a way to force the model to only assign either travelling or foraging to the remaining time intervals that are not specified in knownStates?

This analysis is essentially the same as what you did in Carter et al., 2020 Oikos: "We fitted a three state multivariate, discrete-time HMM with resting as a ‘known state’, and attributed the remaining time intervals with low step lengths and high turn angles to foraging, and intervals with high step lengths and low turn angles to travelling."

The code I am working with is provided below.

Thank you very much,

Will

nbStates<-3 # 3 States (1 = "Transiting", 2 = "Foraging", 3 = "Resting")

dist<-list(step="gamma", 
             angle="wrpcauchy") 

estAngleMean<-list(angle=TRUE)

circularAngleMean<-list(angle=TRUE) 

DM<-list(step= list(mean=~HW*sex*site,
                               sd=~HW), 
                angle=list(mean=~1, 
                                 concentration=~HW*sex*site)) 

formula<-~HW*sex*site

knownStates<-mHMM$knownStates

checkPar0(mHMM,nbStates=nbStates,
          dist=dist,
          DM=DM,
          mixtures=3,
          knownStates=knownStates,
          estAngleMean=estAngleMean,
          larAngleMean=circularAngleMean)

getPar<-list(step=c(0.89,0.64,0.30,0.13,0.15,0.19),
             angle=c(0,0,0,0.6,0.1,0.05))

(getPar<-getParDM(data=mHMM,
                  nbStates=nbStates,
                  dist=dist,
                  Par=getPar,
                  estAngleMean=estAngleMean,
                  knownStates=knownStates,
                  circularAngleMean=circularAngleMean,
                  DM=DM))

Par0<-list(step=getPar$step,
           angle=getPar$angle)

stateNames<-c("Transiting","Foraging","Resting")

m25k<-fitHMM(data=mHMM,
             nbStates=nbStates,
             dist,
             Par0=Par0,
             DM=DM,
             formula=formula,
             mixtures=3,
             knownStates=knownStates,
             estAngleMean=estAngleMean,
             circularAngleMean=circularAngleMean,
             stateNames=stateNames)