Add lrtest support methods

[palday]

As soon as JuliaStats/StatsModels.jl#162 lands, we should add in the following methods:

• [x] loglikelihood and ~deviance~ for REML-fitted models that throw an interpretable error.
• [x] isnested for GeneralizedLinearMixedModel
• [x] isnested for LinearMixedModel
• [x] refactor likelihoodratiotest to use isnested, perhaps add in deprecation (see previous discussion here)

Currently likelihoodratiotest can do one thing that lrtest cannot: model comparison on REML-fitted models. I'm wondering if we should keep it around for that purpose, give it some other name (suggestions?), or add in a custrom lrtest method that calls objective. There is one additional possibility: have deviance return -2 REML but issue a warning "This is not the deviance you are looking for." Given that some funding providers and standards-givers demand REML-fitted models, I don't think we can just omit support for comparing REML-fitted models.

There is one final option: I have code to do a bootstrap-based likelihood-ratio test. We could make our internal model comparison function do a bootstrapped comparison on the objective function and then it really would be different than the lrtest (and not depend on asymptotic assumptions about the deviance that the REML-criterion might not fulfill). In that case, what should we call it?

@dmbates Which of these is the smallest heresy?

[Nosferican]

Is not like the Wald test can't be used for REML models...

[kleinschmidt]

That StatsModels PR has landed now.

[palday]

I guess I know what I'm doing this weekend ;)

[dmbates]

Bump. @palday Feel free to assign this to me if you wish.

[palday]

One thing I don't like about the show method eventually called on the object returned by StatsModels.lrtest is that it doesn't describe the models. It is like

julia> lrtest(bigmodel, model, nullmodel)
  Likelihood-ratio test: 3 models fitted on 12 observations
  ──────────────────────────────────────────────
       DOF  ΔDOF  Deviance  ΔDeviance  p(>Chisq)
  ──────────────────────────────────────────────
  [1]    4         14.0571                      
  [2]    2    -2   15.9559     1.8988     0.3870
  [3]    1    -1   16.3006     0.3447     0.5571
  ──────────────────────────────────────────────

which is not quite as informative as I would like.

Originally posted by @dmbates in https://github.com/JuliaStats/MixedModels.jl/issues/366#issuecomment-691690702

I agree. This is what I'm thinking is a sensible course of action:

1. Add support for lrtest.
2. Refactor slightly so that our internal API takes advantage of the lrtest API / vice versa (e.g. isnested, etc.)
3. For now keep likelihoodratiotest in there, but add a note to its docstring that it may go away in the future.

Keeping likelihoodratiotest in there doesn't really cost us anything, especially since the safety checks are the same as those required by lrtest and all of the other computations (difference in likelihoods, etc.) are used elsewhere, so we're not even adding that much more code complexity. Right before the big 3.0 release, we can revisit this decision.

Also, while I'm on this: should I make isnested return false for REML models with different fixed-effects structures? I haven't yet checked how hard this would be to do programmatically (inspecting the formulae), but technically REML models aren't nested at that point (which of course is just one of many problems with the REML criterion).

I'll get this done in the next few days.

[nalimilan]

Yes, it would make sense to print the formula of each model (when it's not too long). But currently we don't have an API in StatsModels for that AFAIK. That's something we should add while getting rid of TableRegressionModel (see my old nl/tableregressionmodel branch). If that's useful for you right now, you could also just add the following stub in StatsBase or in StatsModels, and use it in a try... catch block in lrtest:

formula(x::StatisticalModel) =
    throw(ArgumentError("formula not implemented for $(nameof(typeof(x)))"))