vignettes/examples/buffaloExample.R

library(momentuHMM)
library(raster)
library(ctmcmove)

## download buffalo data
load(url("https://github.com/henryrscharf/Hooten_et_al_EL_2018/raw/master/data/buffalo/buffalo_Cilla.RData"))

## download distance to water covariate raster
load(url("https://github.com/henryrscharf/Hooten_et_al_EL_2018/raw/master/data/buffalo/dist2sabie.RData"))
names(dist2sabie) <- "dist2sabie"

## standardize dist2sabie based on slope of gradient
dist2sabie_scaled <- dist2sabie / mean(values(raster::terrain(dist2sabie, opt = "slope")), na.rm = T)

# calculate gradient
D_scaled <- ctmcmove::rast.grad(dist2sabie_scaled)

## W (recharge function covariates)
# near_sabie = indicator for <500m from water
intercept <- raster(dist2sabie)
values(intercept) <- 1
W <- stack(list("intercept" = intercept,
                "near_sabie" = dist2sabie < 0.5e3))
W_names <- names(W)

## orthogonalize W based on locations ----
W_ortho <- W
W_path <- extract(x = W, y = matrix(buffalo_proj@coords, ncol = 2))
obstimes <- as.numeric(buffalo_proj$POSIX) / 3600 # convert time to numeric hours
W_tilde <- apply(W_path * c(0, diff(obstimes)), 2, cumsum)
W_tilde_svd <- svd(W_tilde)
W_tilde_proj_mat <- W_tilde_svd$v %*% diag(W_tilde_svd$d^(-1))
W_mat <- as.matrix(W)
W_mat_proj <- W_mat %*% W_tilde_proj_mat
for(layer in 1:ncol(W_mat)){
  values(W_ortho[[layer]]) <- W_mat_proj[, layer]
  names(W_ortho[[layer]]) <- paste0("svd", layer)
}

lnError <- crawl::argosDiag2Cov(50,50,0) # assume 50m isotropic error ellipse
buffaloData <- data.frame(ID=1,
                          time=obstimes,
                          x=buffalo_proj@coords[, 1],
                          y=buffalo_proj@coords[, 2],
                          ln.sd.x=lnError$ln.sd.x, 
                          ln.sd.y = lnError$ln.sd.y, 
                          error.corr= lnError$error.corr)

# plot observed track and distance to water raster
plotSpatialCov(buffaloData,dist2sabie)

crwOut <- crawlWrap(buffaloData,theta=c(6.5,-.1),fixPar=c(1,1,NA,NA),
                    err.model = list(x=~ln.sd.x-1,y=~ln.sd.y-1,rho=~error.corr),
                    timeStep=0.25, # predict at 15 min time steps
                    attempts=10)

spatialCovs <- list(W_intercept=W_ortho$svd1,
                    W_near_sabie=W_ortho$svd2,
                    dist2sabie=dist2sabie,
                    D.x=D_scaled$rast.grad.x,
                    D.y=D_scaled$rast.grad.y)

hmmData <- prepData(crwOut, spatialCovs = spatialCovs, altCoordNames = "mu")

nbStates <- 2
stateNames <- c("charged","discharged")
dist <- list(mu="rw_mvnorm2") # bivariate normal random walk (mu.x, mu.y)

DM <- list(mu=matrix(c("mu.x_tm1",         0,    0,0,0,0,
                       "mu.x_tm1",         0,"D.x",0,0,0,
                                0,"mu.y_tm1",    0,0,0,0,
                                0,"mu.y_tm1","D.y",0,0,0,
                                0,         0,    0,1,0,0,
                                0,         0,    0,0,1,0,
                                0,         0,    0,0,0,1,
                                0,         0,    0,0,0,1,
                                0,         0,    0,1,0,0,
                                0,         0,    0,0,1,0),
                     5*nbStates,
                     6,byrow=TRUE,
                     dimnames=list(c(paste0("mean.",rep(c("x_","y_"),each=nbStates),1:nbStates),
                                     paste0("sigma.",rep(c("x_","xy_","y_"),each=nbStates),1:nbStates)),
                                   c("x:x_tm1",
                                     "y:y_tm1",
                                     "xy:D",
                                     "sigma_1:(Intercept)",
                                     "sigma_2:(Intercept)",
                                     "sigma_12:(Intercept)"))))

