moveHMM-starting-values.R
In moveHMM: Animal Movement Modelling using Hidden Markov Models

## ----local-max, echo = FALSE, out.width='.49\\linewidth', fig.width=5, fig.height=5, fig.show='hold', fig.cap = 'Example of convergence issue. The optimisation was started from two different points, and the red line shows the steps of the optimiser in both cases. On the left, the starting value was too far from the maximum likelihood estimate, and the optimiser got stuck in a local maximum of the function.'----
# Likelihood function
f <- function(x) {
    res <- dnorm(x, 0, 1) + dnorm(x, 5, 0.5)
    # cat(x, ", ", res, ",\n", sep = "") # uncomment to get optimiser steps
    return(res)
}

# Negative likelihood function (for numerical minimisation)
g <- function(x) {
    return(- f(x))
}

# First scenario: starting value = - 3
x0 <- -3
m1 <- nlm(f = g, p = x0)

# Second scenario: starting value = 7
x0 <- 7
m2 <- nlm(f = g, p = x0)

# Optimiser steps in first scenario
steps1 <- matrix(c(-3, 0.004431848,
                   -3, 0.004431848,
                   -2.999997, 0.004431888,
                   -2.986704, 0.004611786,
                   -2.986701, 0.004611827,
                   -2.97293, 0.004805012,
                   -2.972927, 0.004805054,
                   -2.958645, 0.005012956,
                   -2.958642, 0.005013,
                   -2.943814, 0.005237254,
                   -2.943811, 0.005237299,
                   -2.928396, 0.005479779,
                   -2.928393, 0.005479826,
                   -2.912349, 0.005742692,
                   -2.912346, 0.005742741,
                   -2.895624, 0.006028491,
                   -2.895621, 0.006028541,
                   -2.878168, 0.006340079,
                   -2.878165, 0.006340132,
                   -2.85992, 0.00668085,
                   -2.859917, 0.006680905,
                   -2.840813, 0.007054788,
                   -2.840811, 0.007054845,
                   -2.820772, 0.007466599,
                   -2.820769, 0.007466658,
                   -2.79971, 0.007921875,
                   -2.799707, 0.007921937,
                   -2.777531, 0.008427303,
                   -2.777528, 0.008427368,
                   -2.754124, 0.008990936,
                   -2.754121, 0.008991005,
                   -2.729362, 0.009622544,
                   -2.729359, 0.009622616,
                   -2.703098, 0.01033407,
                   -2.703096, 0.01033415,
                   -2.675164, 0.01114025,
                   -2.675162, 0.01114033,
                   -2.645362, 0.01205942,
                   -2.64536, 0.01205951,
                   -2.613461, 0.01311463,
                   -2.613458, 0.01311472,
                   -2.579186, 0.01433518,
                   -2.579183, 0.01433528,
                   -2.542213, 0.01575873,
                   -2.54221, 0.01575883,
                   -2.502151, 0.01743428,
                   -2.502148, 0.01743439,
                   -2.458527, 0.01942651,
                   -2.458525, 0.01942663,
                   -2.410766, 0.02182201,
                   -2.410764, 0.02182214,
                   -2.358159, 0.0247385,
                   -2.358156, 0.02473864,
                   -2.299821, 0.02833869,
                   -2.299819, 0.02833884,
                   -2.234647, 0.03285136,
                   -2.234645, 0.03285152,
                   -2.161236, 0.03860365,
                   -2.161234, 0.03860383,
                   -2.077804, 0.04607092,
                   -2.077802, 0.04607112,
                   -1.982077, 0.05595239,
                   -1.982075, 0.05595261,
                   -1.871175, 0.06928082,
                   -1.871173, 0.06928107,
                   -1.741539, 0.08756124,
                   -1.741537, 0.0875615,
                   -1.589047, 0.112875,
                   -1.589046, 0.1128753,
                   -1.409683, 0.1477044,
                   -1.409682, 0.1477047,
                   -1.201467, 0.1938444,
                   -1.201466, 0.1938446,
                   -0.9685691, 0.2495736,
                   -0.9685681, 0.2495738,
                   -0.7268399, 0.3063317,
                   -0.7268389, 0.306332,
                   2.094733, 0.04447224,
                   0.2689029, 0.3847764,
                   0.2689039, 0.3847763,
                   -0.04701344, 0.3985016,
                   -0.04701244, 0.3985017,
                   0.001419394, 0.3989419,
                   0.001420394, 0.3989419,
                   -2.02258e-06, 0.3989423,
                   -1.02258e-06, 0.3989423),
                 ncol = 2, byrow = TRUE)

