Objective: demonstrate how to calculate standard errors for functions of the movement parameters using the delta method. Throughout this document, we will use the data set from Chatterjee et al. (2024), which has been archived in the Data Repository of the University of Minnesota at the following link.

Setup

Load libraries and data

library(raster)
library(ggplot2) # plotting data
library(cowplot) # multi-panel plots
library(here) # for getting the relative filepaths
library(msm) # for delta-method calculation
library(mixedSSA) # for updating distribution parameters
library(dplyr)

Source helper functions to fit the models

source(here("Delta_method/helper_functions.R"))
## Warning: package 'glmmTMB' was built under R version 4.3.3
## Warning: package 'knitr' was built under R version 4.3.3
### Load the data set
### You can load it directly from the online repository
snapper_dat <- readRDS(url("https://conservancy.umn.edu/bitstreams/afee8481-cc62-4729-a574-7563b75815da/download")) |>
     rename(  # rename the columns to match appendix
         hardbottom_start = reefS4, # dummy variable for start of step reef class
         hardbottom_end = reefE4, # dummy variable for end of step reef class
         distance_to_reefedge = newdist
       )

### Alternatively, you can also load it from your local system 
snapper_dat <- readRDS("Snapper_final_data_6min.rds") |>
  rename(  # rename the columns to match appendix
    hardbottom_start = reefS4, # dummy variable for start of step reef class
    hardbottom_end = reefE4, # dummy variable for end of step reef class
    distance_to_reefedge = newdist
  )
glimpse(snapper_dat) # Data set used in the ISSA analysis...should have 12 columns
## Rows: 465,903
## Columns: 16
## $ id                   <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ burst_               <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5…
## $ sl_                  <dbl> 9.9313927, 0.3451452, 16.2750112, 10.6542033, 7.4…
## $ ta_                  <dbl> -1.82151359, 1.03694629, -0.46531514, -1.30649150…
## $ dt_                  <drtn> 5.283333 mins, 5.283333 mins, 5.283333 mins, 5.2…
## $ step_id_             <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2…
## $ case_                <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, …
## $ dist_edge            <dbl> 4.472136, 5.000000, 0.000000, 0.000000, 11.180340…
## $ reefStart            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4…
## $ reefEnd              <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 1, 4, 1…
## $ log_sl_              <dbl> 2.2957007, -1.0637900, 2.7896309, 2.3659545, 2.00…
## $ cos_ta_              <dbl> -0.24809886, 0.50885143, 0.89368020, 0.26123830, …
## $ step_id              <chr> "2_1", "2_1", "2_1", "2_1", "2_1", "2_1", "2_1", …
## $ hardbottom_start     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1…
## $ hardbottom_end       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0…
## $ distance_to_reefedge <dbl> 4.472136, 5.000000, 0.000000, 0.000000, 11.180340…

Model 1 - Categorical covariate: Benthic habitat class

The benthic habitat, or reef, class is represented as a dummy variable indicating whether the habitat is Sand (0) or Hardbottom (1).

Fit the model

model1 <- fit_model(
  data = snapper_dat,
  model_name = "snapper_issa_6min_categorical_new",
  interaction_var_name = "hardbottom_start", # benthic habitat class at the start of each step
  save_model = TRUE,
  from_cache = TRUE # switch to FALSE, if you have not yet fit the model
)

summary(model1)
##  Family: poisson  ( log )
## Formula:          
## case_ ~ sl_ + log_sl_ + cos(ta_) + hardbottom_end + distance_to_reefedge +  
##     hardbottom_start:sl_ + hardbottom_start:log_sl_ + hardbottom_start:cos(ta_) +  
##     (1 | step_id) + (0 + sl_ + log_sl_ + cos(ta_) | id) + (0 +  
##     distance_to_reefedge | id) + (0 + hardbottom_end | id)
## Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
## 182377.9 182576.8 -91170.9 182341.9   465479 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name                 Variance  Std.Dev.  Corr        
##  step_id (Intercept)          1.266e-08 0.0001125             
##  id      sl_                  3.949e-04 0.0198722             
##          log_sl_              3.062e-02 0.1749752 -0.95       
##          cos(ta_)             3.309e-02 0.1818966  0.25 -0.19 
##  id.1    distance_to_reefedge 4.556e-04 0.0213459             
##  id.2    hardbottom_end       4.330e-02 0.2080790             
## Number of obs: 465488, groups:  step_id, 45915; id, 35
## 
## Conditional model:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -3.116218   0.019086 -163.27  < 2e-16 ***
## sl_                       -0.004719   0.003620   -1.30  0.19239    
## log_sl_                    0.107612   0.033393    3.22  0.00127 ** 
## cos(ta_)                  -0.854549   0.036955  -23.12  < 2e-16 ***
## hardbottom_end             0.246946   0.043437    5.69 1.31e-08 ***
## distance_to_reefedge      -0.013068   0.004001   -3.27  0.00109 ** 
## sl_:hardbottom_start      -0.014503   0.001205  -12.04  < 2e-16 ***
## log_sl_:hardbottom_start   0.052904   0.012030    4.40 1.09e-05 ***
## cos(ta_):hardbottom_start  0.036143   0.020255    1.78  0.07435 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Update Step-length (gamma) distribution using model coefficients