# starting values
Par0=list(mu=c(1, 1, 0, log(85872.66), log(37753.53), 0))
g0 <- 0 # recharge function at time 0
theta <- c(0,0,0) # recharge function parameters

## specify recharge formula
# note that formula for theta requires an 'intercept' term, so we fix it to zero using fixPar
formula <- ~ recharge(g0=~1,theta=~W_intercept+W_near_sabie)

## remove Markov property
betaRef <- c(1,1) # make state 1 the reference state
betaCons <- matrix(c(1,2),2,2) # 1 -> 1 = 2 -> 1 and 1 -> 2 = 2 -> 2

## set fixed parameters
# no intercept or estimated coefficients for t.p.m based entirely on recharge model
# recharge coefficient set to -1 because charged state (state 1) is the reference state
# initial distribution doesn't affect likelihood when Markov property is removed
fixPar <- list(mu = c(Par0$mu[1:2],NA,NA,NA,Par0$mu[6]),
               beta = matrix(c(0,-1,0,-1),2,2),
               delta = c(0.5,0.5),
               theta = c(0,NA,NA)) # fix extra 'intercept' term to zero

# check recharge model specification
checkPar0(hmmData,nbStates=nbStates,dist=dist,formula=formula,
          Par0=Par0,beta0=list(beta=fixPar$beta,g0=g0,theta=theta),delta0=fixPar$delta,
          fixPar=fixPar,DM=DM,betaRef=betaRef,betaCons=betaCons,stateNames=stateNames)

# fit to best predicted path from crawl to get starting values for multiple imputation
buffaloFit <- fitHMM(hmmData,nbStates=nbStates,dist=dist,formula=formula,
                     Par0=Par0,beta0=list(g0=g0,theta=theta),
                     fixPar=fixPar,DM=DM,betaRef=betaRef,betaCons=betaCons,stateNames=stateNames,
                     mvnCoords="mu",optMethod="Nelder-Mead",control=list(maxit=1000))

# multiple imputation fits
bestPar <- getPar(buffaloFit)
set.seed(1,kind="Mersenne-Twister",normal.kind="Inversion")
buffaloFits <- MIfitHMM(crwOut, nSims=28, ncores=4,
                        spatialCovs = spatialCovs, 
                        mvnCoords="mu", altCoordNames = "mu",
                        nbStates=nbStates, dist=dist, formula=formula,
                        Par0=bestPar$Par, beta0=bestPar$beta,
                        fixPar=fixPar, DM=DM, betaRef=betaRef, betaCons=betaCons, stateNames=stateNames,
                        retryFits=3, retrySD=list(mu=c(0,0,3,0,0,0),g0=1,theta=c(0,1,1)),
                        optMethod="Nelder-Mead",control=list(maxit=100000))
buffaloFits
plot(buffaloFits,plotCI=TRUE,ask=FALSE)
plotSpatialCov(buffaloFits,dist2sabie)

trProbs <- getTrProbs(buffaloFits, getCI=TRUE)
# plot estimates and CIs for Pr(discharged) at each time step
plot(trProbs$est[1,2,],type="l", 
     ylim=c(0,1), ylab="Pr(discharged)", xlab="t", col=c("#E69F00", "#56B4E9")[buffaloFits$miSum$Par$states])
arrows(1:dim(trProbs$est)[3],
       trProbs$lower[1,2,],
       1:dim(trProbs$est)[3],
       trProbs$upper[1,2,],
       length=0.025, angle=90, code=3, col=c("#E69F00", "#56B4E9")[buffaloFits$miSum$Par$states], lwd=1.3)
abline(h=0.5,lty=2)

# proportion of entire time series spent in each state
buffaloFits$miSum$Par$timeInStates

# histograms of distance to water by state
par(mfrow=c(2,1))
hist(buffaloFits$miSum$data$dist2sabie[which(buffaloFits$miSum$Par$states==1)],main=stateNames[1],xlab="distance to water (m)")
hist(buffaloFits$miSum$data$dist2sabie[which(buffaloFits$miSum$Par$states==2)],main=stateNames[2],xlab="distance to water (m)")

save.image("buffaloExample.RData")

Try the momentuHMM package in your browser

Any scripts or data that you put into this service are public.

momentuHMM documentation built on Oct. 19, 2022, 1:07 a.m.