# Optimiser steps for second scenario
steps2 <- matrix(c(7, 0.0002676605,
                   7, 0.0002676605,
                   7.000007, 0.0002676455,
                   6.997859, 0.0002722824,
                   6.997866, 0.0002722672,
                   6.995683, 0.0002770557,
                   6.99569, 0.0002770403,
                   6.993471, 0.0002819877,
                   6.993478, 0.000281972,
                   6.991223, 0.0002870862,
                   6.99123, 0.0002870702,
                   6.988936, 0.0002923595,
                   6.988943, 0.0002923433,
                   6.98661, 0.0002978165,
                   6.986617, 0.0002978,
                   6.984244, 0.0003034667,
                   6.984251, 0.0003034498,
                   6.981835, 0.0003093201,
                   6.981842, 0.000309303,
                   6.979383, 0.0003153877,
                   6.97939, 0.0003153703,
                   6.976886, 0.0003216811,
                   6.976893, 0.0003216633,
                   6.974343, 0.0003282126,
                   6.97435, 0.0003281945,
                   6.971751, 0.0003349957,
                   6.971758, 0.0003349773,
                   6.969109, 0.0003420447,
                   6.969116, 0.000342026,
                   6.966415, 0.0003493752,
                   6.966422, 0.000349356,
                   6.963667, 0.0003570037,
                   6.963674, 0.0003569841,
                   6.960863, 0.0003649482,
                   6.96087, 0.0003649283,
                   6.958, 0.0003732282,
                   6.958007, 0.0003732078,
                   6.955077, 0.0003818646,
                   6.955084, 0.0003818439,
                   6.952091, 0.0003908804,
                   6.952098, 0.0003908591,
                   6.949039, 0.0004003001,
                   6.949046, 0.0004002784,
                   6.945918, 0.0004101507,
                   6.945925, 0.0004101285,
                   6.942726, 0.0004204615,
                   6.942733, 0.0004204388,
                   6.939458, 0.0004312643,
                   6.939465, 0.0004312411,
                   6.936113, 0.0004425942,
                   6.93612, 0.0004425704,
                   6.932685, 0.0004544893,
                   6.932692, 0.000454465,
                   6.929172, 0.0004669916,
                   6.929179, 0.0004669667,
                   6.925568, 0.0004801472,
                   6.925575, 0.0004801216,
                   6.92187, 0.0004940069,
                   6.921877, 0.0004939806,
                   6.918073, 0.0005086267,
                   6.918079, 0.0005085997,
                   6.91417, 0.0005240688,
                   6.914177, 0.0005240411,
                   6.910158, 0.0005404021,
                   6.910165, 0.0005403736,
                   6.906029, 0.0005577032,
                   6.906036, 0.0005576738,
                   6.901777, 0.0005760576,
                   6.901784, 0.0005760273,
                   6.897395, 0.000595561,
                   6.897402, 0.0005955298,
                   6.892875, 0.0006163207,
                   6.892882, 0.0006162885,
                   6.888209, 0.0006384575,
                   6.888216, 0.0006384243,
                   6.883387, 0.0006621081,
                   6.883394, 0.0006620737,
                   6.878399, 0.000687427,
                   6.878406, 0.0006873915,
                   6.873234, 0.0007145903,
                   6.873241, 0.0007145535,
                   6.86788, 0.0007437992,
                   6.867886, 0.0007437611,
                   6.862322, 0.0007752844,
                   6.862329, 0.0007752448,
                   6.856547, 0.0008093118,
                   6.856554, 0.0008092706,
                   6.850537, 0.0008461893,
                   6.850544, 0.0008461464,
                   6.844274, 0.0008862754,
                   6.844281, 0.0008862307,
                   6.837736, 0.0009299898,
                   6.837743, 0.000929943,
                   6.8309, 0.0009778265,
                   6.830906, 0.0009777776,
                   6.823739, 0.001030371,
                   6.823745, 0.00103032,
                   6.816222, 0.001088322,
                   6.816229, 0.001088268,
                   6.808316, 0.00115252,
                   6.808323, 0.001152464,
                   6.79998, 0.001223984,
                   6.799986, 0.001223924,
                   6.791167, 0.001303958,
                   6.791174, 0.001303894,
                   6.781825, 0.00139398,
                   6.781832, 0.001393913,
                   6.77189, 0.001495972,
                   6.771897, 0.0014959,
                   6.761287, 0.001612358,
                   6.761294, 0.001612281,
                   6.749928, 0.001746242,
                   6.749935, 0.001746159,
                   6.737705, 0.001901654,
                   6.737712, 0.001901565,
                   6.724488, 0.002083916,
                   6.724494, 0.00208382,
                   6.710113, 0.002300186,
                   6.71012, 0.00230008,
                   6.694379, 0.002560295,
                   6.694386, 0.002560178,
                   6.677027, 0.002878082,
                   6.677034, 0.002877953,
                   6.657721, 0.003273586,
                   6.657728, 0.003273442,
                   6.636015, 0.003776796,
                   6.636021, 0.003776632,
                   6.6113, 0.004434395,
                   6.611306, 0.004434206,
                   6.58272, 0.005322598,
                   6.582726, 0.005322376,
                   6.549023, 0.006573303,
                   6.54903, 0.006573036,
                   6.508295, 0.008432167,
                   6.508302, 0.008431835,
                   6.457424, 0.0114021,
                   6.45743, 0.01140167,
                   6.390954, 0.0166508,
                   6.39096, 0.01665021,
                   6.298313, 0.0274051,
                   6.29832, 0.0274042,
                   6.155994, 0.05510944,
                   6.156, 0.05510788,
                   5.901172, 0.1572351,
                   5.901178, 0.1572318,
                   5.334393, 0.6379931,
                   5.334398, 0.6379885,
                   4.481027, 0.465604,
                   4.989382, 0.7977062,
                   4.989387, 0.7977064,
                   5.002551, 0.7978756,
                   5.002556, 0.7978756,
                   4.999995, 0.797886,
                   5, 0.797886,
                   4.999995, 0.797886,
                   5, 0.797886),
                 ncol = 2, byrow = TRUE)

