Skip to contents

Likelihood-ratio comparison of drmTMB fits

Source: R/methods.R
anova.drmTMB.Rd

anova(m0, m1, ...) tests nested drmTMB() fits with a likelihood-ratio test. The fits are ordered by their number of parameters and each model is tested against the one before it, as in lme4 and glmmTMB. The statistic is \(2 (\ell_{\mathrm{full}} - \ell_{\mathrm{reduced}})\) with df(full) - df(reduced) degrees of freedom, and the table has columns npar, logLik, deviance, Chisq, Df and Pr(>Chisq).

Usage

# S3 method for class 'drmTMB'
anova(object, ..., test = NULL)

Arguments

object

A drmTMB fit.

...

Further drmTMB fits to compare with object.

test

NULL (the default), "Chisq" or "LRT"; only the likelihood-ratio test is offered.

Value

An object of class c("anova", "data.frame") with one row per fit, in order of increasing number of parameters.

Details

Guards, all of which abort with a message saying what to do next:

  • the fits must be on the same observations;

  • ML and REML fits are never mixed, and REML fits must share their fixed-effect (mean) structure, because restricted likelihoods are comparable only across variance structures;

  • penalized (MAP) fits abort, as their variance components are shrunk by the prior;

  • experimental MSPL fits abort, as for every likelihood method;

  • engine = "julia" fits are not supported, and neither are models with the same number of parameters (compare those with AIC()).

Boundary tests. When the larger model adds a random-effect standard deviation, the null value 0 sits on the edge of the parameter space and the usual chi-square reference is conservative. A drmTMB_anova_boundary_warning is signalled for such a pair. If exactly one or two standard deviations are added and nothing else, Pr(>Chisq) is the chi-bar-square mixture from lrt_boundary(); the attribute "p.method" records "chibar" or "chisq" per row. Otherwise the naive chi-square p-value is reported and the warning says it is conservative.

A negative statistic (the smaller model fitted better) is reported as is and gives a p-value of 1 (or the boundary value); it signals non-nested models or a fit that did not converge.

See also

lrt_boundary(), AIC(), model-comparison

Examples

# \donttest{
set.seed(1)
n <- 80
dat <- data.frame(x = rnorm(n))
dat$y <- 0.2 + 0.6 * dat$x + rnorm(n)
m0 <- drmTMB(bf(y ~ 1, sigma ~ 1), family = gaussian(), data = dat)
m1 <- drmTMB(bf(y ~ x, sigma ~ 1), family = gaussian(), data = dat)
anova(m0, m1)
#> Likelihood-ratio comparison of 2 drmTMB fits (ML)
#> Each model is tested against the one before it, ordered by number of parameters.
#>    npar  logLik deviance  Chisq Df Pr(>Chisq)   
#> m0    2 -106.71   213.42                        
#> m1    3 -103.13   206.25 7.1648  1   0.007435 **
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# }