"Error in if (varx == 0) { : missing value where TRUE/FALSE needed" with fit$loo() but not loo::loo(fit$draws("log_lik"))

[mhollanders]

Hi,

Unfortunately I don't have a reproducible example because this is just popping up with a few versions of a big model. When I fit the model with cmdstanr, the following works:

> loo::loo(fit_mvn$draws("log_lik"))

Computed from 400 by 590 log-likelihood matrix.

         Estimate     SE
elpd_loo -31300.8 3489.9
p_loo       323.5   14.9
looic     62601.5 6979.8
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.62]   (good)     367   62.2%   37      
   (0.62, 1]   (bad)      177   30.0%   <NA>    
    (1, Inf)   (very bad)  46    7.8%   <NA>    
See help('pareto-k-diagnostic') for details.

But this doesn't:

> fit_mvn$loo()
Error in if (varx == 0) { : missing value where TRUE/FALSE needed

I checked and there's no NAs anywhere:

> fit_mvn$draws("log_lik") |> anyNA()
[1] FALSE

Does anyone have any idea? Sorry I can't be more helpful with reproducible code. I'd be happy to share the Stan program with simulation file.

Thanks,

Matt

[jgabry]

Thanks for reporting this, that's strange. I'm glad at least loo::loo(fit_mvn$draws("log_lik")) works so this doesn't prevent you from using loo. I'm not sure where that error is coming from. I don't think loo or cmdstanr has a line of code containing if (varx == 0) so I'm guessing this is coming from another package that's being used internally? Do you by any chance have the traceback() available so we can see more detail about the source of the error message?

[mhollanders]

Hey, thanks for the help. traceback() gives me the following:

> fit$loo()
Error in if (varx == 0) { : missing value where TRUE/FALSE needed
> traceback()
8: posterior::autocovariance(sims[, i])
7: FUN(X[[i]], ...)
6: lapply(1:chains, FUN = function(i) posterior::autocovariance(sims[, 
       i]))
5: FUN(array(newX[, i], d.call, dn.call), ...)
4: apply(x, 3, ess_rfun)
3: relative_eff.array(exp(LLarray), cores = r_eff_cores)
2: loo::relative_eff(exp(LLarray), cores = r_eff_cores)
1: fit$loo()

Also, I think the same error is causing the following problem with the priorsense packages:

> priorsense::powerscale_plot_dens(fit, "psi_a")
Error in if (k < 1) { : missing value where TRUE/FALSE needed
In addition: Warning message:
Can't fit generalized Pareto distribution because all tail values are the same. 
> traceback()
12: ps_min_ss(k)
11: .pareto_smooth_extra_diags(k, S)
10: pareto_smooth.default(exp(as.numeric(log_ratios - max(log_ratios))), 
        r_eff = NULL, return_k = TRUE, extra_diags = TRUE, verbose = FALSE)
9: posterior::pareto_smooth(exp(as.numeric(log_ratios - max(log_ratios))), 
       r_eff = NULL, return_k = TRUE, extra_diags = TRUE, verbose = FALSE)
8: powerscale.priorsense_data(x = x, variable = variable, component = scaled_component, 
       alpha = alpha_seq[i], moment_match = moment_match, k_treshold = k_threshold, 
       resample = resample, transform = transform, prediction = prediction, 
       selection = likelihood_selection, ...)
7: powerscale(x = x, variable = variable, component = scaled_component, 
       alpha = alpha_seq[i], moment_match = moment_match, k_treshold = k_threshold, 
       resample = resample, transform = transform, prediction = prediction, 
       selection = likelihood_selection, ...)
6: powerscale_sequence.priorsense_data(psd, lower_alpha = lower_alpha, 
       upper_alpha = upper_alpha, length = length, variable = variable, 
       component = component, moment_match = moment_match, k_threshold = k_threshold, 
       resample = resample, transform = transform, prediction = prediction, 
       auto_alpha_range = auto_alpha_range, symmetric = symmetric, 
       prior_selection = prior_selection, likelihood_selection = likelihood_selection, 
       ...)
5: powerscale_sequence(psd, lower_alpha = lower_alpha, upper_alpha = upper_alpha, 
       length = length, variable = variable, component = component, 
       moment_match = moment_match, k_threshold = k_threshold, resample = resample, 
       transform = transform, prediction = prediction, auto_alpha_range = auto_alpha_range, 
       symmetric = symmetric, prior_selection = prior_selection, 
       likelihood_selection = likelihood_selection, ...)
4: powerscale_sequence.default(x, length = length, ...)
3: powerscale_sequence(x, length = length, ...)
2: powerscale_plot_dens.default(fit, "psi_a")
1: priorsense::powerscale_plot_dens(fit, "psi_a")

Is it likely to be an error I've made with the log_lik variable in the generated quantities?

[jgabry]

Ah, I think this might be an issue with calling loo::relative_eff (it might be returning Inf or NaN or something like that in this case). Do you still get an error if you do fit_mvn$loo(r_eff = FALSE)?

[jgabry]

Ah, I think this might be an issue with calling loo::relative_eff (it might be returning Inf or NaN or something like that in this case). Do you still get an error if you do fit_mvn$loo(r_eff = FALSE)?

@avehtari In cmdstanr's loo method we have:

r_eff <- loo::relative_eff(exp(LLarray), cores = r_eff_cores)

If the issue is that the log ratios are too small, should we change it to something like this?

r_eff <- loo::relative_eff(exp(LLarray + max(-LLarray)), cores = r_eff_cores)

[mhollanders]

Hey @jgabry, I don't get the error with r_eff = FALSE. Thanks so much!

[mhollanders]

FWIW, when I try to do stacking I also get the following error:

> loo_list <- list(loo(fit1$draws("log_lik")), loo(fit2$draws("log_lik")))
Warning messages:
1: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
 
2: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
 
> loo_model_weights(loo_list)
Error in optim(theta.old, fun, gradient, control = control, method = method,  : 
  initial value in 'vmmin' is not finite

loo_compare() doesn't give any errors.

[jgabry]

It seems to have something to do with the initial values used for the optimization step in the stacking algorithm (I guess you could try changing the optimization method that's used, there's an argument for that, but not sure it will make a difference). This might be hard for me to debug without an example I can play with. Or maybe @avehtari or @yao-yl have an idea.

[avehtari]

If the issue is that the log ratios are too small, should we change it to something like this?

Yes

It seems all the reported problems come from at least for one of the LOO folds, the importance ratios under- or overflow (fixed by subtracting the max) and the one importance ratio dominating (subtracting the max is not enough). It should be possible to add additional checks to avoid errors, but the resulting comparisons and model weights are still unreliable with that bad importance ratio distributions. @mhollanders can you share the log_lik draws?

[mhollanders]

Hi @avehtari, sorry for taking so long with this. I've attached the draws, hopefully a .csv is fine (I wasn't able to send the draws_matrix as .rds).
log_lik.csv

[avehtari]

@mhollanders a fix was already merged, can you try by installing loo from github?

[mhollanders]

Hey @avehtari, I installed loo from Github and also updated cmdstanr to the development version. I'm still getting the following:

> fit1$loo()
Error in if (varx == 0) { : missing value where TRUE/FALSE needed

and:

> loo1 <- loo(fit1$draws("log_lik"))
Warning message:
Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
 
> loo2 <- loo(fit2$draws("log_lik"))
Warning message:
Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
 
> loo::loo_model_weights(list(loo1, loo2))
Error in optim(theta.old, fun, gradient, control = control, method = method,  : 
  initial value in 'vmmin' is not finite

[avehtari]

@mhollanders is there any way you could share the fit1 object (Dropbox or something?) so that I could test fit1$loo()?

With the log_lik.csv you provided I was able to find the problem in stacking_weights(), and I have a fix that changes the error to warning, but that doesn't change the fact that you have way too many very high Pareto-k's and the model weights are thus not useful.

[mhollanders]

Hey @avehtari, sorry for taking so long to respond. I have a link here that you should be able to download it with.

Re: the high Pareto-k's, I'm finding these to be super high even where the model recovers the DGP input when there's site-level random effects. I didn't realise you couldn't still use the stacking to get the best predictive performance, irrespective of the Pareto-k.

[mhollanders]

Hey @avehtari, sorry to re-hash this, but it's still an issue where loo::loo() works fine but fit$loo() does not. I've re-attached a csv of posterior draws, as well as a figure of the site-level log_liks separated by region to visualise.

log_lik.csv