A function for running momentuHMM

[EliGurarie]

getHMMStates <- function(df, saveme = FALSE){

dir <- "../YNS_RSF/results/hmm/"
filename <-  with(df[1,], paste0(id, season, year, ".rda"))

if(!(filename %in% list.files(dir)) & nrow(df) > 30){

  with(df[1,], cat(c(id, year, season, "\n")))
  df <- df %>% mutate(z = x + 1i*y, z.mod = c(NA, Mod(diff(z))),
                      dt = c(NA, difftime(time[-1], time[-length(time)], unit = "hours"))) 
  
  df.prep <- prepData(df %>% as.data.frame, type = "UTM",
                      coordNames = c("x", "y"))  %>%  
    mutate(step = z.mod/dt, 
           step = ifelse(dt > 5.5, NA, step),   # NEED TO ALSO FILTER OUT VERY SMALL TIMESTEPS
           angle = ifelse(dt > 5.5, NA, angle))
  
  df.prep$step[df.prep$step == 0] <- 1
  
  {
    stepMean0 <- c(m1 = 100, m2 = 4000)
    stepSD0 <- c(sd1 = 100, sd2 = 1000)
    angleCon0 <- c(rho1  = 0.1, rho2 = .5)
  
    stateNames = c("resident","transit")
    
    dist = list(step = "gamma", angle = "wrpcauchy")
    Par0 = list(step=c(stepMean0, stepSD0),
                angle=c(angleCon0))
    
    df.fit = fitHMM(data = df.prep, nbStates = 2, dist = dist, 
                    Par0 = Par0, stateNames = stateNames)
    
  }
  state <- try(viterbi(df.fit))
  if(inherits(state, "try-error")) state <- NA
  
  M <- df.fit$mle$gamma
  step <- df.fit$mle$step
  angle <- df.fit$mle$angle[2,] 
  
  probs = c(s2 = M[1,2]/(M[1,2] + M[2,1]), 
            s1 = M[2,1]/(M[1,2] + M[2,1]))
  
  r <- list(state = state, M = M, probs = probs, step = step, angle = angle)
  if(saveme) save(r, file = paste0(dir, filename))
  return(r)
}

}