We will use the update_dist function in the mixedSSA package, which will handle the updating for us.

# We will need to pass the tentative sl distribution
tentative_sl_dist <- amt::fit_distr(snapper_dat$sl_[snapper_dat$case_ == "TRUE"],
  "gamma",
  na.rm = T
)


# Call the function to get the updated gamma distribution parameter values
model1$frame <- model1$frame |>
  mutate(
    hardbottom_start = as.factor(hardbottom_start)
  )

updated_param_sl_cat <- mixedSSA::update_dist(
  model = model1,
  dist_name = "gamma",
  beta_sl = "sl_", # the name of the step length coefficient in our model
  beta_log_sl = "log_sl_", # the name of the log(step length) coefficient in our model
  interaction_var_name = "hardbottom_start", # benthic habitat class at the start of each step
  tentative_dist = tentative_sl_dist,
  quantiles = NULL
)


# mutate the returned object for ease of plotting legend correctly
updated_param_sl_cat@updated_parameters <- updated_param_sl_cat@updated_parameters |>
  mutate(
    grouping = case_when(
      grouping == "tentative" ~ "tentative",
      grouping == 0 ~ "Sand",
      grouping == 1 ~ "Hardbottom"
    )
  )

# print the updated parameters
updated_param_sl_cat@updated_parameters
##     grouping random_effect   beta_sl beta_log_sl    shape    scale
## 1  tentative            NA        NA          NA 1.090138 18.23989
## 2       Sand            NA -0.004719    0.107612 1.197750 16.79434
## 3 Hardbottom            NA -0.019222    0.160516 1.250654 13.50503

Plot updated gamma distributions

To do this, we will use the plot_updated_dist function in the mixedSSA package, which will handle plotting for us.

plot_sl_cat <- mixedSSA::plot_updated_dist(
  updated_dist_params_obj = updated_param_sl_cat,
  include_tentative = FALSE,
  print_plot = FALSE,
  xlim = 100
) +
  scale_color_manual( # add some custom plotting
    name = "Habitat class",
    values = c("Sand" = "brown", "Hardbottom" = "orange")
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.8, 0.8),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 10)
  )

Update turn-angle (von Mises) distribution using model coefficients

# for turn-angle distributions in mixedSSA, we need to pass the tentative distribution
tentative_ta_dist <- amt::fit_distr(snapper_dat$ta_[snapper_dat$case_ == "TRUE"],
  "vonmises",
  na.rm = T
)

updated_param_ta_cat <- mixedSSA::update_dist(
  model = model1,
  dist_name = "vonmises",
  beta_cos_ta = "cos(ta_)", # the name of the cos(ta) coefficient in our model
  interaction_var_name = "hardbottom_start",
  tentative_dist = tentative_ta_dist
)

# mutate the returned object for ease of plotting legend correctly
updated_param_ta_cat@updated_parameters <- updated_param_ta_cat@updated_parameters |>
  mutate(
    grouping = case_when(
      grouping == "tentative" ~ "tentative",
      grouping == 0 ~ "Sand",
      grouping == 1 ~ "Hardbottom"
    )
  )

updated_param_ta_cat@updated_parameters
##     grouping random_effect beta_cos_ta     kappa mu
## 1  tentative            NA          NA  0.348415  0
## 2       Sand            NA   -0.854549 -0.506134  0
## 3 Hardbottom            NA   -0.818406 -0.469991  0

Plot updated von Mises distributions

plot_ta_cat <- mixedSSA::plot_updated_dist(
  updated_dist_params_obj = updated_param_ta_cat,
  vonmises_mu = pi,
  include_tentative = FALSE,
  print_plot = FALSE
) +
  scale_color_manual(
    name = "Habitat class",
    values = c("Sand" = "brown", "Hardbottom" = "orange")
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.5, 0.8),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 10)
  )

Delta method to estimate SE for mean speed for each category of the hardbottom_start variable

The msm package used for the delta method calculation requires:

  • a vector of estimated coefficients (which we save as estmean below)
  • an estimate of their variance-covariance matrix (which we save as estvar below)
  • the function of the parameters for which we want to estimate the standard error (we will discuss this in more detail below)
estmean <- summary(model1)$coefficients$cond[, 1]
estvar <- vcov(model1)$cond

For the gamma distribution, the mean is equal to the product of the scale and shape parameters. Thus, the function we need to pass to msm = scale x shape, but where the scale and shape parameters are the updated versions involving the regression coefficients from the fitted model.

  • updated scale parameters require the coefficients associated with the step lengths, sl_, and their interactions
  • updated shape parameters require the coefficients associated with the log(step lengths), log_sl_, and their interactions

We need to tell the deltamethod function in the msm package which coefficients from the fitted model are needed to calculate the updated scale and shape parameters, and thus, the mean of the updated gamma distribution. We do this using xi, where the i is replaced with the index needed to pull off the correct term in the estimated coefficient vector. To understand this, let’s look at the coefficients from the model:

estmean
##               (Intercept)                       sl_                   log_sl_ 
##              -3.116217777              -0.004718976               0.107612099 
##                  cos(ta_)            hardbottom_end      distance_to_reefedge 
##              -0.854548616               0.246946149              -0.013068445 
##      sl_:hardbottom_start  log_sl_:hardbottom_start cos(ta_):hardbottom_start 
##              -0.014502606               0.052904020               0.036142682

We can use x2 to refer to the coefficient associated with step length (sl_) and x7 to refer to the coefficient associated with the interaction between step length (sl_) and hardbottom_start. These coefficients are used to update the scale parameter. We will also need to refer to x3 and x8, coefficients associated with log(step length) (log_sl_) and its interaction with hardbottom_start, when updating the shape parameter.

### tentative distribution parameters
cat_sl_dist <- updated_param_sl_cat@updated_parameters
tentative_shape <- cat_sl_dist[1, "shape"]
tentative_scale <- cat_sl_dist[1, "scale"]

Calculating SEs: Note that the deltamethod calculates the variance of the function of our parameters. We then take the sqrt to get the SEs. The functions of the model parameters that we want to estimate (e.g. ~ (tentative_shape + x3) * (1 / (1 / tentative_scale - x2))) in the delta method calculation uses the corrected equations 9, 10, 12, and 13 of the updated version of the manuscript with all random effect terms (i.e., the \(b\)’s) set to 0. These equations were originally incorrectly specified in Appendix 1.

# filter out "tentative", we won't calculate a SE for the tentative parameters
cat_sl_dist <- cat_sl_dist |>
  filter( 
    grouping != "tentative"
    )

# SE for Sand Habitat (hardbottom_start == 0)
cat_sl_dist$se[1] <- sqrt(
  deltamethod(
    ~ (tentative_shape + x3) * (1 / (1 / tentative_scale - x2)),
    estmean, 
    estvar
    )
  )


# SE for Hardbottom Habitat (hardbottom_start == 1)
cat_sl_dist$se[2] <- sqrt(
  deltamethod(
    ~ (tentative_shape + x3 + x8) * (1 / (1 / tentative_scale - (x2 + x7))), 
    estmean, 
    estvar
    )
  )

### check the final output
cat_sl_dist
##     grouping random_effect   beta_sl beta_log_sl    shape    scale        se
## 1       Sand            NA -0.004719    0.107612 1.197750 16.79434 0.8569052
## 2 Hardbottom            NA -0.019222    0.160516 1.250654 13.50503 0.6727631
### Plot the results
plot_data <- cat_sl_dist |>
  mutate(
    mean = shape * scale,
    ci_lwr = mean - 1.96 * se,
    ci_upr = mean + 1.96 * se
  )

### Check results for mean speed
plot_data
##     grouping random_effect   beta_sl beta_log_sl    shape    scale        se
## 1       Sand            NA -0.004719    0.107612 1.197750 16.79434 0.8569052
## 2 Hardbottom            NA -0.019222    0.160516 1.250654 13.50503 0.6727631
##       mean   ci_lwr   ci_upr
## 1 20.11542 18.43589 21.79495
## 2 16.89012 15.57151 18.20874
### Plot
plot_delta_cat <- ggplot(plot_data) +
  geom_point(aes(x = grouping, y = mean), size = 2.5, col = "blue") +
  geom_linerange(aes(
    x = grouping,
    ymin = ci_lwr,
    ymax = ci_upr
  ), linewidth = 1) +
  labs(x = "Benthic habitat class", y = "Mean step length") +
  theme_bw() +
  ylim(c(15, 25)) +
  theme(
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 10)
  )

