update buffalo recharge example in vignette

Author: brettmcclintock-NOAA

DESCRIPTION

@@ -2,7 +2,7 @@ Package: momentuHMM
 Type: Package
 Title: Maximum Likelihood Analysis of Animal Movement Behavior Using Multivariate Hidden Markov Models
 Version: 2.0.0
-Date: 2024-02-02
+Date: 2024-06-12
 Depends:
     R (>= 2.10)
 Author: Brett McClintock, Theo Michelot

NEWS

@@ -1,4 +1,4 @@
-momentuHMM 2.0.0 2024-02-02
+momentuHMM 2.0.0 2024-06-12
 -----------------
 NEW FEATURES
 

publications/momentuHMM.bib

@@ -1,6 +1,5 @@
 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"))

publications/publications.md

@@ -2783,7 +2783,7 @@ with recharge function
 where $w_j$ is the distance to nearest water indicator covariate at location ${\boldsymbol \mu}_j$. Thus the probability of being in the ``discharged'' state decreases as the recharge function $(g_t)$ increases. Note that there are no t.p.m. working parameter coefficients in Eq. \ref{eq:rechargeTPM1} and, because $\gamma_{11}=\gamma_{21}$ and $\gamma_{12}=\gamma_{22}$, the state switching in this model is not Markov (i.e., the state at time $t$ does not depend on the state at time $t-1$). 
 
 In order to fit this model, we must first load the data and format the covariate rasters:
-<<recharge-1, echo=TRUE, eval=TRUE, cache=TRUE, message=FALSE>>=
+<<recharge-1, echo=TRUE, eval=TRUE, cache=TRUE, message=FALSE, warning=FALSE>>=
 library(raster)
 
 ## download buffalo data