# Grid of parameter values, for plot
grid <- seq(-5, 10, length = 1e3)

# Plot optimiser steps for scenario 1
plot(grid, f(grid), type = "l", xlab = "parameter", ylab = "likelihood",
     main = "Local maximum")
points(steps1[seq(1, nrow(steps1), by = 1), ], type = "o",
       col = "firebrick2", pch = 20, cex = 0.7)
points(steps1[1,1], steps1[1,2], pch = "+",
       col = "seagreen", cex = 2.5)
points(steps1[nrow(steps1),1], steps1[nrow(steps1),2],
       pch = "+", col = "royalblue", cex = 2.5)
legend("topleft", legend = c("Starting value", "Estimate"),
       col = c("seagreen", "royalblue"), pch = "+", bty = "n", pt.cex = 2)

# Plot optimiser steps for scenario 2
plot(grid, f(grid), type="l", xlab = "parameter", ylab = "likelihood",
     main = "Global maximum")
points(steps2[seq(1, nrow(steps2), by = 1), ], type = "o",
       col = "firebrick2", pch = 20, cex = 0.7)
points(steps2[1,1], steps2[1,2], pch = "+",
       col = "seagreen", cex = 2.5)
points(steps2[nrow(steps2),1], steps2[nrow(steps2),2], pch = "+",
       col = "royalblue", cex = 2.5)
legend("topleft", legend = c("Starting value", "Estimate"),
       col = c("seagreen", "royalblue"), pch = "+", bty = "n", pt.cex = 2)

## ----readdat, cache = TRUE, results = 'hide', message = FALSE-----------------
# Load package
library(moveHMM)

# Derive step lengths and turning angles
hmmdata <- prepData(haggis_data, type = "UTM")

# Display first rows of data set
head(hmmdata)

## ----stephist2, out.width='.6\\linewidth', fig.width=6, fig.height=5, fig.align="center"----
# Plot histogram of step lengths
hist(hmmdata$step, xlab = "step length", main = "")

## ----steppar------------------------------------------------------------------
# Starting values for the step length parameters
stepMean0 <- c(1, 5) # initial means (one for each state)
stepSD0 <- c(1, 5) # initial standard deviations (one for each state)
stepPar0 <- c(stepMean0, stepSD0)

## ----zeromass-----------------------------------------------------------------
# Indices of steps of length zero
whichzero <- which(hmmdata$step == 0)

# Proportion of steps of length zero in the data set
length(whichzero)/nrow(hmmdata)

## ----anglehist2, out.width='.6\\linewidth', fig.width=6, fig.height=5, fig.align="center"----
# Plot histogram of turning angles
hist(hmmdata$angle, breaks = seq(-pi, pi, length = 15), xlab = "angle", main = "")

