summarylr: Summarizes a glm, adding a column of GLRT or Rao score P-values

Description

summarylr is an improved summary function for standard glm (stats package) adding LRT or Rao score P-values

Usage

summarylr(
  object,
  dispersion = NULL,
  correlation = FALSE,
  symbolic.cor = FALSE,
  ...,
  force = FALSE,
  debuglevel = level_warning,
  method = "LRT",
  keep.wald = FALSE
)

Arguments

object

glm object; as obtained by calling stats::glm or MASS::glm.nb.

dispersion

the dispersion parameter for the family used. Either a single numerical value or NULL (the default), when it is inferred from object (see stats::summary.glm).

correlation

logical; if TRUE, the correlation matrix of the estimated parameters is returned and printed (see stats::summary.glm).

symbolic.cor

logical; if TRUE, print the correlations in a symbolic form (see symnum) rather than as numbers (see stats::summary.glm).

...

Additional arguments to be passed to stats::summary.glm

force

logical; if TRUE, force computation of P-values in case of convergence problems.

debuglevel

integer value; set to 0 (default) to disable warnings, 1 to enable warnings and 2 to enable warnings and notes.

method

NULL or character vector of length 0, 1 or 2; may be code"LRT" or "Rao" or c("LRT", "Rao") to compute specified P-values. You can set method=NULL to compute no additional P-values.

keep.wald

logical; if TRUE, the standard Wald's P-values are kept in the output of print.summary.glmglrt. Even if keep.wald=FALSE, the standard wald P-values are not erased from the summary.glmglrt object. They are only hidden by print.summary.glmglrt.

Value

It returns a summary object of type summary.glmglrt that gets pretty printed by link[glmglrt:print.summary.glmglrt]{print.summary.glmglrt} The return value is an S3 object compatible with stats::summary.glm but with an additional field $extra field having sub-fields. $extra$pvalues is a numeric matrix with columns "LRT P-value" and/or "Rao P-value", containing the relevant P-values. As new columns may be added in future, you should rely on column names rather than column indexes. Only P-values of methods requested in the method parameter are stored in this matrix. $extra$debuglevel is equal to the debuglevel passed to summarylr. $extra$keep.wald is equal to the keep.wald passed to summarylr. In case of convergence problems, the field $extra$problem_of_convergence will be added. It will be a character string with the value "general" (because model$converged = FALSE), "all" (because all coefficients have huge variances) or "specific" (because at least one coefficient has a huge variance). Other problem strings may be added in the future. If weights are specified in a way that make P-values invalid, the field $extra$problem_weights will be added as a character string describing the problem. Actually, the only known problem is "non-constant".

Details

This function works the same as the standard summary.glm function but provides additionnal parameters The core parameter method="LRT" makes summarylr adds a column LRT P-value to the output. This P-value is computed by repeatdly fitting the model dropping one coefficient at a time and using the anova.glm(test="Chisq") function to perform generalized likelihood ratio test by approximation of the deviance difference to a chi-square distribution. This provides P-values less biased than the standard Wald P-values that summary provides. Moreover, this LRT method is consistent with the profile likelihood confidence intervals that confint.glm provides. The option method="Rao" generates Rao's score P-values. method="Chisq" is synonymous to method="LRT". For exhaustivity, the option method="Wald" (synonym "wald", "F") generates Wald's P-values. Several methods can be used, e.g. method=c("LRT","Rao") computes both LRT and Rao P-values. New methods may be added in the future.

Extra parameters are passed-through to the summary.glm function.

Examples

Run this code

# NOT RUN {
summarylr(glm(family="binomial", cbind(5,3)~1))
data(mtcars)
# do not properly converge (warnings)
mtcars$outcome = mtcars$disp > median(mtcars$disp)
mod=glm(family=binomial(log), data=mtcars,outcome ~ 0+qsec+wt,start=c(-0.1,0.3))
summarylr(mod) # warns that P-values are not computed because model did not converge
summarylr(mod, force=TRUE) # compute P-values anyway !
# also works with negative binomial models
summarylr(MASS::glm.nb(data=mtcars, I(cyl*gear) ~ 1+wt,link="sqrt"),test="LRT")
# }