VGAM (version 1.1-6)

lrt.stat: Likelihood Ratio Test Statistics Evaluated at the Null Values

Description

Generic function that computes likelihood ratio test (LRT) statistics evaluated at the null values (consequently they do not suffer from the Hauck-Donner effect).

Usage

lrt.stat(object, ...)
lrt.stat.vlm(object, values0 = 0, subset = NULL, omit1s = TRUE,
          all.out = FALSE, trace = FALSE, ...)

Arguments

object, values0, subset

Same as in wald.stat.vlm.

omit1s, all.out, trace

Same as in wald.stat.vlm.

…

Ignored for now.

Value

By default, a vector of signed square root of the LRT statistics; these are asymptotically standard normal under the null hypotheses. If all.out = TRUE then a list is returned with the following components: lrt.stat the signed LRT statistics, pvalues the 2-sided p-values, Lrt.stat2 the usual LRT statistic, values0 the null values.

Warning

See wald.stat.vlm.

Details

When summary() is applied to a vglm object a 4-column Wald table is produced. The corresponding p-values are generally viewed as inferior to those from a likelihood ratio test (LRT). For example, the Hauck and Donner (1977) effect (HDE) produces p-values that are biased upwards (see hdeff). Other reasons are that the Wald test is often less accurate (especially in small samples) and is not invariant to parameterization. By default, this function returns p-values based on the LRT by deleting one column at a time from the big VLM matrix and then restarting IRLS to obtain convergence (hopefully). Twice the difference between the log-likelihoods (or equivalently, the difference in the deviances if they are defined) is asymptotically chi-squared with 1 degree of freedom. One might expect the p-values from this function therefore to be more accurate and not suffer from the HDE. Thus this function is a recommended alternative (if it works) to summaryvglm for testing for the significance of a regression coefficient.

See Also

score.stat, wald.stat, summaryvglm, anova.vglm, vglm, lrtest, confintvglm, pchisq, profilevglm, hdeff.

Examples

# NOT RUN {
set.seed(1)
pneumo <- transform(pneumo, let = log(exposure.time),
                            x3 = rnorm(nrow(pneumo)))
fit <- vglm(cbind(normal, mild, severe) ~ let, propodds, pneumo)
cbind(coef(summary(fit)),
      "signed LRT stat" = lrt.stat(fit, omit1s = FALSE))
summary(fit, lrt0 = TRUE)  # Easy way to get it
# }