4141# ' standard approach of comparing differences of deviances to a Chi-squared
4242# ' distribution, a practice derived for Gaussian linear models or
4343# ' asymptotically, and which only applies to nested models in any case.
44- # '
45- # ' The values in `p_worse` column are computed using the normal
44+ # '
45+ # ' The values in the `p_worse` column are computed using the normal
4646# ' approximation and values from the columns `elpd_diff` and
4747# ' `se_diff`. Sivula et al. (2025) discuss the conditions when the
48- # ' normal approximation used for SE and `se_diff` is good. , and
48+ # ' normal approximation used for SE and `se_diff` is good, and the
4949# ' column `diag_pnorm` contains possible diagnostic messages: 1)
5050# ' small data (N < 100), 2) similar predictions (|elpd_diff| < 4),
5151# ' or 3) possible outliers (khat > 0.5).
5858# ' selection process. In that case users are recommended to avoid model
5959# ' selection based on LOO-CV, and instead to favor model averaging/stacking or
6060# ' projection predictive inference.
61- # '
61+ # '
6262# ' @seealso
6363# ' * The [FAQ page](https://mc-stan.org/loo/articles/online-only/faq.html) on
6464# ' the __loo__ website for answers to frequently asked questions.
@@ -122,24 +122,28 @@ loo_compare.default <- function(x, ...) {
122122 diffs <- mapply(FUN = elpd_diffs , loos [ord [1 ]], loos [ord ])
123123 elpd_diff <- apply(diffs , 2 , sum )
124124 se_diff <- apply(diffs , 2 , se_elpd_diff )
125+
125126 # compute probabilities that a model has worse elpd than the best model
126127 # using a normal approximation (Sivula et al., 2025)
127128 p_worse <- stats :: pnorm(0 , elpd_diff , se_diff )
128- p_worse [elpd_diff == 0 ] <- NA
129- N <- nrow( diffs )
129+ p_worse [elpd_diff == 0 ] <- NA
130+
130131 # diagnostics to assess whether the normal approximation can be trusted
131- if (N < 100 ) {
132+ N <- nrow(diffs )
133+ if (N < 100 ) {
132134 # small N (Sivula et al., 2025)
133135 diag_pnorm <- rep(" N < 100" , length(elpd_diff ))
134- diag_pnorm [elpd_diff == 0 ] = " "
136+ diag_pnorm [elpd_diff == 0 ] <- " "
135137 } else {
136138 diag_pnorm <- rep(" " , length(elpd_diff ))
137139 # similar predictions (Sivula et al., 2025)
138- diag_pnorm [elpd_diff > - 4 & elpd_diff != 0 ] <- " similar predictions"
139- # possible outliers in differences (Sivula et al., 2025;
140- # Vehtari et al., 2024)
140+ diag_pnorm [elpd_diff > - 4 & elpd_diff != 0 ] <- " similar predictions"
141+ # possible outliers in differences (Sivula et al., 2025; Vehtari et al., 2024)
141142 khat_diff <- rep(NA , length(elpd_diff ))
142- khat_diff [elpd_diff != 0 ] <- apply(diffs [,elpd_diff != 0 , drop = FALSE ], 2 , \(x ) ifelse(length(unique(x ))< = 20 , NA , posterior :: pareto_khat(x , tail = " both" )))
143+ khat_diff [elpd_diff != 0 ] <- apply(
144+ diffs [, elpd_diff != 0 , drop = FALSE ], 2 ,
145+ \(x ) ifelse(length(unique(x )) < = 20 , NA , posterior :: pareto_khat(x , tail = " both" )
146+ ))
143147 diag_pnorm [khat_diff > 0.5 ] <- paste0(" khat_diff > 0.5" )
144148 }
145149 rownames(comp ) <- rnms
@@ -159,29 +163,25 @@ loo_compare.default <- function(x, ...) {
159163# ' @export
160164# ' @param digits For the print method only, the number of digits to use when
161165# ' printing.
162- # ' @param simplify For the print method only, should only the essential columns
163- # ' of the summary matrix be printed? The entire matrix is always returned, but
164- # ' by default only the most important columns are printed.
165- # ' @param pnorm For the print method only, should we include the normal
166+ # ' @param p_worse For the print method only, should we include the normal
166167# ' approximation based probability of each model having worse performance than
167- # ' the best model?
168- print.compare.loo <- function (x , ... , digits = 1 , simplify = TRUE , p_worse = TRUE ) {
168+ # ' the best model? The default is `TRUE`.
169+ print.compare.loo <- function (x , ... , digits = 1 , p_worse = TRUE ) {
169170 xcopy <- x
170- if (inherits(xcopy , " old_compare.loo" )) {
171- if (NCOL(xcopy ) > = 2 && simplify ) {
172- patts <- " ^elpd_|^se_diff|^p_|^waic$|^looic$"
173- xcopy <- xcopy [, grepl(patts , colnames(xcopy ))]
174- }
175- } else if (NCOL(xcopy ) > = 2 && simplify ) {
176- xcopy <- xcopy [, c(" elpd_diff" , " se_diff" )]
171+ if (NCOL(xcopy ) > = 2 ) {
172+ xcopy <- xcopy [, c(" elpd_diff" , " se_diff" )]
177173 }
178174 if (p_worse ) {
179- print(cbind(.fr(xcopy , digits ), p_worse = .fr(x [," p_worse" ],2 ), diag_pnorm = x [, " diag_pnorm" ]), quote = FALSE )
180- invisible (x )
175+ print(
176+ cbind(.fr(xcopy , digits ),
177+ p_worse = .fr(x [, " p_worse" ], 2 ),
178+ diag_pnorm = x [, " diag_pnorm" ]),
179+ quote = FALSE
180+ )
181181 } else {
182182 print(cbind(.fr(xcopy , digits )), quote = FALSE )
183- invisible (x )
184183 }
184+ invisible (x )
185185}
186186
187187
0 commit comments