Fitting a HMM with multiple known states and TPM constraints

[MattIDCarter]

Hi Brett, I hope you're doing well!
I'm working on fitting HMMs to seal tracking data. The basic idea is the same as what we did in the 2020 Oikos paper without the covariate effects on state-dependent parameters... i.e. a known state "rest" determined by the proportion of time diving, and forage and transit states determined by the model.
I have had success fitting this in momentuHMM and constraining the TPM by setting beta0 accordingly to avoid transition into or out of the known state when this state is not active.
This is all pretty similar to Will's question here: https://github.com/bmcclintock/momentuHMM/discussions/68

However, I want to extend this model a little and add complexity. I would like to break the "rest" state down into "haulout" and at-sea rest ("surface") based on majority rule for time spent in haulout. So, effectively a 4 state model with two known states: haulout (state1) and surface (state2), as well as ARS (state3) and transit (state4). Of course I could just fit one known state and assign to at-sea or haulout post-hoc, but I am interested in the transition probabilities from the other states, and whether or not distinguishing between these two resting states may impact the inference of ARS and transit states.

I'm having a little trouble parameterising beta0 to constrain the TPM and avoid switching into or out of either known state when they are not active (i.e. I don't want the model to assign state 1 when haulout == 0, and likewise it should not assign state 2 when surface == 0). Playing around with it, it seems beta0 expects a 3-row matrix, with the first row being (Intercept), the second being haulout1, and the third being surface1. I have been working on the assumption that (Intercept) is when neither known state is active, but I'm not sure this is correct. I am still getting instances where transitions into haulout or surface are happening from either ARS or transit when the known states are not active.

I note that in your response to Will, you suggest setting workbounds as extra insurance... I have not found this to be necessary in my case to get the desired effect with one known state, however I thought it may help here. I tried to set this up but it seems the model wants a delta0 object of length 9, which confused me somewhat as I'm not sure what it's expecting; 9 would imply a matrix of 3x3... so I'm not sure what the columns and rows should correspond to.

Any advice greatly appreciated!
Cheers,
Matt

Code below:

nbStates <- 4
stateNames <- c("haulout", "surface", "ARS", "travel")  
stepDist <- "gamma"
angleDist <- "wrpcauchy"
nbStates <- 4

stepPar0 <- c(0.05, 0.5, 0.5, 1, 0.01, 0.2, 0.2, 1, 0.03, 1e-10, 1e-10, 1e-10)
anglePar0 <- c(0.99, 0.7, 0.4, 0.95)

knownStates <- rep(NA, nrow(d))
knownStates[d$haulout == 1] <- 1
knownStates[d$surface == 1] <- 2

formula <- ~ 0 + as.factor(haulout) + as.factor(surface)

betaRef <- c(3, 3, 3, 4)
beta0 = matrix(c(-100, -100, NA, -100, -100, NA, -100, -100, NA, -100, -100, NA,
                 100, -100, -100, 100, 0, 0, 100, 0, 0, 100, 0, 0,
                 0, 100, 0, -100, 100, -100, 0, 100, 0, 0, 100, 0), nrow = 3, ncol = nbStates*(nbStates-1), byrow=TRUE,
               dimnames = list(c("(Intercept)", "haulout1", "surface1"),
                               c("1 -> 1", "1 -> 2", "1 -> 4", "2 -> 1", "2 -> 2", "2 -> 4", "3 -> 1", "3 -> 2", "3 -> 4", "4 -> 1", "4 -> 2", "4 -> 3")))
fixPar <- list(beta=beta0)

test <- fitHMM(d, nbStates = nbStates, dist = list(step = stepDist, angle = angleDist), 
               Par0 = list(step = stepPar0, angle = anglePar0), 
               formula = formula, betaRef = betaRef, fixPar = fixPar,
               knownStates = knownStates, stateNames = stateNames)

[bmcclintock]

You're close, but the -100s should only be in the first row for beta0. Here's an extension of the #68 example with two known states, I think it will do what you want but please be sure to play around with it and verify:

#### Simulate some data
nbObs <- 1000
nbStates <- 4
stateNames <- c("haulout","surface","foraging","transit")
dist <- list(step="gamma",angle="vm",dive="beta",dry="beta")

