When to include a spatially-varying covariate in the formula argument?

[Lewis-Barnett-NOAA]

This is based on a question received by email and some personal experimentation.

In our sdmTMB workshop/teaching slides, we indicate that when fitting a model with a spatially-varying covariate, one will likely want to include the spatially-varying term in the formula argument in addition to the spatial_varying argument. However, this can result in a singularity error:

solve.default(h, g) :
system is computationally singular

Can we give more specific guidance on when to include the spatially varying term in the formula? If it is yet unclear, I can put together a MWE to better understand when this occurs.

[seananderson]

I think I'd need a MWE here. In general, I think you'd nearly always want a covariate in both since the spatially varying field is assumed to be drawn from a mean-zero distribution.

[Lewis-Barnett-NOAA]

That makes sense. I'm working on an example, but I think the issue is more to do with how the intercepts are specified than whether the spatially varying term is in the formula. If I can make a clean example maybe we can add it to the documentation, as the only examples we have now are for the special case of the spatially-varying coefficient of time.

[seananderson]

I will admit how the intercept in the spatially varying formula is handled is a mess right now and I've been confused myself. It may need some rethinking.

[Lewis-Barnett-NOAA]

Yes, I just got confused about handling intercepts in the spatially varying formula. In playing with examples, I haven't been able to reproduce the error that solve() was throwing in a more complicated case, which was:

Error in solve.default(h, g) :
Lapack routine dgesv: system is exactly singular: U[82,82] = 0

where there were 41 years included in a delta model.

However, I think that this error is caused by the Hessian not being invertible, which I can reproduce with an improperly structured model (m3 below). However, this does not occur in models 1 and 2 without the additional covariate, as these produce the same model fit regardless of the improper handling of the intercepts in model 2. Also, in this case solve() does not throw the error, but the Hessian is not positive definite, as reported by sanity(). Bottom line I guess is that there are a lot of ways to improperly specify these spatially varying coefficient models!

pcod_spde <- make_mesh(pcod, c("X", "Y"), cutoff = 12)

d <- pcod
d$scaled_year <- (pcod$year - mean(pcod$year)) / sd(pcod$year)

m1 <- sdmTMB(
density ~ 0 + scaled_year, data = d,
spatial_varying = ~ 0 + scaled_year,
mesh = pcod_spde,
family = delta_gamma(type = "poisson-link"),
time = "year",
extra_time = c(2006, 2008, 2010, 2012, 2014, 2016),
spatial = "on",
spatiotemporal = "ar1"
)

m2 <- sdmTMB(
density ~ 0 + scaled_year, data = d,
spatial_varying = ~ scaled_year,
mesh = pcod_spde,
family = delta_gamma(type = "poisson-link"),
time = "year",
extra_time = c(2006, 2008, 2010, 2012, 2014, 2016),
spatial = "on",
spatiotemporal = "ar1"
)

m3 <- sdmTMB(
density ~ 0 + depth_scaled + scaled_year, data = d,
spatial_varying = ~ 0 + scaled_year,
mesh = pcod_spde,
family = delta_gamma(type = "poisson-link"),
time = "year",
extra_time = c(2006, 2008, 2010, 2012, 2014, 2016),
spatial = "on",
spatiotemporal = "ar1"
)

[ericward-noaa]

Model m1 and m2 are identical, because the 0 in the spatial_varying formula is ignored.

One potential problem in folks interpreting m1 is that without a global intercept the mean of the spatial field is != 0 because it is compensating for that. Here it's just

pred_m1 <- predict(m1)
mean(pred_m1$omega_s1) # -0.24

Alternatively, one could fit the same model estimating the global intercept (with a goal of making the mean of omega_s 0)

m0 <- update(m1, formula = density ~ 1 + scaled_year)

But the mean of those spatial fields are still not 0, which is weird. The mean of omega_s is also non-zero when (1) the spatiotemporal field is iid, and (2) when the spatiotemporal fields are set to off

pred_m0 <- predict(m1)
mean(pred_m0$omega_s1) 

For model3, I think there's some similar issues arising with the spatiotemporal fields. If they're set to iid, they should all have mean 0. But if they're rw or ar1 like here, the means of each are != 0, and subject to a non-linear trend.

The bigger thing I'm not totally getting is why the means of omega_s and epsilon_st are never 0, even on a regularized grid. Using this example from the help files,

m <- sdmTMB(
  data = pcod,
  formula = density ~ 0 + as.factor(year),
  spatiotemporal = "iid",
  time = "year", mesh = pcod_spde, family = tweedie(link = "log")
)
nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
predictions <- predict(m, newdata = nd)

predictions %>%
  group_by(year) %>%
  summarize(m = mean(epsilon_st))

predictions %>%
  filter(year==2003) %>%
  summarize(m = mean(omega_s))

[seananderson]

The fields are drawn from MVN(0, Sigma) but any given realized field will not have mean zero, just like how random intercepts come from N(0, sigma), but any one random intercept value will not have a value of 0. So, even IID random fields will have different means through time, but they're penalized towards mean 0. And a given spatial field will be penalized towards zero but may not have a calculated mean zero.

[ericward-noaa]

Yeah, I get that -- I should have been more specific. I don't expect the random field means to be exactly 0. I'd expect them to be pretty tiny -- like 0.001. But in some of the examples above, the means are drifting around between 0.05 - 0.1. Maybe that's just to be expected.

These examples from the glmmTMB help have simple random effects in the intercept and have the kind of small means I was looking for

library(glmmTMB)
#> Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
#> glmmTMB was built with TMB version 1.9.6
#> Current TMB version is 1.9.10
#> Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
(m1 <- glmmTMB(count ~ mined + (1|site),
               zi=~mined,
               family=poisson, data=Salamanders))
#> Formula:          count ~ mined + (1 | site)
#> Zero inflation:         ~mined
#> Data: Salamanders
#>       AIC       BIC    logLik  df.resid 
#> 1908.4695 1930.8080 -949.2348       639 
#> Random-effects (co)variances:
#> 
#> Conditional model:
#>  Groups Name        Std.Dev.
#>  site   (Intercept) 0.28    
#> 
#> Number of obs: 644 / Conditional model: site, 23
#> 
#> Fixed Effects:
#> 
#> Conditional model:
#> (Intercept)      minedno  
#>      0.0879       1.1419  
#> 
#> Zero-inflation model:
#> (Intercept)      minedno  
#>       1.139       -1.736
mean(ranef(m1)$cond$site[[1]])
#> [1] 0.004552813

(m2 <- glmmTMB(count ~ spp + mined + (1|site),
               zi=~spp + mined,
               family=nbinom2, data=Salamanders))
#> Formula:          count ~ spp + mined + (1 | site)
#> Zero inflation:         ~spp + mined
#> Data: Salamanders
#>       AIC       BIC    logLik  df.resid 
#> 1670.2990 1750.7176 -817.1495       626 
#> Random-effects (co)variances:
#> 
#> Conditional model:
#>  Groups Name        Std.Dev.
#>  site   (Intercept) 0.3799  
#> 
#> Number of obs: 644 / Conditional model: site, 23
#> 
#> Dispersion parameter for nbinom2 family (): 1.51 
#> 
#> Fixed Effects:
#> 
#> Conditional model:
#> (Intercept)        sppPR        sppDM      sppEC-A      sppEC-L     sppDES-L  
#>     -0.6104      -0.9637       0.1707      -0.3871       0.4879       0.5895  
#>       sppDF      minedno  
#>     -0.1133       1.4294  
#> 
#> Zero-inflation model:
#> (Intercept)        sppPR        sppDM      sppEC-A      sppEC-L     sppDES-L  
#>      0.9100       1.1614      -0.9393       1.0424      -0.5623      -0.8930  
#>       sppDF      minedno  
#>     -2.5398      -2.5630
mean(ranef(m2)$cond$site[[1]])
#> [1] -0.002161247

(m3 <- glmmTMB(count ~ spp + mined + (1|site),
               zi=~spp + mined,
               family=truncated_poisson, data=Salamanders))
#> Formula:          count ~ spp + mined + (1 | site)
#> Zero inflation:         ~spp + mined
#> Data: Salamanders
#>       AIC       BIC    logLik  df.resid 
#> 1790.1669 1866.1177 -878.0834       627 
#> Random-effects (co)variances:
#> 
#> Conditional model:
#>  Groups Name        Std.Dev.
#>  site   (Intercept) 0.2306  
#> 
#> Number of obs: 644 / Conditional model: site, 23
#> 
#> Fixed Effects:
#> 
#> Conditional model:
#> (Intercept)        sppPR        sppDM      sppEC-A      sppEC-L     sppDES-L  
#>    -0.06702     -0.52093      0.22458     -0.19548      0.64672      0.60514  
#>       sppDF      minedno  
#>     0.04602      1.01447  
#> 
#> Zero-inflation model:
#> (Intercept)        sppPR        sppDM      sppEC-A      sppEC-L     sppDES-L  
#>      1.7556       1.6785      -0.4269       1.1046      -0.4269      -0.6716  
#>       sppDF      minedno  
#>     -0.4269      -2.4038
mean(ranef(m3)$cond$site[[1]])
#> [1] 0.005181326

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

[seananderson]

I think the difference is that the glmmTMB examples are means across many random effect values, similar to having a bunch of IID spatiotemporal fields and then taking the mean of all those. I'd expect those to converge towards zero. But, I think the equivalent here would be having 1 draw from a N(0, sigma) random intercept distribution, which could be much further from zero.

[Lewis-Barnett-NOAA]

Thanks for helping me think through this fellas.

Now that I've received an example from the collaborator, I don't think the singularity error is related to the spatially varying coefficients at all and I haven't been able to reproduce the error given this example. I'm thinking it has to be a package version issue, as this collaborator was able to fit the same models and data previously. I'll update this if I can sort out the source of the problem.