LRT for comparing non-mixed to mixed models

[palday]

Closes #491.

This only allows one non-mixed model and it has to be the first argument. That's fairly restrictive, but generally speaking, I want to avoid testing a massive series of increasingly complicated fixed-effects models and then adding random-effects onto that series.

[palday]

@dmbates I can feel you looking over my shoulder and shaking your head when I write something like "the deviance isn't really the deviance". 🙂 Do you have a better way of handling / labelling this?

[dmbates]

The important point is that the difference in -2 log likelihood is the same as the difference in deviances would be if they could be calculated. Of course, now that I have written that I wonder if it makes sense when comparing a LinearModel and a LinearMixedModel. We know what the saturated model for a LinearModel is but we don't for a LinearMixedModel. My statement seems to indicate that the deviance for the saturated model in a LinearMixedModel would then have to be the same as that of a LinearModel but I am not sure that makes sense. I am starting to confuse myself.

[palday]

I'm glad I'm not the only one who struggles with some of this stuff. 😄

I think the relevant bit here is actually that the test statistic for LRT is -2 log(L_1/L2) which then is just the difference of -2 log likelihood. Using the true deviance was always just a convenience because the additive constants cancel out, when both models have the same constant.

All that said, my question is: should we change the label of that column and/or the internal fields? I would prefer to keep the internal fieldnames, but changing the show method's output is trivial.

[dmbates]

Hmm. For a linear mixed model it would be more accurate to label the column as -2 loglik. For a GLMM I think it should be deviance. I believe that for a Binomial GLMM that is not also Bernoulli (e.g. cbpp) the quantity in that column will not be -2 loglik

[palday]

Ah, yes, that is true. And indeed the loglikelihood and deviance for GLMM are defined independently.

Let me see what it would take to add a conditional to the show method --I think we'll need to add something to the LikelihoodRatioTest struct because right now there is no notion of which model it came from.

[codecov]

Codecov Report

Merging #508 (cdd1369) into main (715fac5) will increase coverage by 0.14%.
The diff coverage is 98.63%.

❗ Current head cdd1369 differs from pull request most recent head dae16e2. Consider uploading reports for the commit dae16e2 to get more accurate results

@@            Coverage Diff             @@
##             main     #508      +/-   ##
==========================================
+ Coverage   93.99%   94.14%   +0.14%     
==========================================
  Files          26       26              
  Lines        2164     2236      +72     
==========================================
+ Hits         2034     2105      +71     
- Misses        130      131       +1     

Impacted Files	Coverage Δ
src/linearmixedmodel.jl	95.47% <96.15%> (+0.02%)	⬆️
src/likelihoodratiotest.jl	98.47% <100.00%> (+0.79%)	⬆️
src/optsummary.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 715fac5...dae16e2. Read the comment docs.

[kleinschmidt]

I think I understand why true nesting isn't checked (it's hard to do in general) but it feels like a potential footgun...I suggested a few places where you could be clearer in the docstrings and/or add some additional heuristics for non-formula-bearing GLMs.

[dmbates]

I think I understand why true nesting isn't checked (it's hard to do in general) but it feels like a potential footgun...I suggested a few places where you could be clearer in the docstrings and/or add some additional heuristics for non-formula-bearing GLMs.

There is a way of checking nesting - approximately checking due to floating point, of course. The columns of the model matrix for the linear model should be in the column space of m.X where m is the linear mixed model. If you take a QR decomposition of m.X and multiply the linear model's model matrix by Q' then in each of the columns of the product, the squared length of the last n - rank(m.X) positions should be negligible compared to the squared length of the first rank(m.X) positions.

[palday]

There is a way of checking nesting - approximately checking due to floating point, of course. The columns of the model matrix for the linear model should be in the column space of m.X where m is the linear mixed model. If you take a QR decomposition of m.X and multiply the linear model's model matrix by Q' then in each of the columns of the product, the squared length of the last n - rank(m.X) positions should be negligible compared to the squared length of the first rank(m.X) positions.

I like this idea (and that's what the tolerance kwarg is for, I suppose), but I think that should be a different PR because we can add that to the nesting checks for two MixedModels.

[dmbates]

The CI failures on Julia v1.4 are because the newly added code calls contains which doesn't exist in v1.4.

Should we continue to test on 1.4? Is it because that is currently the LTS release?

[kleinschmidt]

I didn't look closely at the matrix checking stuff because it's a bit over my head but otherwise looks good.

[palday]

@dmbates The current LTS is 1.0! I'm happy to drop 1.4 and 1.5 with MixedModels 4.0, if for no other reason than CI speed. If/when 1.6 is declared LTS, my preference would be to have to LTS to change out tiered support architecture to "LTS" and "current release".