## ----vonmises, echo = FALSE, message = FALSE, out.width='.6\\linewidth', fig.width=6, fig.height=5, fig.align="center"----
# Load package for von Mises distribution
library(CircStats)

# Colours for line plots below
cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00")

# Grid of turning angle values
anglegrid <- seq(-pi, pi, length = 1000)

# Concentration parameters
kappas <- c(0, 1, 2, 5, 10, 20)

# Plot the six density functions
plot(anglegrid, dvm(anglegrid, mu = 0, kappa = kappas[1]), type = "l",
     col = cols[1], xlab = "Turning angle (radians)", ylab = "Density",
     ylim = c(0, 2), xaxt = "n")
for(i in 2:length(kappas)) {
    kappa <- kappas[i]
    points(anglegrid, dvm(anglegrid, mu = 0, kappa = kappa), type = "l", col = cols[i])
}
axis(side = 1, at = c(-pi, -pi/2, 0, pi/2, pi),
     labels = expression(-pi, -pi/2, 0, pi/2, pi))
legend("topleft", legend = kappas, col = cols, lty = 1,
       title = expression(kappa), bty = "n")

## ----anglepar-----------------------------------------------------------------
# Starting values for the turning angle parameters
angleMean0 <- c(pi, 0) # initial means (one for each state)
angleCon0 <- c(1, 10) # initial concentrations (one for each state)
anglePar0 <- c(angleMean0, angleCon0)

## ----fithmm-------------------------------------------------------------------
m <- fitHMM(data = hmmdata, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0)

## ----random-------------------------------------------------------------------
# For reproducibility
set.seed(12345)

# Number of tries with different starting values
niter <- 25

# Save list of fitted models
allm <- list()

for(i in 1:niter) {
    # Step length mean
    stepMean0 <- runif(2,
                       min = c(0.5, 3),
                       max = c(2, 8))

    # Step length standard deviation
    stepSD0 <- runif(2,
                     min = c(0.5, 3),
                     max = c(2, 8))

    # Turning angle mean
    angleMean0 <- c(0, 0)

    # Turning angle concentration
    angleCon0 <- runif(2,
                       min = c(0.5, 5),
                       max = c(2, 15))

    # Fit model
    stepPar0 <- c(stepMean0, stepSD0)
    anglePar0 <- c(angleMean0, angleCon0)
    allm[[i]] <- fitHMM(data = hmmdata, nbStates = 2, stepPar0 = stepPar0,
                        anglePar0 = anglePar0)
}

## ----allnllk------------------------------------------------------------------
# Extract likelihoods of fitted models
allnllk <- unlist(lapply(allm, function(m) m$mod$minimum))
allnllk

## ----mbest--------------------------------------------------------------------
# Index of best fitting model (smallest negative log-likelihood)
whichbest <- which.min(allnllk)

# Best fitting model
mbest <- allm[[whichbest]]
mbest

## ----parallel, eval = FALSE---------------------------------------------------
#  # Package for parallel computations
#  library(parallel)
#  
#  # Create cluster of size ncores
#  ncores <- detectCores() - 1
#  cl <- makeCluster(getOption("cl.cores", ncores))
#  # Export objects needed in parallelised function to cluster
#  clusterExport(cl, list("hmmdata", "fitHMM"))
#  
#  # Number of tries with different starting values
#  niter <- 25
#  
#  # Create list of starting values
#  allPar0 <- lapply(as.list(1:niter), function(x) {
#      # Step length mean
#      stepMean0 <- runif(2,
#                         min = c(0.5, 3),
#                         max = c(2, 8))
#  
#      # Step length standard deviation
#      stepSD0 <- runif(2,
#                       min = c(0.5, 3),
#                       max = c(2, 8))
#  
#      # Turning angle mean
#      angleMean0 <- c(0, 0)
#  
#      # Turning angle concentration
#      angleCon0 <- runif(2,
#                         min = c(0.5, 5),
#                         max = c(2, 15))
#  
#      # Return vectors of starting values
#      stepPar0 <- c(stepMean0, stepSD0)
#      anglePar0 <- c(angleMean0, angleCon0)
#      return(list(step = stepPar0, angle = anglePar0))
#  })
#  
#  # Fit the niter models in parallel
#  allm_parallel <- parLapply(cl = cl, X = allPar0, fun = function(par0) {
#      m <- fitHMM(data = hmmdata, nbStates = 2, stepPar0 = par0$step,
#                  anglePar0 = par0$angle)
#      return(m)
#  })
#  
#  # Then, we can extract the best-fitting model from allm_parallel
#  # as before