TheoMichelot / moveHMM

R package - Animal movement modelling using hidden Markov models
https://cran.r-project.org/package=moveHMM
38 stars 11 forks source link

Feature request: use fitted model to decode new data #30

Closed FlukeAndFeather closed 2 years ago

FlukeAndFeather commented 3 years ago

With very large data sets, fitting the HMM can take a prohibitively long time. Is it possible to fit a HMM to a subset of the data, then decode the remainder using that fitted model? Something like the following. I think the implementation would be fairly straightforward; add a newdata parameter to viterbi() (replacing viterbi.R#L30).

library(dplyr)
library(moveHMM)

# Simulate data
nbAnimals <- 10
nbStates <- 2
nbCovs <- 2
mu <- c(15, 50)
sigma <- c(10, 20)
angleMean <- c(pi, 0)
kappa <- c(0.7, 1.5)
stepPar <- c(mu, sigma)
anglePar <- c(angleMean, kappa)
stepDist <- "gamma"
angleDist <- "vm"
zeroInflation <- FALSE
obsPerAnimal <- 100

data <- simData(nbAnimals = nbAnimals,
                nbStates = nbStates,
                stepDist = stepDist,
                angleDist = angleDist,
                stepPar = stepPar,
                anglePar = anglePar,
                nbCovs = nbCovs,
                zeroInflation = zeroInflation,
                obsPerAnimal = obsPerAnimal)

# Split train-decode
train_subset <- filter(data, ID %in% unique(data$ID[1:5]))
decode_subset <- anti_join(data, train_subset, by = "ID")

# Fit model
m <- fitHMM(data = train_subset, 
            nbStates = nbStates, 
            stepPar0 = stepPar,
            anglePar0 = kappa,
            stepDist = stepDist,
            angleDist = angleDist,
            angleMean = angleMean)

# This would be useful!
# viterbi(m, decode_subset)
dibeyp commented 3 years ago

Would probably be useful to have smoothed backwards probability as well as viterbi states.

Just a thought Cheers Toby

On Wed, Nov 10, 2021, 05:03 Max Czapanskiy @.***> wrote:

With very large data sets, fitting the HMM can take a prohibitively long time. Is it possible to fit a HMM to a subset of the data, then decode the remainder using that fitted model? Something like the following. I think the implementation would be fairly straightforward; add a newdata parameter to viterbi() (replacing viterbi.R#L30 https://github.com/TheoMichelot/moveHMM/blob/7253976472fc1a98d8b594e8aab9a4fa8af960bb/R/viterbi.R#L30 ).

library(dplyr) library(moveHMM)

Simulate datanbAnimals <- 10nbStates <- 2nbCovs <- 2mu <- c(15, 50)sigma <- c(10, 20)angleMean <- c(pi, 0)kappa <- c(0.7, 1.5)stepPar <- c(mu, sigma)anglePar <- c(angleMean, kappa)stepDist <- "gamma"angleDist <- "vm"zeroInflation <- FALSEobsPerAnimal <- 100

data <- simData(nbAnimals = nbAnimals, nbStates = nbStates, stepDist = stepDist, angleDist = angleDist, stepPar = stepPar, anglePar = anglePar, nbCovs = nbCovs, zeroInflation = zeroInflation, obsPerAnimal = obsPerAnimal)

Split train-decodetrain_subset <- filter(data, ID %in% unique(data$ID[1:5]))decode_subset <- anti_join(data, train_subset, by = "ID")

Fit modelm <- fitHMM(data = train_subset,

        nbStates = nbStates,
        stepPar0 = stepPar,
        anglePar0 = kappa,
        stepDist = stepDist,
        angleDist = angleDist,
        angleMean = angleMean)

This would be useful!# viterbi(m, decode_subset)

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/TheoMichelot/moveHMM/issues/30, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACORVDTVJEURTV6EEJRPMD3ULFO7XANCNFSM5HV5ROCA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

TheoMichelot commented 2 years ago

There is some support for what you are asking, I think, with the fit = FALSE argument of fitHMM. When that argument is used, a model is created with the parameters fixed to the initial values provided (rather than estimated from the data). Then, you can call viterbi on that object to get the output you're interested in.

In your example, you would need the following after fitting the first model:

m2 <- fitHMM(data = decode_subset,
             nbStates = nbStates,
             stepPar0 = as.vector(t(m$mle$stepPar)),
             anglePar0 = as.vector(t(m$mle$anglePar))[3:4],
             beta0 = m$mle$beta,
             stepDist = stepDist,
             angleDist = angleDist,
             angleMean = angleMean,
             fit = FALSE)

viterbi(m2)

One advantage of this approach is that you can also call other moveHMM functions on this new model, e.g., stateProbs (as suggested by Toby) or plotPR.

On the other hand, your suggestion for viterbi doesn't require much additional code, and might be more intuitive to some users, so I think I will integrate it. Thanks!