Reporting of model comparison (from performance::compare_performance)

[DominiqueMakowski]

Reporting model comparisons indices textually has always been a pain in the a***, but I got an initial draft (only tested for LMs for now) that could be useful (at least it is for me, as I can copypasta the list of indices).

library(report)
#> report is in alpha - help us improve by reporting bugs on github.com/easystats/report/issues

m1 <- lm(Sepal.Length ~ Petal.Length * Species, data = iris)
m2 <- lm(Sepal.Length ~ Petal.Length + Species, data = iris)
m3 <- lm(Sepal.Length ~ Petal.Length, data = iris)

x <- performance::compare_performance(m1, m2, m3)
r <- report(x)
r
#> Regarding the explanatory power, m1 (R2 = 0.84, adj. R2 = 0.83, AIC = 106.77, BIC = 127.84, RMSE = 0.33, Sigma = 0.34) is not significantly different from m2 (R2 = 0.84, adj. R2 = 0.83, AIC = 106.23, BIC = 121.29, RMSE = 0.33, Sigma = 0.34, p = 0.186) and significantly different from m3 (R2 = 0.76, adj. R2 = 0.76, AIC = 160.04, BIC = 169.07, RMSE = 0.40, Sigma = 0.41, p < .001). Regarding the predictive power, there is strong evidence (BF = 26.52) in favour of the superiority of m2 compared to m1 and extreme evidence (BF < 1/1000) against the superiority of m3 compared to m1
summary(r)
#> Regarding the explanatory power, m1 (adj. R2 = 0.83, BIC = 127.84) is not significantly different from m2 (adj. R2 = 0.83, BIC = 121.29, p = 0.186) and significantly different from m3 (adj. R2 = 0.76, BIC = 169.07, p < .001). Regarding the predictive power, there is strong evidence (BF = 26.52) in favour of the superiority of m2 compared to m1 and extreme evidence (BF < 1/1000) against the superiority of m3 compared to m1

# Specific reports
report_statistics(x)
#> R2 = 0.84, adj. R2 = 0.83, AIC = 106.77, BIC = 127.84, RMSE = 0.33, Sigma = 0.34
#> R2 = 0.84, adj. R2 = 0.83, AIC = 106.23, BIC = 121.29, RMSE = 0.33, Sigma = 0.34
#> R2 = 0.76, adj. R2 = 0.76, AIC = 160.04, BIC = 169.07, RMSE = 0.40, Sigma = 0.41
report_performance(x)
#> m2 (R2 = 0.84, adj. R2 = 0.83, AIC = 106.23, BIC = 121.29, RMSE = 0.33, Sigma = 0.34) has a not significantly different explanatory power from m1 (R2 = 0.84, adj. R2 = 0.83, AIC = 106.77, BIC = 127.84, RMSE = 0.33, Sigma = 0.34, p = 0.186), and there is strong evidence (BF = 26.52) in favour of the hypothesis that m2 has a stronger predictive power than m1
#> m3 (R2 = 0.76, adj. R2 = 0.76, AIC = 160.04, BIC = 169.07, RMSE = 0.40, Sigma = 0.41) has a significantly different explanatory power from m1 (R2 = 0.84, adj. R2 = 0.83, AIC = 106.77, BIC = 127.84, RMSE = 0.33, Sigma = 0.34, p < .001), and there is extreme evidence (BF < 1/1000) against the hypothesis that m3 has a stronger predictive power than m1

^{Created on 2020-12-22 by the reprex package (v0.3.0)}

• Following this crossvalidated answer, I essentially decided to split the indices into two "groups", one which would cover "explanatory power" (most of which are absolute indices like R2, RMSE, RSS etc) and the other one "predictive power" (AIC, BIC etc.). This makes it rather nice because for these two groups we have a comparison statistics, the p-value from performance_lrt and the BIC-based BF.
• I will give a look at performance_lrt because it would be good to maybe retrieve (and store as attribute and pass it to compare_performance) the actual test statistics rather than just the p-value.
• Also, it'd be good in this case in which report_table() doesn't do anything (as it's already a table), to also have the printing of report tables to print the attributes like titles, notes etc.

[mattansb]

If you are splitting the explanatory from predictive power, shouldn't the AIC/BIC be in the predictive power part? (Also, unlike R2, R2adj is also a predictive measure 🤷‍♂️)

Also, I think you should add the F(df1,df2) to the p-value. (And maybe even cohens_f(model, model2) if you'd like...)

[DominiqueMakowski]

If you are splitting the explanatory from predictive power, shouldn't the AIC/BIC be in the predictive power part? (Also, unlike R2, R2adj is also a predictive measure 🤷‍♂️)

Well my rationale was that essentially you know you can roughly cluster indices at least in two groups that correlate more with one another, like BIC / AIC together and R2s RSS RMSE together, and these two groups have their related "comparison measure" (that compares the models to the first one), so it makes it a convenient breakdown. But additionally in the crossvalidated answer above, it's what they seem to say, basically R2s & co a measure of "explained variance" in the data whether the BIC / AIC are more related to the quality of the fit and its "generalizability" (?) in new data, hence predictive power.

And it's true that there are often cases where BIC/AIC rankings "contradicts" R2s/RMSE ones, so having a clear conceptual separation would facilitate reporting and provide the basis to strongly encourage people thinking about what these indices mean instead of just dumping out the whole table

Also, I think you should add the F(df1,df2) to the p-value. (And maybe even cohens_f(model, model2) if you'd like...)

Yeah, but these should probably be returned by performance::performance_lrt() first and then passed to compare_performance attributes.

[mattansb]

whether the BIC / AIC are more related to the quality of the fit and its "generalizability" (?) in new data, hence predictive power.

This is also what R2_adjusted does 🙃 - it is an estimate of R2 in the population (not in the data). That is, it generalizes to unobserved data.

[DominiqueMakowski]

Note to myself: drop the R2s from the report of model_comparison, because it can be misleading:

• Spiess, A.-N. & Neumeyer, N. An evaluation of R2 as an inadequate measure for nonlinear models in pharmacological and biochemical research: a Monte Carlo approach. BMC Pharmacol 10, 6–6 (2010).

And the R2s are reported anyway when the model is reported.

[DominiqueMakowski]

Or maybe not, it's still informative to eyeball the R2s

[DominiqueMakowski]

I asked this question related to performance_lrt(), could be useful:
https://stats.stackexchange.com/questions/502042/how-to-report-the-results-of-a-likelihood-ratio-test-lrt-from-the-anova-r-fu

[DominiqueMakowski]

library(report)
#> report is in alpha - help us improve by reporting bugs on github.com/easystats/report/issues

m1 <- lm(Sepal.Length ~ Petal.Length * Species, data = iris)
m2 <- lm(Sepal.Length ~ Petal.Length + Species, data = iris)
m3 <- lm(Sepal.Length ~ Petal.Length, data = iris)
report(performance::compare_performance(m1, m2, m3))
#> We compared three linear models; m1 (R2 = 0.84, adj. R2 = 0.83, AIC = 106.77, BIC = 127.84, RMSE = 0.33, Sigma = 0.34), m2 (R2 = 0.84, adj. R2 = 0.83, AIC = 106.23, BIC = 121.29, RMSE = 0.33, Sigma = 0.34) and m3 (R2 = 0.76, adj. R2 = 0.76, AIC = 160.04, BIC = 169.07, RMSE = 0.40, Sigma = 0.41).

^{Created on 2021-01-16 by the reprex package (v0.3.0)}

I simplified this now, as the tricky bits will go into report.test_performance()