@@ -2802,9 +2802,6 @@ dist2sabie_scaled <- dist2sabie / mean(values(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)
@@ -2858,16 +2855,16 @@ crwOut <- crawlWrap(buffaloData,
 @
 
 Now we're ready to specify the recharge model. We'll first fit the model to the best predicted track from \verb|crawlWrap| and then use this model fit to specify starting values for a multiple imputation analysis:
-<<recharge-3, echo=TRUE, eval=TRUE, cache=TRUE, message=FALSE>>=
+<<recharge-3, echo=TRUE, eval=TRUE, cache=TRUE, message=FALSE, warning=FALSE>>=
 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)
+                    D = dist2sabie_scaled)
 
 # best predicted track data
 hmmData <- prepData(crwOut, 
                     spatialCovs = spatialCovs, 
+                    gradient = TRUE,
                     altCoordNames = "mu")
 head(hmmData[,c("ID","time", "mu.x", "mu.y",
                 "W_intercept","W_near_sabie","dist2sabie",
@@ -2968,7 +2965,7 @@ Now that we have our starting values (``bestPar''), let's fit $28$ imputations o
 <<recharge-4, echo=-1, eval=FALSE>>=
 set.seed(1,kind="Mersenne-Twister",normal.kind="Inversion")
 buffaloFits <- MIfitHMM(crwOut, nSims=28,
-                        spatialCovs = spatialCovs, 
+                        spatialCovs = spatialCovs, gradient = TRUE,
                         mvnCoords="mu", altCoordNames = "mu",
                         nbStates=nbStates, dist=dist, formula=formula,
                         Par0=bestPar$Par, beta0=bestPar$beta,
@@ -2987,7 +2984,7 @@ plotSpatialCov(buffaloFits,dist2sabie)
 # plot estimates and CIs for Pr(discharged) at each time step
 trProbs <- getTrProbs(buffaloFits, getCI=TRUE)
 plot(trProbs$est[1,2,],type="l", ylim=c(0,1), 
-     ylab="Pr(discharged)", xlab="t", 
+     ylab="Pr(discharged)", xlab="time step", 
      col=c("#E69F00", "#56B4E9")[buffaloFits$miSum$Par$states])
 arrows(1:dim(trProbs$est)[3],
        trProbs$lower[1,2,],

vignettes/codDist.pdf

@@ -4359,9 +4359,6 @@ \subsection{African buffalo recharge dynamics}
                                                       \hlkwc{opt} \hlstd{=} \hlstr{"slope"}\hlstd{)),}
                                        \hlkwc{na.rm} \hlstd{= T)}
 
-\hlcom{# calculate gradient}
-\hlstd{D_scaled} \hlkwb{<-} \hlstd{ctmcmove}\hlopt{::}\hlkwd{rast.grad}\hlstd{(dist2sabie_scaled)}
-
 \hlcom{## W (recharge function covariates)}
 \hlcom{# near_sabie = indicator for <500m from water}
 \hlstd{intercept} \hlkwb{<-} \hlkwd{raster}\hlstd{(dist2sabie)}
@@ -4427,32 +4424,32 @@ \subsection{African buffalo recharge dynamics}
 \hlstd{spatialCovs} \hlkwb{<-} \hlkwd{list}\hlstd{(}\hlkwc{W_intercept} \hlstd{= W_ortho}\hlopt{$}\hlstd{svd1,}
                     \hlkwc{W_near_sabie} \hlstd{= W_ortho}\hlopt{$}\hlstd{svd2,}
                     \hlkwc{dist2sabie} \hlstd{= dist2sabie,}
-                    \hlkwc{D.x} \hlstd{= D_scaled}\hlopt{$}\hlstd{rast.grad.x,}
-                    \hlkwc{D.y} \hlstd{= D_scaled}\hlopt{$}\hlstd{rast.grad.y)}
+                    \hlkwc{D} \hlstd{= dist2sabie_scaled)}
 
 \hlcom{# best predicted track data}
 \hlstd{hmmData} \hlkwb{<-} \hlkwd{prepData}\hlstd{(crwOut,}
                     \hlkwc{spatialCovs} \hlstd{= spatialCovs,}
+                    \hlkwc{gradient} \hlstd{=} \hlnum{TRUE}\hlstd{,}
                     \hlkwc{altCoordNames} \hlstd{=} \hlstr{"mu"}\hlstd{)}
 \hlkwd{head}\hlstd{(hmmData[,}\hlkwd{c}\hlstd{(}\hlstr{"ID"}\hlstd{,}\hlstr{"time"}\hlstd{,} \hlstr{"mu.x"}\hlstd{,} \hlstr{"mu.y"}\hlstd{,}
                 \hlstr{"W_intercept"}\hlstd{,}\hlstr{"W_near_sabie"}\hlstd{,}\hlstr{"dist2sabie"}\hlstd{,}
                 \hlstr{"D.x"}\hlstd{,}\hlstr{"D.y"}\hlstd{)])}
 \end{alltt}
 \begin{verbatim}
-##   ID     time     mu.x     mu.y   W_intercept W_near_sabie dist2sabie
-## 1  1 313344.4 384560.8 -2770752 -0.0001873081 -0.003616337   1323.424
-## 2  1 313344.7 384560.3 -2770751 -0.0001873081 -0.003616337   1323.424
-## 3  1 313344.9 384559.8 -2770749 -0.0001873081 -0.003616337   1323.424
-## 4  1 313345.2 384559.3 -2770747 -0.0001873081 -0.003616337   1323.424
-## 5  1 313345.4 384559.1 -2770744 -0.0001873081 -0.003616337   1323.424
-## 6  1 313345.7 384559.1 -2770742 -0.0001873081 -0.003616337   1323.424
-##          D.x       D.y
-## 1 -0.8709203 0.2576337
-## 2 -0.8709203 0.2576337
-## 3 -0.8709203 0.2576337
-## 4 -0.8709203 0.2576337
-## 5 -0.8709203 0.2576337
-## 6 -0.8709203 0.2576337
+##   ID     time     mu.x     mu.y   W_intercept W_near_sabie dist2sabie       D.x
+## 1  1 313344.4 384560.8 -2770752 -0.0001873081 -0.003616337   1323.424 -1.228246
+## 2  1 313344.7 384560.3 -2770751 -0.0001873081 -0.003616337   1323.424 -1.227881
+## 3  1 313344.9 384559.8 -2770749 -0.0001873081 -0.003616337   1323.424 -1.227418
+## 4  1 313345.2 384559.3 -2770747 -0.0001873081 -0.003616337   1323.424 -1.226851
+## 5  1 313345.4 384559.1 -2770744 -0.0001873081 -0.003616337   1323.424 -1.226157
+## 6  1 313345.7 384559.1 -2770742 -0.0001873081 -0.003616337   1323.424 -1.225450
+##         D.y
+## 1 0.3505877
+## 2 0.3779253
+## 3 0.3777812
+## 4 0.3776649
+## 5 0.3776060
+## 6 0.3775941
 \end{verbatim}
 \begin{alltt}
 \hlstd{nbStates} \hlkwb{<-} \hlnum{2}
@@ -4582,7 +4579,7 @@ \subsection{African buffalo recharge dynamics}
 \definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
 \begin{alltt}
 \hlstd{buffaloFits} \hlkwb{<-} \hlkwd{MIfitHMM}\hlstd{(crwOut,} \hlkwc{nSims}\hlstd{=}\hlnum{28}\hlstd{,}
-                        \hlkwc{spatialCovs} \hlstd{= spatialCovs,}
+                        \hlkwc{spatialCovs} \hlstd{= spatialCovs,} \hlkwc{gradient} \hlstd{=} \hlnum{TRUE}\hlstd{,}
                         \hlkwc{mvnCoords}\hlstd{=}\hlstr{"mu"}\hlstd{,} \hlkwc{altCoordNames} \hlstd{=} \hlstr{"mu"}\hlstd{,}
                         \hlkwc{nbStates}\hlstd{=nbStates,} \hlkwc{dist}\hlstd{=dist,} \hlkwc{formula}\hlstd{=formula,}
                         \hlkwc{Par0}\hlstd{=bestPar}\hlopt{$}\hlstd{Par,} \hlkwc{beta0}\hlstd{=bestPar}\hlopt{$}\hlstd{beta,}
@@ -4601,7 +4598,7 @@ \subsection{African buffalo recharge dynamics}
 \hlcom{# plot estimates and CIs for Pr(discharged) at each time step}
 \hlstd{trProbs} \hlkwb{<-} \hlkwd{getTrProbs}\hlstd{(buffaloFits,} \hlkwc{getCI}\hlstd{=}\hlnum{TRUE}\hlstd{)}
 \hlkwd{plot}\hlstd{(trProbs}\hlopt{$}\hlstd{est[}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{,],}\hlkwc{type}\hlstd{=}\hlstr{"l"}\hlstd{,} \hlkwc{ylim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{1}\hlstd{),}
-     \hlkwc{ylab}\hlstd{=}\hlstr{"Pr(discharged)"}\hlstd{,} \hlkwc{xlab}\hlstd{=}\hlstr{"t"}\hlstd{,}
+     \hlkwc{ylab}\hlstd{=}\hlstr{"Pr(discharged)"}\hlstd{,} \hlkwc{xlab}\hlstd{=}\hlstr{"time step"}\hlstd{,}
      \hlkwc{col}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"#E69F00"}\hlstd{,} \hlstr{"#56B4E9"}\hlstd{)[buffaloFits}\hlopt{$}\hlstd{miSum}\hlopt{$}\hlstd{Par}\hlopt{$}\hlstd{states])}
 \hlkwd{arrows}\hlstd{(}\hlnum{1}\hlopt{:}\hlkwd{dim}\hlstd{(trProbs}\hlopt{$}\hlstd{est)[}\hlnum{3}\hlstd{],}
        \hlstd{trProbs}\hlopt{$}\hlstd{lower[}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{,],}
@@ -4614,8 +4611,8 @@ \subsection{African buffalo recharge dynamics}
 \end{alltt}
 \end{kframe}
 \end{knitrout}
-\noindent As in \cite{HootenEtAl2019}, we found that the buffalo spent a majority of time steps in the discharged state ($73$\%, 95\% CI: $70-76\%$) and thus needed to recharge regularly near water resources. With estimated $g_0=-0.27$ (95\% CI: $-0.75-0.2$), $\theta_1=1.43$ (95\% CI: $-0.61-3.48$), and $\theta_2=2.52$ (95\% CI: $1.47-3.57$), the estimated recharge function and transition probabilities (Figure \ref{fig:recharge}) look very similar to those reported by \cite{HootenEtAl2019}. However, \cite{HootenEtAl2019} found some evidence that the buffalo orients
-toward surface water when in the discharged state, but our discrete-time formulation did not find evidence of such biased movement ($\beta^\mu=0.88$, 95\% CI: $-3.2-4.96$). This difference could be attributable to several factors, including our formulation being in discrete time (instead of continuous time), our use of a 2-stage multiple imputation approach based on the CTCRW (instead of a single-stage model), and the absence of prior distributions in our non-Bayesian model. Nevertheless, inferences about recharge and state-switching dynamics are essentially the same between our discrete-time formulation and the continuous-time model of \cite{HootenEtAl2019}.
+\noindent As in \cite{HootenEtAl2019}, we found that the buffalo spent a majority of time steps in the discharged state ($73$\%, 95\% CI: $70-76\%$) and thus needed to recharge regularly near water resources. With estimated $g_0=-0.27$ (95\% CI: $-0.75-0.2$), $\theta_1=1.43$ (95\% CI: $-0.62-3.47$), and $\theta_2=2.52$ (95\% CI: $1.46-3.57$), the estimated recharge function and transition probabilities (Figure \ref{fig:recharge}) look very similar to those reported by \cite{HootenEtAl2019}. However, \cite{HootenEtAl2019} found some evidence that the buffalo orients
+toward surface water when in the discharged state, but our discrete-time formulation did not find evidence of such biased movement ($\beta^\mu=0.64$, 95\% CI: $-2.16-3.44$). This difference could be attributable to several factors, including our formulation being in discrete time (instead of continuous time), our use of a 2-stage multiple imputation approach based on the CTCRW (instead of a single-stage model), and the absence of prior distributions in our non-Bayesian model. Nevertheless, inferences about recharge and state-switching dynamics are essentially the same between our discrete-time formulation and the continuous-time model of \cite{HootenEtAl2019}.
 \begin{figure}[htbp]
   \centering
   \includegraphics[width=0.9\textwidth]{plot_buffaloStates.pdf}\\

vignettes/codStates.pdf

@@ -35,9 +35,9 @@ dev.off()
 for(plt in seq(1,17)[-c(1,2,9,10,12,13,15,16,17)])
   unlink(paste0("plot_elephantResults0",ifelse(plt>9,"","0"),plt,".pdf"))
 
-png(filename="elephant_plotSat.png",width=5,height=5,units="in",res=80)
-plotSat(data.frame(x=rawData$lon,y=rawData$lat),zoom=8,location=c(median(rawData$lon),median(rawData$lat)),ask=FALSE)
-dev.off()
+#png(filename="elephant_plotSat.png",width=5,height=5,units="in",res=80)
+#plotSat(data.frame(x=rawData$lon,y=rawData$lat),zoom=8,location=c(median(rawData$lon),median(rawData$lat)),ask=FALSE)
+#dev.off()
 
 rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])
 
@@ -240,15 +240,15 @@ tmpData<-spTransform(tmpData,CRS="+proj=longlat +ellps=WGS84")
 png(file=paste0(getwd(),"/plot_harbourSeal.png"),width=7.25,height=5,units="in",res=84)
 par(mfrow=c(1,2))
 for(id in c(8,1)){
-  freqs<-miSum.ind$Par$states[hsData$ID==id]
-  x<-tmpData$x[which(hsData$ID==id)]
-  y<-tmpData$y[which(hsData$ID==id)]
+  freqs<-miSum.ind$Par$states[which(hsData$ID==id)]
+  x<-coordinates(tmpData)[which(hsData$ID==id),1]
+  y<-coordinates(tmpData)[which(hsData$ID==id),2]
   errorEllipse<-miSum.ind$errorEllipse[which(hsData$ID==id)]
   errorEllipse<-lapply(errorEllipse,function(x){x<-as.data.frame(x);
-                                                coordinates(x)<-c("x","y");
-                                                proj4string(x)<-CRS("+init=epsg:27700 +units=m");
-                                                x<-spTransform(x,CRS="+proj=longlat +ellps=WGS84");
-                                                as.data.frame(x)})
+  coordinates(x)<-c("x","y");
+  proj4string(x)<-CRS("+init=epsg:27700 +units=m");
+  x<-spTransform(x,CRS="+proj=longlat +ellps=WGS84");
+  as.data.frame(x)})
 
   #png(file=paste0(getwd(),"/plot_harbourSeal",id,".png"),width=5,height=7.25,units="in",res=90)
   plot(x,y,type="o",pch=20,cex=.5,xlab="longitude",ylab="latitude",cex.lab=0.8,cex.axis=.7)
@@ -287,9 +287,9 @@ append.RData(timeIn1,file=paste0(getwd(),"/vignette_inputs.RData"))
 append.RData(timeIn2,file=paste0(getwd(),"/vignette_inputs.RData"))
 #append.RData(m2,file=paste0(getwd(),"/vignette_inputs.RData"))
 #append.RData(newProj,file=paste0(getwd(),"/vignette_inputs.RData"))
-png(file=paste0(getwd(),"/plot_northernFulmarExample.png"),width=7.25,height=5,units="in",res=80)
-plotSat(m2,zoom=7,shape=c(17,1,17,1,17,1),size=2,col=rep(c("#E69F00", "#56B4E9", "#009E73"),each=2),stateNames=c("sea ARS","sea Transit","boat ARS","boat Transit","colony ARS","colony Transit"),projargs=newProj,ask=FALSE)
-dev.off()
+#png(file=paste0(getwd(),"/plot_northernFulmarExample.png"),width=7.25,height=5,units="in",res=80)
+#plotSat(m2,zoom=7,shape=c(17,1,17,1,17,1),size=2,col=rep(c("#E69F00", "#56B4E9", "#009E73"),each=2),stateNames=c("sea ARS","sea Transit","boat ARS","boat Transit","colony ARS","colony Transit"),projargs=newProj,ask=FALSE)
+#dev.off()
 
 rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])
 
@@ -326,11 +326,11 @@ fitmix1_Par <- getPar(fitmix1)
 append.RData(fitmix1_Par,file=paste0(getwd(),"/vignette_inputs.RData"))
 Par0_mix2 <- getPar0(fitmix1,mixtures=2)
 Par0_mix2$beta$beta[1,] <- c(-2.26, -3.93, -0.58, 
-                              0.03, -2.25, -0.26, 
+                             0.03, -2.25, -0.26, 
                              -3.38, -4.79, -2.82, 
                              -1.06,  -3.3, -3.43)
 Par0_mix2$beta$beta[2,] <- c(-2.51, -3.32, -2.63, 
-                              0.03, -1.26, -0.12, 
+                             0.03, -1.26, -0.12, 
                              -96.8, -3.62, -1.75, 
                              -1.76, -2.14, -1.38)
 Par0_mix2$beta$pi <- c(0.73, 0.27)
@@ -409,7 +409,7 @@ rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])
 ### buffalo recharge example
 ###################################################
 #source(paste0(getwd(),"/examples/buffaloExample.R"))
-load(paste0(example_wd,"buffaloExample.RData"))
+load(paste0(example_wd,"buffaloExample_gradCT.RData"))
 bufPar <- buffaloFits$miSum$Par$beta[c("mu","g0","theta")]
 bufPar$timeInStates <- buffaloFits$miSum$Par$timeInStates
 append.RData(bufPar,file=paste0(getwd(),"/vignette_inputs.RData"))
@@ -430,14 +430,38 @@ for(plt in c(1:10,12))
 pdf(file=paste0(getwd(),"/plot_buffaloResults.pdf"),width=8,height=3)
 par(mar=c(5,4,4,2)-c(0,0,2,1)) # bottom, left, top, right
 plot(trProbs$est[1,2,],type="l", 
-     ylim=c(0,1), ylab="Pr(discharged)", xlab="t", col=c("#E69F00", "#56B4E9")[buffaloFits$miSum$Par$states])
+     ylim=c(0,1), ylab="Pr(discharged)", xlab="time step", 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)
 dev.off()
+
+#pdf(file=paste0(getwd(),"/plot_buffaloExampleCT%03d.pdf"),width=8,height=3,onefile = FALSE)
+#plot(buffaloFitsCT,plotCI=TRUE,legend.pos="bottom",ask=FALSE)
+#dev.off()
+
+#pdf(file="plot_buffaloStatesCT.pdf",width=7.5,height=5)
+#plotSpatialCov(buffaloFitsCT,dist2sabie)
+#dev.off()
+
+#for(plt in c(1:10,12))
+#  unlink(paste0("plot_buffaloExampleCT0",ifelse(plt>9,"","0"),plt,".pdf"))
+
+#pdf(file=paste0(getwd(),"/plot_buffaloResultsCT.pdf"),width=8,height=3)
+#par(mar=c(5,4,4,2)-c(0,0,2,1)) # bottom, left, top, right
+#plot(cumsum(diff(buffaloFitsCT$miSum$data$time)),trProbsCT$est[1,2,-1],type="l", 
+#     ylim=c(0,1), ylab="Pr(discharged)", xlab="t", col=c("#E69F00", "#56B4E9")[buffaloFitsCT$miSum$Par$states])
+#arrows(cumsum(diff(buffaloFitsCT$miSum$data$time)),
+#       trProbsCT$lower[1,2,-1],
+#       cumsum(diff(buffaloFitsCT$miSum$data$time)),
+#       trProbsCT$upper[1,2,-1],
+#       length=0.025, angle=90, code=3, col=c("#E69F00", "#56B4E9")[buffaloFitsCT$miSum$Par$states], lwd=1.3)
+#abline(h=0.5,lty=2)
+#dev.off()
+
 rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])
 
 ###################################################
@@ -455,8 +479,10 @@ rm(list=ls()[-which(ls()=="example_wd")])
 ### Langevin example
 ###################################################
 load(paste0(example_wd,"langevinExample.RData"))
-library(patchwork)
+library(raster)
+library(viridis)
 library(ggplot2)
+library(patchwork)
 covnames <- c("bathy", "slope", "d2site")
 covplot <- list()
 for(i in 1:3) {