Merge branch 'master' into staging

# Conflicts:
#	DESCRIPTION

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: aniMotum
 Title: Fit Continuous-Time State-Space and Latent Variable Models for Quality Control of Argos Satellite (and Other) Telemetry Data and for Estimating Changes in Animal Movement
-Version: 1.2-1.9001
-Date: 2023-11-16
+Version: 1.2-03.9001
+Date: 2024-02-06
 Authors@R: 
     c(
       person(given = "Ian",
@@ -33,7 +33,7 @@ RoxygenNote: 7.2.3
 Imports: 
     tibble (>= 2.1.3),
     ggplot2 (>= 3.4.0),
-    TMB (>= 1.9.2),
+    TMB (>= 1.9.6),
     sf (>= 1.0-14),
     tidyr,
     dplyr (>= 1.0.0),

R/plot.sim_fit.R

@@ -81,9 +81,16 @@ plot.sim_fit <- function(x,
       st_make_valid()
   }
 
-  if(!zoom) bounds <- st_bbox(wm.sf)
+  if(!zoom & ortho) {
+    message("Cannot render globe with ortho projection, setting zoom = TRUE...")
+    #bounds <- st_bbox(wm.sf)
+    bounds <- st_bbox(pos.sf)
+    } else if(!zoom & !ortho) {
+      bounds <- st_bbox(wm.sf)
+      }
+
   if(zoom) bounds <- st_bbox(pos.sf)
-
+  
   ## do plots
   p <- lapply(x$sims, function(x) {
 
@@ -112,15 +119,8 @@ plot.sim_fit <- function(x,
         size = 1
       ) + 
       xlab(element_blank()) +
-      ylab(element_blank())
-
-    if(zoom) {
-      m <- m + 
-      theme_minimal()
-    } else {
-      m <- m +
-        theme_void()
-    }
+      ylab(element_blank()) +
+      theme_void()
 
   })
   ## arrange plots

R/sim_fit.R

@@ -197,16 +197,21 @@ sim_fit <-
                                    sigma = Sigma * dt[i],
                                    lower = vmin,
                                    upper = vmax)
+
                  ## wrap x values, reflect y values
                  mu1 <- wrap_x(mu[i-1,] + v[i,] * dt[i], ex[1:2])
+
                  mu1 <- reflect_y(mu1, ex[3:4])
                  if(!is.null(grad)) {
                    pv <- as.numeric(c(extract(grad[[1]], rbind(mu1))[1],
                            extract(grad[[2]], rbind(mu1))[1]))
+
+                   if(all(is.na(pv))) pv <- c(0,0) # unsure why NA's can creep in here
                    mu[i,] <- mu1 + pv * beta
                  } else {
                    mu[i,] <- mu1
                  }
+
                }
                df <- data.frame(
                  rep = j,

data-raw/potential_fn_gradients.R

@@ -1,6 +1,6 @@
 ## code to prepare `potential_fn_gradients` dataset goes here
 
-usethis::use_data(potential_fn_gradients, overwrite = TRUE)
+#usethis::use_data(potential_fn_gradients, overwrite = TRUE)
 
 ## IDJ - 
 
@@ -10,9 +10,11 @@ proj <- "+proj=merc +units=km +datum=WGS84"
 ## create world land raster
 wm <- rnaturalearth::ne_countries(scale = 10, returnclass = "sf")
 
+## set resolution to ensure a sufficiently small file size for the R pkg
+##  finer resolutions may be more useful for custom applications
 y <- terra::rast(crs = sf::st_crs(wm)$wkt,
           vals = 1,
-          resolution = c(0.5, 0.5),
+          resolution = c(0.25, 0.25),  
           xmin = -180,
           xmax = 180,
           ymin = -86,
@@ -27,13 +29,15 @@ x[is.na(x)] <- -1
 x[x == 1] <- NA
 x[x == -1] <- 1
 
-## calculate gradient rasters - these are needed to keep fish off land
+## calculate gradient rasters - these are needed to keep animals off land
 dist <- terra::distance(x)
 x1 <- terra::terrain(dist, v = "slope", unit = "radians")
 y1 <- terra::terrain(dist, v = "aspect", unit = "radians")
 grad.x <- -1 * x1 * cos(0.5 * pi - y1)
 grad.y <- -1 * x1 * sin(0.5 * pi - y1)
 grad <- c(grad.x, grad.y)
 
+## need to use terra::wrap() to export the rasters & terra::unwrap() to import 
+##     from an exported (wrapped) file
 grad <- terra::wrap(grad)
 save(grad, file = "inst/extdata/grad.rda", compress = "xz")

inst/TMB-version

@@ -1 +1 @@
-1.9.6
+1.9.10

vignettes/Mapping.Rmd

@@ -29,14 +29,6 @@ fit <- fit_ssm(se2,
 
 map(fit, what = "predicted")
 ```
-
-```{r map 1a, echo=FALSE, fig.height=4, fig.width=7, warning=FALSE, message=FALSE}
-se2 <- subset(sese, id %in% unique(id)[1:2])
-fit <- fit_ssm(se2, 
-               model = "mp", 
-               time.step = 24, 
-               control = ssm_control(verbose = 0))
-```
 ![](images/mapping/map_1a.jpg){width=95%}
 
 Here, SSM-predicted locations are coloured by the estimated move persistence index ($\gamma_t$). To aid visualisation of high and low move persistence regions along each track, the $\gamma_t$ estimates are normalised separately for each track to span the interval 0,1. This can be turned off with `normalise = FALSE`, or applied to the group of tracks with `group = TRUE`.

vignettes/Move_persistence_models.Rmd

@@ -16,7 +16,7 @@ require(aniMotum)
 require(ggplot2)
 ```
 `aniMotum` provides a few different ways to fit a move persistence model to obtain estimates of a continuous-valued behavioural index along individual tracks (Auger-Méthé et al. 2017; Jonsen et al. 2019). In `aniMotum` versions prior to 1.0-5, move persistence was estimated via `fit_mpm()` and, when fitting to Argos or other location error-prone data types, required a two-stage approach of first calling `fit_ssm()` and then calling `fit_mpm()`. The first step fits a state-space model to account for location error in the data:
-```{r fit ssm, warning=FALSE, message=FALSE}
+```{r fit ssm, eval=FALSE, warning=FALSE, message=FALSE}
 d <- subset(sese, id == "ct109-186-14")
 fit <- fit_ssm(d, 
                vmax = 3,
@@ -26,7 +26,7 @@ fit <- fit_ssm(d,
 ```
 
 The second step estimates move persistence from the SSM-estimated locations:
-```{r fit mpm, warning=FALSE, message=FALSE}
+```{r fit mpm, eval=FALSE, warning=FALSE, message=FALSE}
 fmp <- fit_mpm(fit, 
                what = "predicted", 
                model = "mpm",
@@ -48,7 +48,7 @@ map(fit, fmp, what = "predicted", silent = TRUE)
 `aniMotum` now provides a move persistence model in state-space form to allow simultaneous estimation of move persistence and location states, and this is invoked via `fit_ssm()` by specifying `model = "mp"`. Currently, this model can only be fit to individuals separately.
 
 Here, we fit the move persistence SSM to the same southern elephant seal track as above:
-```{r fit_ssm, warning=FALSE, message=FALSE}
+```{r fit_ssm, eval=FALSE, warning=FALSE, message=FALSE}
 fit <- fit_ssm(d, 
                vmax = 3, 
                model = "mp", 

vignettes/Overview.Rmd

@@ -186,7 +186,7 @@ Original variable names are not preserved in the output object `fit` but rather
 
 `fit_ssm` can be applied to single or multiple tracks, without modification. The specified SSM is fit to each individual separately and the resulting output is a compound `tibble` with rows corresponding to each individual `ssm_df` fit object. The `converged` column indicates whether each model fit converged successfully.
 
-```{r multi-fits, message=FALSE}
+```{r multi-fits, message=FALSE, eval=FALSE}
 ## fit to data with two individuals
 fit <- fit_ssm(sese2,
                model = "crw",
@@ -196,6 +196,11 @@ fit <- fit_ssm(sese2,
 ## list fit outcomes for both seals
 fit
 ```
+```{r multi-fits2, echo=FALSE}
+load("data/overview/fit.rda")
+fit
+```
+
 
 Individual `id` is displayed in the 1st column, all fit output resides in a list (`ssm`) in the 2nd column, `convergence` status (whether the optimizer found a global minimum) of each model fit is displayed in the 3rd column, whether the Hessian matrix was positive-definite and could be solved to obtain parameter standard errors (`pdHess`) is displayed in the 4th column, and the specified process model (`rw`, `crw`, or `mp`) in the 5th column. In some cases, the optimizer will converge but the Hessian matrix is not positive-definite, which typically indicates the optimizer converged on a local minimum. In this case, some standard errors may be calculated but not all. One possible solution is to try specifying a longer `time.step` or set `time.step = NA` to turn off predictions and return only fitted values (location estimates at the pre-filtered observation times). If `pdHess = FALSE` persists then careful inspection of the supplied data is warranted to determine if suspect observations not identified by `prefilter` are present. The excellent [glmmTMB troubleshooting vignette](https://CRAN.R-project.org/package=glmmTMB/vignettes/troubleshooting.html) may also provide hints at solutions. Convergence failures should be examined for potential data issues, however, in some cases changes to the optimization parameters via `ssm_control()` (see `?fit_ssm` and `?ssm_control` on usage) may overcome mild issues (see `?nlminb` or `?optim` for details on optimization control parameters).
 

vignettes/SSM_fitting.Rmd

@@ -26,7 +26,7 @@ Fitting SSM's to tracking data in `aniMotum` is intended to be as simple as poss
 4. using [`trip::sda`](https://CRAN.R-project.org/package=trip) to identify potentially extreme observations based on speed, distance, and angle (Freitas et al. 2008).  
 
 Several arguments to `fit_ssm` provide a measure of control over data processing. The `vmax`, `ang` and `distlim` arguments provide control over the `trip::sda` filter (see `?fit_ssm` for details). The `sda` filter can occasionally fail (possibly when calculating angles between successive displacements), in this case we fall back to using `trip::speedfilter` and a warning issued. The `trip` filter can be turned off altogether with `spdf = FALSE`, although this is usually not advisable if working with older Least-Squares Argos data, the lower prevalence of extreme observations in newer Kalman filtered Argos data may not always require speed-filtering prior to fitting an SSM. The minimum time interval allowed between observations is set by `min.dt`. Finally, the processed data can be returned, e.g. for careful inspection, without fitting an SSM with the `pf = TRUE` argument:
-```{r pre-processing, message=FALSE}
+```{r pre-processing, eval=FALSE, message=FALSE}
 fit_ssm(ellie, 
         vmax = 3, 
         pf = TRUE)
@@ -52,8 +52,8 @@ The default emf values are given by `emf()`:
 emf()
 ```
 These can be adjusted by modifying the data.frame and passing this to `fit_ssm` via `emf.
-```{r emf mod, message=FALSE, fig.width=7, fig.height=6}
-emf2 <- mutate(emf(), emf.x = c(0.1, 1, 2, 4, 20, 8, 15, 15))
+```{r emf mod, eval=FALSE, message=FALSE}
+emf2 <- dplyr::mutate(emf(), emf.x = c(0.1, 1, 2, 4, 20, 8, 15, 15))
 
 fit1 <- fit_ssm(subset(sese2, id == unique(id)[1]), 
                 model = "rw", 
@@ -74,6 +74,7 @@ hist(grab(fit1, "f")$x - grab(fit2, "f")$x, breaks = 20,
 hist(grab(fit1, "f")$y - grab(fit2, "f")$y, breaks = 20, 
      main = "fit1$y - fit2$y", xlim = c(-20,20), xlab = "Distance (km)")
 ```
+![](images/ssm_fitting/histos.jpeg){width="98%"}
 
 In this contrived example, re-specified `emf.x` values result in differences in the predicted locations of up to about 20 km in the x-direction and a more subtle 4 km in the y-direction, even though the `emf.y` values were not edited. The latter is due to a model-estimated correlation between the x and y errors.
 
@@ -83,7 +84,7 @@ Now that CLS Argos provides locations estimated by their Kalman filter algorithm
 The `map` argument allows some model parameters to by turned off by fixing them at 0 on the log-scale. This can be useful in cases where model simplification can aid in convergence or simply provide a better fit to the data. For example, the `crw` has a parameter `psi` that re-scales the semi-major axis of Argos KF error ellipses to account for apparent bias in the error ellipses (Jonsen et al. 2020). Occasionally, there may be insufficient information to reliably estimate this parameter, or one might wish to compare model fits with and without error ellipse bias correction.
 
 Here we fit the `crw` model with (`fit1`) and without (`fit2`) the bias-correction parameter `psi`, by supplying a named list to the `map` argument. Where `factor(NA)` is a `TMB` convention that ensures `psi` is fixed during optimization.
-```{r demonstrate map use, message=FALSE, warning=FALSE}
+```{r demonstrate map use, eval=FALSE, message=FALSE, warning=FALSE}
 fit1 <- fit_ssm(ellie, 
                 model = "crw", 
                 time.step = 24,
@@ -97,6 +98,10 @@ fit2 <- fit_ssm(ellie,
 
 c(fit1$ssm[[1]]$AICc, fit2$ssm[[1]]$AICc)
 ```
+```{r demonstrate map use2, message=FALSE, echo=FALSE}
+load("data/ssm_fitting/fits_aic.rda")
+c(fit1$ssm[[1]]$AICc, fit2$ssm[[1]]$AICc)
+```
 
 Comparison of AIC~c~ from the model fits implies the model including the `psi` parameter `fit1` provides a slightly better fit.
 

vignettes/SSM_validation.Rmd

@@ -17,30 +17,33 @@ require(aniMotum)
 
 Just like any other statistical model fitting exercise, validating SSM fits to tracking data is an essential part of any analysis. SSM fits can be visualized quickly using the generic `plot` function on model fit objects:
 
-```{r visualise SSM fit1, eval=TRUE, fig.width=7, fig.height=5}
+```{r visualise SSM fit1, eval=FALSE}
 fit.rw <- fit_ssm(ellie, 
                   model = "rw", 
                   time.step = 24, 
                   control = ssm_control(verbose = 0))
-
+```
+```{r fit_rw_f, eval=FALSE, fig.width=7, fig.height=5}
 plot(fit.rw, what = "fitted")
-
+```
+![](images/ssm_validation/fit_rw_f.jpeg){width="98%"}
+```{r fit_rw_p, eval=FALSE, fig.width=7, fig.height=5}
 plot(fit.rw, what = "predicted")
 ```
-
+![](images/ssm_validation/fit_rw_p.jpeg){width="98%"}
 
 Both fitted and predicted locations (gold) can be plotted as 1-D time-series on top of the observations (blue) to visually detect any lack of fit. Observations that failed to pass the `prefilter` stage prior to SSM fitting (black x's) can be included (default) or omitted with the `outlier` argument. Uncertainty is displayed as a $\pm$ 2 SE envelope (pale gold) around estimates. A rug plot along the x-axis aids detection of data gaps in the time-series. Note, in second plot, the larger standard errors for predicted locations through small data gaps.
 
 The SSM fits can also be visualised as 2-D tracks via the `type` argument:
-```{r visualise as track, eval=TRUE, fig.width=7, fig.height=6}
+```{r visualise as track, eval=FALSE}
 plot(fit.rw, "p", type = 2, alpha = 0.1)
 ```
-
+![](images/ssm_validation/fit_rw_track.jpeg){width="98%"}
 
 This option provides an intuitive view of the estimated track (gold) through the observations (blue), along with a 2-dimensional representation of the location uncertainty (pale gold 95% confidence ellipses). Here, we use the `alpha` argument to increase the transparency of the confidence ellipses to aid visualization in regions where many overlay one another (upper left and lower right in figure).
 
 Residual plots are important for validating models, but classical Pearson residuals, for example, are not appropriate for state-space models. Instead, a one-step-ahead prediction residual, provides a useful if computationally demanding alternative. In `aniMotum`, prediction residuals from state-space model fits are calculated using the `osar` function and can be visualized as time-series plots, Q-Q plots, or autocorrelation functions:
-```{r osar plots1, eval=TRUE, warning=FALSE, message=FALSE, fig.width=7, fig.height=6}
+```{r osar plots1, eval=FALSE, warning=FALSE, message=FALSE}
 # use patchwork package to arrange plot.osar options
 require(patchwork)
 # calculate & plot residuals
@@ -49,12 +52,12 @@ res.rw <- osar(fit.rw)
 (plot(res.rw, type = "ts") | plot(res.rw, type = "qq")) / 
   (plot(res.rw, type = "acf") | plot_spacer())
 ```
-
+![](images/ssm_validation/osar_rw_plots.jpg){width="98%"}
 
 Here, the 3 residual plots highlight a poor fit of the `rw` SSM to the x-coordinate. The time-series residual plot implies a bias in the `rw` SSM fit over the latter half of the track, the Q-Q plot highlights a distinct departure from normality in the x residuals, and the x residuals are positively autocorrelated up to about lag 5.
 
-Fitting the `crw` SSM to the same data, we can see the prediction residual plots imply a less biased fit with approximately normal residuals that have little autocorrelation in both coordinates. Additionally, the `crw` model fit has a lower AIC~c~ than the `rw`. Note that AIC statistics can be misleading for time-series models and should not be used as the sole criterion for preferring one model fit over another. Here, at least, the AIC~c~ values are in agreement with the prediction residual diagnostics.
-```{r fit crw and plot osar1, eval=TRUE, warning=FALSE, message=FALSE, fig.width=7, fig.height=6}
+Fitting the `crw` SSM to the same data, we can see the prediction residual plots imply a less biased fit with approximately normal residuals that have little autocorrelation in both coordinates. 
+```{r fit crw and plot osar1, eval=FALSE, warning=FALSE, message=FALSE}
 fit.crw <- fit_ssm(ellie, 
                    model = "crw", 
                    time.step = 24, 
@@ -64,10 +67,16 @@ res.crw <- osar(fit.crw)
 
 (plot(res.crw, type = "ts") | plot(res.crw, type = "qq")) / 
   (plot(res.crw, type = "acf") | plot_spacer())
-
+```
+![](images/ssm_validation/osar_crw_plots.jpg){width="98%"}
+Additionally, the `crw` model fit has a lower AIC~c~ than the `rw` model. Note that AIC statistics can be misleading for time-series models and should not be used as the sole criterion for preferring one model fit over another. Here, at least, the AIC~c~ values are in agreement with the prediction residual diagnostics.
+```{r AICa, eval=FALSE}
+c(fit.rw$ssm[[1]]$AICc, fit.crw$ssm[[1]]$AICc)
+```
+```{r AICb, echo=FALSE}
+load("data/ssm_validation/fits.rda")
 c(fit.rw$ssm[[1]]$AICc, fit.crw$ssm[[1]]$AICc)
 ```
-
 
 As calculation of prediction residuals can be computationally demanding, typically requiring more time than fitting the model, especially for multiple individual tracks, the `osar` function is automatically implemented in parallel when calculating residuals for more than 2 tracks. 
 

vignettes/Track_simulation.Rmd

@@ -25,7 +25,7 @@ The `sim` function is useful for situations where simulated random tracks are re
 
 #### Simulate tracks from each of the 3 process models
 
-```{r sim}
+```{r sim, eval=FALSE}
 set.seed(pi) 
 sim.rw <- sim(N = 200, model = "rw")
 sim.crw <- sim(N = 200, model = "crw", D = 0.5)
@@ -44,8 +44,7 @@ sim.mp <- sim(N = 200, model = "mp", sigma_g = 0.3)
 ![](images/track_simulation/plot_sim1.jpg){width="95%"}
 
 #### Fit SSM to simulated lon,lat locations (with Argos error)
-
-```{r SSM, warning=FALSE, message=FALSE, fig.width=7, fig.height=6}
+```{r SSM, eval=FALSE, warning=FALSE, message=FALSE}
 # coerce simulated RW data to format expected by fit_ssm
 d <- with(sim.rw, data.frame(id = 1, date, lc, lon, lat))