Par <- list(step=c(1,1,1,10,0.5,0.5,0.5,2),
            angle=c(0,0,0,0,0.1,0.1,0.1,10),
            dive=c(1,1,10,10,10,10,1,1),
            dry=c(10,1,1,1,1,10,10,10))

beta <- matrix(-1.5,1,nbStates*(nbStates-1))

simDat <- simData(obsPerAnimal = nbObs,nbStates=nbStates,dist=dist,Par=Par,beta=beta,stateNames=stateNames,states=TRUE)

fitDat <- simDat[-nrow(simDat),] # can't have missing covariates
threshold <- 0.44
fitDat$rest <- factor(ifelse(fitDat$dive<threshold,ifelse(fitDat$dry > threshold,1,2),0),levels=c(0,1,2))

fitDat$dive <- fitDat$dry <- NULL

#### Set up model
betaRef <- c(3,3,3,4) # make state 3 (foraging) the tpm reference for state 1 (haulout) and state 2 (surface)
beta0 <- matrix(c(-100, -100,  0.5, -100, -100,  0.5, -100, -100, -1.5, -100, -100, -1.5,
                   100,    0,    0,  100,    0,    0,  100,    0,    0,  100,    0,    0,
                     0,  100,    0,    0,  100,    0,    0,  100,    0,    0,  100,    0),3,nbStates*(nbStates-1),byrow=TRUE,dimnames=list(c("rest0","rest1","rest2"),c("1 -> 1", "1 -> 2","1 -> 4", "2 -> 1", "2 -> 2","2 -> 4","3 -> 1", "3 -> 2", "3 -> 4","4 -> 1", "4 -> 2", "4 -> 3")))

delta0 <- matrix(c(   0,   100,  100,
                   -100,  -100, -100,
                    100,     0,    0),3,nbStates-1,byrow=TRUE,dimnames=list(c("rest0","rest1","rest2"),stateNames[-1]))

formula <- ~ 0 + rest
formulaDelta <- ~ 0 + rest

#### check tpm
trProbs <- getTrProbs(fitDat[,c("ID","rest"),drop=FALSE],nbStates=nbStates,beta=beta0,formula=formula,betaRef=betaRef,stateNames=stateNames)
trProbs[,,which.min(fitDat$rest)] # tpm when next step has rest=0
trProbs[,,which(fitDat$rest==1)[1]] # tpm when next step has rest=1
trProbs[,,which.max(fitDat$rest)] # tpm when next step has rest=2

Par0 <- Par[c("step","angle")]
Par0$angle <- Par$angle[-(1:nbStates)]

#### Set knownStates and fixPar$beta for extra insurance (for even extra insurance could also set workBounds$beta to be [-Inf,<<100] or [>>-100,Inf] as appropriate for estimated tpm parameter)
knownStates <- rep(NA,nbObs-1)
knownStates[fitDat$rest==1] <- 1
knownStates[fitDat$rest==2] <- 2
fixPar <- list(beta=beta0)
fixPar$beta[which(beta0 > -100 & beta0 < 100 & beta0!=0)] <- NA
workBounds <- list(delta=matrix(c(-Inf,-Inf,100-1.e-6,100-1.e-6,-Inf,-Inf,100-1.e-6,-Inf,-Inf,Inf,-100+1.e-6,Inf,Inf,-100+1.e-6,Inf,Inf,-100+1.e-6,Inf),nrow=length(delta0),ncol=2))
fixPar$delta <- momentuHMM:::nw2w(delta0,workBounds$delta)
if(fitDat$rest[1]==0){
  fixPar$delta[1,] <- NA
}

#### fit the model
fit <- fitHMM(fitDat,nbStates=nbStates,dist=dist[c("step","angle")],Par0=Par0,beta0=beta0,formula=formula,formulaDelta=formulaDelta,betaRef=betaRef,fixPar=fixPar,workBounds=workBounds,knownStates=knownStates,stateNames=stateNames)
fit
plot(fit,plotCI=TRUE,ask=FALSE)
unique(viterbi(fit)[which(fitDat$rest==0)])
unique(viterbi(fit)[which(fitDat$rest==1)])
unique(viterbi(fit)[which(fitDat$rest==2)])

[MattIDCarter]

Thanks Brett- worked perfectly!
I had been setting up the known states as two separate covariates rather than one - that's why I couldn't figure out the structure of the beta0 matrix.
Thanks again for your help!
Matt