Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Deviances to LRT output #806

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Add Deviances to LRT output #806

wants to merge 4 commits into from

Conversation

jmgirard
Copy link
Contributor

@jmgirard jmgirard commented Mar 26, 2025

This helps make it clearer which model is best and aids in teaching. Could put the Log Likelihood in there instead, but it is the deviances that the LRT is comparing. Could also move the column left of Chi2 to emphasize that the Chi2 are the differences in deviances.

m1 <- suppressMessages(lme4::lmer(Sepal.Length ~ Petal.Width + (1 | Species), data = iris))
m2 <- suppressMessages(lme4::lmer(Sepal.Length ~ Petal.Width + Petal.Length + (1 | Species), data = iris))
m3 <- suppressMessages(lme4::lmer(Sepal.Length ~ Petal.Width * Petal.Length + (1 | Species), data = iris))
performance::test_likelihoodratio(m1, m2, m3)
#> # Likelihood-Ratio-Test (LRT) for Model Comparison (ML-estimator)
#> 
#> Name |   Model |    Dev | df | df_diff |  Chi2 |      p
#> -------------------------------------------------------
#> m1   | lmerMod | 202.22 |  4 |         |       |       
#> m2   | lmerMod | 116.96 |  5 |       1 | 85.26 | < .001
#> m3   | lmerMod | 116.12 |  6 |       1 |  0.84 |  0.359
#> 
#> Model 'm3' seems to have the best model fit.

p <- stats::pchisq(chi2, abs(dfs_diff), lower.tail = FALSE)

out <- data.frame(
Dev = devs,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd call it Deviance no?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure

@DominiqueMakowski
Copy link
Member

While we're at it, could we add a couple of sentences to the function description about what it does, what it's comparing etc. (because its documentation is pretty bare bones for now)

@jmgirard
Copy link
Contributor Author

jmgirard commented Mar 26, 2025

While we're at it, could we add a couple of sentences to the function description about what it does, what it's comparing etc. (because its documentation is pretty bare bones for now)

It is pulling in the docs from test_performance(), which is more elaborated. But I can add more if you think it's still too short. Or did you mean comments in the code?

@DominiqueMakowski
Copy link
Member

No I meant for the general docs. I forgot it took some from the general function, up to you then if there's any details you think are worth mentioning if not it's good to go

@jmgirard
Copy link
Contributor Author

I think we can push it then.

@strengejacke
Copy link
Member

hm, why are snapshots removed? I can add them back later.
Would you mind adding a test? Are there models where we are certain that no Dev column is returned, to that that case?

@jmgirard
Copy link
Contributor Author

hm, why are snapshots removed? I can add them back later. Would you mind adding a test? Are there models where we are certain that no Dev column is returned, to that that case?

Not sure. I cloned the main branch and just got those errors. I can write a few tests tonight or tomorrow.

@jmgirard
Copy link
Contributor Author

jmgirard commented Mar 27, 2025

I added deviance and tests for MLMs. We can probably do the same thing for OLS models but it's a bit tricky because "deviance" has a different meaning there. The chi-square is still based on -2*logLik but it isn't called a deviance in that context.

@bwiernik
Copy link
Contributor

It's not? I've always heard it called deviance in that context too

@jmgirard
Copy link
Contributor Author

jmgirard commented Mar 28, 2025

It's not? I've always heard it called deviance in that context too

I believe in OLS "deviance" is defined as the residual sum of squares.

# Fit model in OLS
fit <- lm(Sepal.Width ~ 1 + Species, data = iris)

# Get deviance in base R
deviance(fit)
#> [1] 16.962

# Get deviance in easystats
insight::get_deviance(fit)
#> [1] 16.962

# Get RSS for comparison
sum(resid(fit)^2)
#> [1] 16.962

# Get -2logLik in base R
-2 * logLik(fit)
#> 'log Lik.' 98.7326 (df=4)

# Get -2logLik in easystats
-2 * insight::get_loglikelihood(fit)
#> 'log Lik.' 98.7326 (df=4)

@jmgirard
Copy link
Contributor Author

Do we think it would be clearer to call this column (a) Deviance for GLM and MLM but -2LL or N2LL for OLS or (b) -2LL or N2LL for all models?

@DominiqueMakowski
Copy link
Member

How about something more generic like Difference

@jmgirard
Copy link
Contributor Author

jmgirard commented Mar 29, 2025

How about something more generic like Difference

Difference would work well for GLM where the deviance is a comparison to a saturated model, but that doesn't fit as well for OLS or MLM. I think I personally prefer -2LL with a definition in the documentation, but maybe Criterion could work if a generic term is desired? -2LL is sometimes called the log-likelihood criterion.

@jmgirard
Copy link
Contributor Author

jmgirard commented Mar 29, 2025

I also just found an edge-case where base R disagrees with insight: when fitting a fixed effects model in glmmTMB. Might be worth discussing if this is desired (or I could open a new issue on insight - just let me know):

fit <- glmmTMB::glmmTMB(Sepal.Length ~ 1 + Species, data = iris)

# Base R deviance = RSS
deviance(fit)
#> [1] 38.9562

# Insight deviance = -2LL
insight::get_deviance(fit)
#> [1] 223.452

logLik(fit)*-2
#> 'log Lik.' 223.452 (df=4)

@DominiqueMakowski
Copy link
Member

Why not "Criterion", but "N2LL" seems a bit too cryptic

@strengejacke
Copy link
Member

Thanks, look good! I'll review the next days and then we can merge it.

@jmgirard
Copy link
Contributor Author

jmgirard commented Apr 1, 2025

Thanks, look good! I'll review the next days and then we can merge it.

I'll be adding a little bit more soon. Maybe best to wait before final merge.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants