Bad models can still generate good indices

[ericward-noaa]

I've been having a email discussion with some colleagues and want to highlight a point -- that even models that might be missing important features, like covariates in the main effects, can still generate decent indices of abundance (largely because these models are very flexible when intercepts / spatial and spatiotemporal fields are included). Example below:

library(sdmTMB)

set.seed(123)

# make fake data for locations -- this includes joining 2 data frames, one for observation locations 
# and one for the prediction grid. this doesn't have to be done in one step like this -- but want to 
# make sure spatial effects are consistent across both
n_loc <- 500
n_year <- 6
predictor_dat <- data.frame(
  X = runif(n_loc*n_year), Y = runif(n_loc*n_year),
  year = rep(1:n_year, each = n_loc), type = "obs"
)
# make grid for predicting 
grid_dat <- expand.grid(X = seq(0,1,length.out=100),
                        Y = seq(0,1,length.out=100),
                        year = unique(predictor_dat$year),
                        type = "grid")
all_dat <- rbind(grid_dat, predictor_dat)
mesh <- make_mesh(all_dat, xy_cols = c("X", "Y"), cutoff = 0.1)

year_effects <- rnorm(length(unique(predictor_dat$year)), 3, 1)
all_dat$fyear <- as.factor(all_dat$year)

# make a field for simulated covariate -- this has mean 1
coef_field <- sdmTMB_simulate(
  formula = ~ 1,
  data = all_dat,
  time = "year",
  mesh = mesh,
  family = gaussian(),
  range = 0.1,
  sigma_E = 0.0,
  phi = 0.1,
  sigma_O = 0.2,
  seed = 42,
  B = c(1) # B0 = intercept
)
all_dat$a1 <- coef_field$mu # don't include obs error

sim_dat <- sdmTMB_simulate(
  formula = ~ -1 + fyear + a1 + I(a1^2),
  data = all_dat,
  time = "year",
  mesh = mesh,
  family = lognormal(),#gaussian(),
  range = 0.5,
  sigma_E = 0.1,
  phi = 0.1,
  sigma_O = 0.2,
  seed = 42,
  B = c(year_effects, -0.4, 0.03) # B0 = intercept, B1 = a1 slope
)

all_dat <- dplyr::left_join(all_dat, sim_dat)
#> Joining with `by = join_by(X, Y, year, a1)`

# fit the model to the observed data. first use a small sample of the above
obs_dat <- dplyr::filter(all_dat, type != "grid")
sub_dat <- obs_dat[sample(1:nrow(obs_dat), size=300, replace=F),]
mesh <- make_mesh(sub_dat, xy_cols = c("X", "Y"), cutoff = 0.14)
fit1 <- sdmTMB(observed ~ -1 + as.factor(year),
              time = "year",
              mesh <- mesh,
              spatiotemporal = "iid",
              family = lognormal(),
              data = sub_dat)
# Fit simple model with 'a1' included as a quadratic effect
fit2 <- update(fit1, formula = observed ~ -1 + as.factor(year) + a1 + I(a1^2))

# Model 2 -- with the covariate is better from AIC standpoint
AIC(fit1, fit2)
#>      df      AIC
#> fit1 10 1204.623
#> fit2 12 1184.564

# But when it comes to making predictions, both models generate the same index
# First the model with no covariate
pred1 <- predict(fit1, newdata = dplyr::filter(all_dat, type=="grid"), 
                 return_tmb_object = TRUE)
ind1 <- get_index(pred1, bias_correct = TRUE)
ind1
#>   year        est        lwr        upr  log_est         se
#> 1    1  109446.57  105353.22  113698.96 11.60319 0.01944820
#> 2    2   89404.59   86252.79   92671.56 11.40093 0.01831135
#> 3    3   30202.22   29085.38   31361.95 10.31567 0.01922480
#> 4    4   65816.13   63665.14   68039.79 11.09462 0.01695328
#> 5    5 1783270.96 1723879.94 1844708.12 14.39396 0.01728183
#> 6    6  126249.20  122035.65  130608.23 11.74601 0.01731896

# And now the model with the covariate
pred2 <- predict(fit2, newdata = dplyr::filter(all_dat, type=="grid"), 
                 return_tmb_object = TRUE)
ind2 <- get_index(pred2, bias_correct = TRUE)
ind2
#>   year        est        lwr        upr  log_est         se
#> 1    1  109572.13  105622.53  113669.43 11.60434 0.01873066
#> 2    2   89414.73   86370.34   92566.43 11.40104 0.01767436
#> 3    3   30265.85   29177.32   31394.98 10.31778 0.01868817
#> 4    4   65403.56   63334.63   67540.08 11.08833 0.01640054
#> 5    5 1776144.20 1718771.49 1835432.01 14.38996 0.01675288
#> 6    6  125703.85  121652.17  129890.47 11.74168 0.01671604

The main takeaway being that the indices are perfectly correlated and of the same scale -- both models include year as a fixed effect, but one omits an important quadratic covariate and the other includes it.

^{Created on 2024-02-22 with reprex v2.1.0}