Fix ELPD

[seananderson]

I think ELPD is currently incorrect in sdmTMB. The averaging is happening across data points but the averaging is actually supposed to happen across MCMC samples, of which we have the equivalent of 1. So ELPD reduces to the sum of the log likelihoods, which is fold_elpd and sum_loglik in the sdmTMB_cv() output.

# from the 'loo' package:
pointwise_elpd_calcs <- function(ll){
  elpd <- matrixStats::colLogSumExps(ll) - log(nrow(ll))
  ic <- -2 * elpd
  cbind(elpd, ic)
}

# from the loo::elpd help:
# elpd(matrix): 
# > An S by N matrix, where S is the size of the posterior sample
# > (with all chains merged) and N is the number of data points.

# fake log likelihoods:
log_liks <- c(1, 2, 3, 4, 5, 6)

# a matrix of the log likelihoods, 1 row to represent 1 MCMC draw:
m <- matrix(log_liks, ncol = 6, nrow = 1)
m
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    1    2    3    4    5    6

# pointwise_elpd_calcs takes the mean log likelihood across MCMC samples:
p <- pointwise_elpd_calcs(m)
p
#>      elpd  ic
#> [1,]    1  -2
#> [2,]    2  -4
#> [3,]    3  -6
#> [4,]    4  -8
#> [5,]    5 -10
#> [6,]    6 -12

# which here, or with maximum likelihood, is the same

# then they are summed:
loo:::elpd_object(p, dim(p))
#> 
#> Computed from 6 by 2 log-likelihood matrix using the generic elpd function
#> 
#>      Estimate  SE
#> elpd     21.0 4.6
#> ic      -42.0 9.2
sum(p[,1])
#> [1] 21

# but sdmTMB is doing that first calculation across data points instead
# of samples:
sdmTMB:::log_sum_exp(log_liks) - log(length(log_liks))
#> [1] 4.664434
log(mean(exp(log_liks)))
#> [1] 4.664434

# so, the equivalent to ELPD is actually just
# `fold_elpd` and `sum_loglik`
sum(log_liks)
#> [1] 21

^{Created on 2023-06-23 with reprex v2.0.2}

That's not to say that what we're doing is completely wrong. I believe it should scale the same as the real ELPD and should render higher values for the 'better' model, but the scale is off. It's currently a log average likelihood. It's supposed to be a sum of log likelihoods.

[ericward-noaa]

I'm glad you caught this - I was having some issues recently trying to interpret the difference in ELPD values between models w/simulated data and it didn't make sense.

Is it worth just removing the ELPD calculation entirely? Currently the sum of log likelihoods is being returned, so the ELPD calculation is redundant

[seananderson]

Yes, it might break some people's code, but perhaps that's a good thing because it will make them stop, read the help file, and understand what happened rather than silently having the results for the elpd element change.