Why can't I replicate the lognormal model in sdmTMB with other methods?

[ericward-noaa]

Question that @kellijohnson-NOAA and I were asked this week: I'd like to use sdmTMB for non-spatial index standardization. I was comparing it to fitting a model with Stan (and combining the results myself). The estimates of year effects are off in sdmTMB -- why?

This example with the positive model is using the positive density:

pcod$log_density <- log(pcod$density)

fit_pos <- sdmTMB(
  density ~ 0 + as.factor(year), family = lognormal(),
  time = "year",
  data = pcod[pcod$density>0,],
  spatial = "off", spatiotemporal="off"
) 

And when I fit a model using rstanarm, my code looks like the chunk below

stan_pos <-stan_glm(
  log(density) ~ 0 + as.factor(year),
  family = gaussian,
  data = pcod[pcod$density>0,],
  prior = normal(location = 0, scale = 10),
  prior_intercept = normal(location = 0, scale = 10),
  iter = 2000
)

Note: fitting the model in brms returns equivalent estimates of the year effects. Similarly, using the lognormal() family in brms with an identity link returns the same. The coefficients from the stan_glm model and sdmTMB are below.

<style> </style>

year	stan	tmb_logn
2003.00	3.32	4.30
2004.00	3.54	4.52
2005.00	3.52	4.50
2007.00	2.87	3.85
2009.00	3.13	4.12
2011.00	3.71	4.70
2013.00	3.49	4.47
2015.00	3.57	4.55
2017.00	3.23	4.21

The sdmTMB estimates are perfectly correlated with the Bayesian ones, but off by a scale. So it doesn't affect the trend of the index -- but it does impact the absolute biomass estimates. Why is sdmTMB different?

[ericward-noaa]

The issue is that sdmTMB implements the widely used bias correction factor for the lognormal distribution, -sigma^ 2 / 2.0. This is in the code here:

sdmTMB/src/sdmTMB.cpp

Line 749 in be1ec3a

tmp_ll = sdmTMB::dlnorm(y_i(i,m), log(mu_i(i,m)) - pow(phi(m), Type(2)) / Type(2), phi(m), true);

The above approaches -- which either log-transform then fit a model with a Gaussian response, or use a lognormal family with an identity link, don't include the bias correction term.

You could fit the same model yourself, with the following Stan code

cat("data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
}
parameters {
  vector[K] b;  // population-level effects
  real<lower=0> sigma;  // dispersion parameter
}
model {
    vector[N] mu = X * b - sigma*sigma / 2.0; // penalty
    target += lognormal_lpdf(Y | mu, sigma);
}", file = "logn.stan")

And when we compare the estimates with the custom Stan code, they perfectly match the sdmTMB estimates:

<style> </style>

year	tmb_logn	stan_custom
2003	4.30	4.31
2004	4.52	4.54
2005	4.50	4.51
2007	3.85	3.86
2009	4.12	4.12
2011	4.70	4.71
2013	4.47	4.48
2015	4.55	4.56
2017	4.21	4.22

For very small values of sigma, these results will be more similar -- but for many cases, it's a good idea to include the bias correction term, as is in sdmTMB