Model 2 - Continuous covariate: Distance to reef edge (distance_to_reefedge)

Fit the model

model2 <- fit_model(
  data = snapper_dat,
  model_name = "snapper_issa_6min_continuous_new",
  interaction_var_name = "distance_to_reefedge",
  save_model = TRUE,
  from_cache = TRUE
)

summary(model2)
##  Family: poisson  ( log )
## Formula:          
## case_ ~ sl_ + log_sl_ + cos(ta_) + hardbottom_end + distance_to_reefedge +  
##     distance_to_reefedge:sl_ + distance_to_reefedge:log_sl_ +  
##     distance_to_reefedge:cos(ta_) + (1 | step_id) + (0 + sl_ +  
##     log_sl_ + cos(ta_) | id) + (0 + distance_to_reefedge | id) +  
##     (0 + hardbottom_end | id)
## Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
## 182695.3 182894.3 -91329.7 182659.3   465577 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name                 Variance  Std.Dev.  Corr        
##  step_id (Intercept)          1.266e-08 0.0001125             
##  id      sl_                  4.269e-04 0.0206625             
##          log_sl_              3.274e-02 0.1809387 -0.94       
##          cos(ta_)             3.001e-02 0.1732364  0.32 -0.26 
##  id.1    distance_to_reefedge 5.137e-04 0.0226652             
##  id.2    hardbottom_end       4.110e-02 0.2027392             
## Number of obs: 465586, groups:  step_id, 45931; id, 35
## 
## Conditional model:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -3.047e+00  2.132e-02 -142.93  < 2e-16 ***
## sl_                           -7.333e-03  3.760e-03   -1.95   0.0511 .  
## log_sl_                        7.201e-02  3.466e-02    2.08   0.0378 *  
## cos(ta_)                      -8.926e-01  3.548e-02  -25.16  < 2e-16 ***
## hardbottom_end                 2.346e-01  4.241e-02    5.53 3.19e-08 ***
## distance_to_reefedge          -2.295e-02  4.543e-03   -5.05 4.37e-07 ***
## sl_:distance_to_reefedge      -2.521e-04  4.628e-05   -5.45 5.08e-08 ***
## log_sl_:distance_to_reefedge   6.403e-03  8.876e-04    7.21 5.41e-13 ***
## cos(ta_):distance_to_reefedge  5.594e-03  7.813e-04    7.16 8.11e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Update Step-length (gamma) distribution using model coefficients

We will again use the update_dist function in the mixedSSA package, which will handle the updating for us. Because we are interacting a continuous variable with step length and log step length, we are effectively assuming the shape and scale paraemters of the gamma distribution vary (continuously) with distance from reefedge. To better understand this relationship, we will update the parameters for different quantiles of distance from reefedge.

### The quantiles we will update the parameters with
quantiles <- c(0.05, 0.5, 0.95)

updated_param_sl_cont <- mixedSSA::update_dist(
  model = model2,
  dist_name = "gamma",
  beta_sl = "sl_", # the name of the step lengths coefficient in our model
  beta_log_sl = "log_sl_", # the name of the log(step length) coefficient in our model
  interaction_var_name = "distance_to_reefedge", # distance to reef edge
  quantiles = quantiles, # for continuous interaction variables
  tentative_dist = tentative_sl_dist
)

Plot updated gamma distributions

