Building upon simpler models using getPar0

[CatTFonseca]

Hello,

I am using momentuHMM to analyse the diving behaviour of fin whales, following the example provided on pilot whales.
I already ran my models with no discrete effects or covariates. I ran with 30 different initializations for the initial values (Par0), tested with 2 to 6 hidden states, and then chose the best model based on the negative log-likelihood and AIC values.

Now I would like to introduce random effects to check their impact and I am wondering what is the best approach to do that. I have seen in previous examples that I should build upon previous models.

1. Then, should I use getPar0 to get starting values from each of the previous models (with mixtures = 2 for example) and then fit an HMM again using those initial values and mixtures equal to 2 (code example below)? getPar0 initializes some unknown values as 0, can I use them like that or will that negatively impact my results?

# Number of existing models to use to get starting values
nmods = length(all_mods_unlist)

all_mods_mix <- list() # empty list to save all models
all_params_mix <- list()

for(i in 1:nmods){
    ## get starting values for new model from existing momentuHMM model fit
    Par0_mix <- getPar0(all_mods_unlist[[i]], mixtures=2)
    all_params_mix[[i]] <- Par0_mix
    
    ## fit model with 2 mixtures on TPM
    all_mods_mix[[i]]  <- fitHMM(processedDive_dat, 
                                 nbStates=length(all_mods_unlist[[i]]$stateNames), 
                                 dist=dist,
                                 Par0=Par0_mix$Par, beta0=Par0_mix$beta,
                                 stationary=FALSE,
                                 mixtures=2,
                                 #stateNames=stateNames,
                                 )
    md <- all_mods_mix[[i]] 
  }
}

2. Or instead, should I just build upon the 'best model' and try it with different values of mixtures?

I am also trying to test with individual level fixed effects but it is taking ages to run it only using the 'best model', like this:

#--- With individual level fixed effects ---#
Par0_fix <- getPar0(mbest,formula=~0+ID,mixtures=1)
fitfix <- fitHMM(processedDive_dat, 
                 nbStates=6, 
                 dist=dist,
                 formula = ~0+ID,
                 Par0=Par0_fix$Par, 
                 beta0=Par0_fix$beta,
                 #stateNames=stateNames,
                )
fitfix

[bmcclintock]

This can be as much art as science for complicated models. Either option might work well. I'd suggest you try both and also look into the retryFits argument. The key is to monitor the negative log likelihood value, where lower is better (this is printed when using the retryFits argument or can be accessed via model$mod$minimum (e.g. fitfix$mod$minimum). If you do a bunch of retryFits and at some point the negative log likelihood stops getting lower (and fairly consistently converges to this value thereafter), then this can suggest convergence (but not necessarily to the global maximum). See this {moveHMM} vignette for more on the topic. Also see this paper on the dangers of false convergence for these types of models.

For the individual-level fixed effects model, one trick is to first fit each track individually and then combine these to get starting values for the individual-level fixed effects model (the beta estimates shouldn't change much from the individual model fits). See the harbour seal example in the {momentuHMM} vignette and the accompanying code for the individual-level effects model.

[CatTFonseca]

Thanks for the help! Will try it out.