Error using LOO with: pointwise = TRUE

Foztarz commented 1 year ago

Thanks for this incredibly useful package! I'm running into issues trying to use loo() for model comparison. I have very large models (fitted on imputed data using brm_multiple()) so I'm trying to use pointwise estimates for LOO to avoid reaching maximum memory usage. loo(..., pointwise = FALSE) runs fine on (smaller) ordbetareg models, but with pointwise = TRUE I get the following error: Error in draws$dpars$mu[, i] It looks to me like the draws are not in the format that log_lik_ord_beta_reg() expects them to be in. Depending on how loo() extracts them, the specific draws for mu are probably stored in draws$dpars$mu$fe. I had a look at this in the example below (I'm also interested in effects on phi, which I think may have a similar problem). I really have no idea about the inner workings of LOO though.

    require(ordbetareg)
    # Test data ---------------------------------------------------------------
    set.seed(20230202)
    #data where both mu & phi change as a function of x
    dd <- data.frame(yy = round(plogis(
                            rnorm(n = 101, 
                                  mean = -50:50 / 5, 
                                  sd = seq(from = 5, 
                                           to = 0.1, 
                                           length.out = 101)
                                  )), 2),
                    xx = -50:50 / 5)
    # Fit ordered beta model --------------------------------------------------
    ordb_model<- ordbetareg( formula =
                  bf(yy ~ xx, phi ~ xx),
                  phi_reg = TRUE,
                  data = dd,
                  cores = parallel::detectCores()-1,
                  chains =  parallel::detectCores()-1,
                  iter = 1e3,
                  init = '0',
                  backend = 'cmdstanr'
                )

    # Inspect model predictions ----------------------------------------------
    loo_heavy <- loo(ordb_model)
    ## Computed from 3500 by 4901 log-likelihood matrix
    ## 
    ## Estimate     SE
    ## elpd_loo  -9065.8  524.2
    ## p_loo      8283.9  482.7
    ## looic     18131.7 1048.5
    ## ------
    ##   Monte Carlo SE of elpd_loo is NA.
    ## 
    ## Pareto k diagnostic values:
    ##   Count Pct.    Min. n_eff
    ## (-Inf, 0.5]   (good)     4047  82.6%   209       
    ## (0.5, 0.7]   (ok)        198   4.0%   57        
    ## (0.7, 1]   (bad)       183   3.7%   14        
    ## (1, Inf)   (very bad)  473   9.7%   1                

    loo_point <- loo(ordb_model, pointwise = TRUE)
    ##Error in draws$dpars$mu[, i] : incorrect number of dimensions
    loo_sub_point <- loo_subsample(ordb_model)
    ##Error in draws$dpars$mu[, i] : incorrect number of dimensions

    #indeed the draws do not have this shape
    sub_draws <- prepare_predictions(ordb_model, point_estimate = 'median')
    summary(sub_draws$dpars$mu)
    ## Length Class        Mode   
    ## family 18     customfamily list   
    ## ndraws  1     -none-       numeric
    ## nobs    1     -none-       numeric
    ## fe      2     -none-       list   
    ## sp      0     -none-       list   
    ## cs      0     -none-       list   
    ## sm      0     -none-       list   
    ## gp      0     -none-       list   
    ## re      0     -none-       list   
    ## ac      0     -none-       list  

    #estimates are stored in $fe$b
    head(sub_draws$dpars$mu$fe$b)

    log_lik_ord_beta_reg <-
      function(i, draws) {

        # mu <- draws$dpars$mu[,i]
        mu <- with( draws$dpars$mu$fe, b%*%X[i,]  ) #something like this to give estimate for that datapoint
        phi <- with( draws$dpars$phi$fe, b%*%X[i,] ) #something like this to give estimate for that datapoint
        mu <- plogis(mu) # expected on the response scale?
        phi <- plogis(phi) # expected on the response scale?
        y <- draws$data$Y[i]
        cutzero <- draws$dpars$cutzero
        cutone <- draws$dpars$cutone

        thresh1 <- cutzero
        thresh2 <- cutzero + exp(cutone)

        if(y==0) {
          out <- log(1 - plogis(qlogis(mu) - thresh1))
        } else if(y==1) {
          out <- log(plogis(qlogis(mu) - thresh2))
        } else {
          out <- log(plogis(qlogis(mu)-thresh1) - plogis(qlogis(mu) - thresh2)) + dbeta(y,mu*phi,(1-mu)*phi,log=T)
        }

        out

      }

    sum(log_lik_ord_beta_reg(20,sub_draws))
    ##-1.188826