We will again use the plot_updated_dist function in the 1mixedSSA` package for plotting, passing it our newly updated parameters.

plot_sl_cont <- mixedSSA::plot_updated_dist(
  updated_dist_params_obj = updated_param_sl_cont,
  include_tentative = FALSE,
  print_plot = FALSE,
  xlim = 100
) +
  scale_color_manual(
    name = "Distance to\nreef edge",
    labels = c("Close", "Medium", "Far"),
    values = c("orange", "red", "brown")
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.8, 0.8),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 10)
  )

Update turn-angle (von Mises) distribution using model coefficients

Again, because this is a von Mises distribution we need to pass the tentative distribution of the turn angles. We will pass the same tentative distribution from before.

updated_param_ta_cont <- mixedSSA::update_dist(
  model = model2,
  dist_name = "vonmises",
  beta_cos_ta = "cos(ta_)", # the name of the cos(ta) coefficient in our model
  interaction_var_name = "distance_to_reefedge",
  quantiles = quantiles, # for continuous interaction variables
  tentative_dist = tentative_ta_dist
)

Plot updated von Mises distributions

plot_ta_cont <- mixedSSA::plot_updated_dist(
  updated_dist_params_obj = updated_param_ta_cont,
  vonmises_mu = pi,
  include_tentative = FALSE,
  print_plot = FALSE
) +
  scale_color_manual(
    name = "Distance to\nreef edge",
    labels = c("Close", "Medium", "Far"),
    values = c("orange", "red", "brown")
  ) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(0.5, 0.8),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 10)
  )

Delta method for SE calculation of mean speed

This will follow the same process as before, but we will calculate the mean speed (and its SE) for different quantiles of distance from reefedge. The function of the model parameters that we want to estimate (i.e., ~ (tentative_shape + x3 + x8 * distanceq) * (1 / (1 / tentative_scale - (x2 + x7 * distanceq)))) in the delta method calculation uses equations 17 and 18 in the updated manuscript version, which were originally mis-specified in Appendix 1.

### Estimated coefficients and their var/cov matrix
contmean <- summary(model2)$coefficients$cond[, "Estimate"]
contvar <- vcov(model2)$cond

### tentative distribution parameters
cont_sl_dist <- updated_param_sl_cont@updated_parameters
tentative_shape <- cont_sl_dist[1, "shape"]
tentative_scale <- cont_sl_dist[1, "scale"]

# filter out "tentative", we won't calculate a SE for the tentative parameters
cont_sl_dist <- cont_sl_dist |>
  filter(
    grouping != "tentative"
  ) |>
  mutate(
    # get the actual values of quantiles for delta method
    dist_quant = quantile(snapper_dat$distance_to_reefedge,
                          na.rm = T,
                          probs = as.numeric(grouping)
    )
  )


for (i in 1:nrow(cont_sl_dist)) {
  distanceq <- cont_sl_dist$dist_quant[i] # pull of the distance value
  cont_sl_dist$se[i] <- sqrt(
    deltamethod(
      ~ (tentative_shape + x3 + x8 * distanceq) * (1 / (1 / tentative_scale - (x2 + x7 * distanceq))),
      contmean,
      contvar
      )
    )
}


# check output
cont_sl_dist
##   grouping random_effect   beta_sl beta_log_sl    shape    scale dist_quant
## 1     0.05            NA -0.007333    0.072007 1.162145 16.08812    0.00000
## 2      0.5            NA -0.008089    0.091217 1.181355 15.89469    3.00000
## 3     0.95            NA -0.017374    0.327022 1.417160 13.85057   39.82462
##          se
## 1 0.8101875
## 2 0.8066155
## 3 0.8020715
plot_data <- cont_sl_dist |>
  mutate(
    mean = shape * scale,
    ci_lwr = mean - 1.96 * se,
    ci_upr = mean + 1.96 * se
  )

### Check results for mean speed
plot_data
##   grouping random_effect   beta_sl beta_log_sl    shape    scale dist_quant
## 1     0.05            NA -0.007333    0.072007 1.162145 16.08812    0.00000
## 2      0.5            NA -0.008089    0.091217 1.181355 15.89469    3.00000
## 3     0.95            NA -0.017374    0.327022 1.417160 13.85057   39.82462
##          se     mean   ci_lwr   ci_upr
## 1 0.8101875 18.69673 17.10876 20.28469
## 2 0.8066155 18.77727 17.19630 20.35823
## 3 0.8020715 19.62848 18.05642 21.20054
### plotting the results
plot_delta_sd <- ggplot(plot_data) +
  geom_point(aes(x = dist_quant, y = mean), size = 2.5, col = "blue") +
  geom_linerange(aes(x = dist_quant, ymin = ci_lwr, ymax = ci_upr), linewidth = 1) +
  labs(x = "Distance to reef edge (m)", y = "Mean step length") +
  theme_bw() +
  ylim(c(15, 25)) +
  theme(
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 10)
  )

Final plot used in the publication

This is the final plot with the updated figure with delta method along with the different movement parameter distribution updated with mixedSSA

figure5 <- cowplot::plot_grid(
  plot_sl_cat, plot_delta_cat, plot_ta_cat,
  plot_sl_cont, plot_delta_sd, plot_ta_cont,
  ncol = 3, labels = c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)"),
  hjust = -0.1
)

cowplot::save_plot(
  filename = here("Delta_method/figure5.jpeg"),
  plot = figure5, base_height = 8, dpi=300,
)
knitr::include_graphics(here("Delta_method/figure5.jpeg"))