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

## ----init, echo=FALSE, message=FALSE------------------------------------------
library(MASS)
library(boot)
library(CircStats)
library(sp)
library(Rcpp)
library(moveHMM)

## ----trackData2,size='small'--------------------------------------------------
head(elk_data)

## ----trackData4,size='small'--------------------------------------------------
elk_data$Easting <- elk_data$Easting/1000
elk_data$Northing <- elk_data$Northing/1000

## ----trackData5,size='small'--------------------------------------------------
head(elk_data)

## ----prepData,size='small'----------------------------------------------------
data <- prepData(elk_data,type="UTM",coordNames=c("Easting","Northing"))

## ----data,size='small'--------------------------------------------------------
head(data)

## ----summary,size='small'-----------------------------------------------------
summary(data)

## ----plot_moveData1,size='small',eval=FALSE-----------------------------------
#  plot(data,compact=T)

## ----plot_moveData2,echo=FALSE,results='hide'---------------------------------
pdf(file="plot_moveData.pdf")
plot(data,compact=T,ask=FALSE)
dev.off()

## ----fitHMM,size='small',eval=FALSE-------------------------------------------
#  ## standardize covariate values
#  data$dist_water <-
#      (data$dist_water-mean(data$dist_water))/sd(data$dist_water)
#  
#  ## initial parameters for gamma and von Mises distributions
#  mu0 <- c(0.1,1) # step mean (two parameters: one for each state)
#  sigma0 <- c(0.1,1) # step SD
#  zeromass0 <- c(0.1,0.05) # step zero-mass
#  stepPar0 <- c(mu0,sigma0,zeromass0)
#  angleMean0 <- c(pi,0) # angle mean
#  kappa0 <- c(1,1) # angle concentration
#  anglePar0 <- c(angleMean0,kappa0)
#  
#  ## call to fitting function
#  m <- fitHMM(data=data,nbStates=2,stepPar0=stepPar0,
#              anglePar0=anglePar0,formula=~dist_water)

## ----savemodels,echo=FALSE,eval=FALSE-----------------------------------------
#  ## code to fit and save the 2 and 3-state models loaded below
#  library(moveHMM)
#  trackData <- read.table("http://www.esapubs.org/archive/ecol/E085/072/elk_data.txt",
#                          sep="\t",header=TRUE)[(1:735),c(1,2,3,7)]
#  colnames(trackData)[1] <- "ID"
#  colnames(trackData)[4] <- "dist_water"
#  trackData$Easting <- trackData$Easting/1000
#  trackData$Northing <- trackData$Northing/1000
#  
#  data <- prepData(trackData,type="UTM",coordNames=c("Easting","Northing"))
#  data$dist_water <- (data$dist_water-mean(data$dist_water))/sd(data$dist_water)
#  
#  ## Fit a 2-state HMM
#  mu0 <- c(0.1,1)
#  sigma0 <- c(0.1,1)
#  zeromass0 <- c(0.05,0.01)
#  stepPar0 <- c(mu0,sigma0,zeromass0)
#  angleMean0 <- c(pi,0)
#  kappa0 <- c(1,1)
#  anglePar0 <- c(angleMean0,kappa0)
#  
#  m <- fitHMM(data=data,nbStates=2,stepPar0=stepPar0,anglePar0=anglePar0,formula=~dist_water)
#  
#  ## Fit a 3-state HMM
#  mu0 <- c(0.1,0.5,3)
#  sigma0 <- c(0.05,0.5,1)
#  zeromass0 <- c(0.05,0.01,0.01)
#  stepPar0 <- c(mu0,sigma0,zeromass0)
#  angleMean0 <- c(pi,pi,0)
#  kappa0 <- c(1,1,1)
#  anglePar0 <- c(angleMean0,kappa0)
#  
#  m3 <- fitHMM(data=data,nbStates=3,stepPar0=stepPar0,anglePar0=anglePar0,
#               formula=~dist_water)
#  
#  save(m,m3,file="models.RData")

## ----fitHMM2,echo=FALSE,cache=TRUE--------------------------------------------
load("models.RData")

## ----print_moveHMM,size='small'-----------------------------------------------
m

## ----nlm-error, eval = FALSE, highlight = FALSE-------------------------------
#    Error in nlm(nLogLike, wpar, nbStates, bounds, parSize, data, stepDist,  :
#    non-finite value supplied by 'nlm'

## ----normalize,size='small',eval=FALSE----------------------------------------
#  data$dist_water <-
#      (data$dist_water-mean(data$dist_water))/sd(data$dist_water)

## ----conf,size='small',warning=FALSE------------------------------------------
CI(m)

## ----ci-warn, eval = FALSE, highlight = FALSE---------------------------------
#    Warning message:
#    In CI(m) :
#    Some of the parameter estimates seem to lie close to
#    the boundaries of their parameter space. The associated
#    CIs are probably unreliable (or might not be computable).

## ----plot_moveHMM,size='small',eval=FALSE-------------------------------------
#  plot(m, plotCI=TRUE)

## ----plot_moveHMM2,echo=FALSE,results='hide'----------------------------------
pdf(file="plot_moveHMM.pdf")
plot(m,ask=FALSE,plotCI=TRUE)
dev.off()

## ----states,size='small',cache=TRUE-------------------------------------------
states <- viterbi(m)
states[1:25]

## ----stateProbs,size='small',cache=TRUE---------------------------------------
sp <- stateProbs(m)
head(sp)

## ----plotStates,size='small',eval=FALSE---------------------------------------
#  plotStates(m,animals="elk-115")

## ----plotStates2, fig.width='4in', fig.height='4in', out.width='4in', out.height='4in', fig.pos='ht', fig.align='center', fig.cap="Decoded states sequence (top row), and state probabilities of observations (middle and bottom rows) for elk-115", cache=TRUE, echo=FALSE, results='hide'----
plotStates(m,animals="elk-115",ask=FALSE)

## ----plotStationary, fig.width='4in', fig.height='4in', out.width='4in', out.height='4in', fig.pos='htbp', fig.align='center', fig.cap='Output of \\texttt{plotStationary}. Stationary state probabilities, as functions of the distance to water, with 95\\% confidence intervals.'----
plotStationary(m, plotCI=TRUE)

## ----AIC,size='small',eval=FALSE----------------------------------------------
#  # initial parameters
#  mu0 <- c(0.1,0.5,3)
#  sigma0 <- c(0.05,0.5,1)
#  zeromass0 <- c(0.05,0.0001,0.0001)
#  stepPar0 <- c(mu0,sigma0,zeromass0)
#  angleMean0 <- c(pi,pi,0)
#  kappa0 <- c(1,1,1)
#  anglePar0 <- c(angleMean0,kappa0)
#  
#  # fit the 3-state model
#  m3 <- fitHMM(data=data,nbStates=3,stepPar0=stepPar0,
#               anglePar0=anglePar0,formula=~dist_water)

## ----AIC2,size='small'--------------------------------------------------------
AIC(m,m3)

## ----pseudoRes,size='small',eval=FALSE----------------------------------------
#  # compute the pseudo-residuals
#  pr <- pseudoRes(m)
#  
#  # time series, qq-plots, and ACF of the pseudo-residuals
#  plotPR(m)

## ----plotPR, fig.width='4in', fig.height='5in', out.width='4in', out.height='5in', fig.pos='ht!', fig.align='center', fig.cap="Time series, qq-plots, and autocorrelation functions of the pseudo-residuals of the 2-state model.", cache=TRUE, echo=FALSE, results='hide', message=FALSE----
plotPR(m)