brglm (version 0.5-2)

brglm: Bias reduction in Binomial-response GLMs

Description

Fits binomial-response GLMs using the bias-reduction method developed in Firth (1993) for the removal of the leading ($\mathop{\rm O}(n^{-1})$) term from the asymptotic expansion of the bias of the maximum likelihood estimator. Fitting is performed using pseudo-data representations, as described in Kosmidis (2007, Chapter 5). For estimation in binomial-response GLMs, the bias-reduction method is an improvement over traditional maximum likelihood because:

• the bias-reduced estimator is second-order unbiased and has smaller variance than the maximum likelihood estimator and
• the resultant estimates and their corresponding standard errors arealwaysfinite while the maximum likelihood estimates can be infinite (in situations where complete or quasi separation occurs).

Usage

brglm(formula, family = binomial, data, weights, subset, na.action,
      start = NULL, etastart, mustart, offset,
      control.glm = glm.control1(...), model = TRUE, method = "brglm.fit",
      pl = FALSE, x = FALSE, y = TRUE, contrasts = NULL,
      control.brglm = brglm.control(...), ...)
brglm.fit(x, y, weights = rep(1, nobs), start = NULL, etastart = NULL,
          mustart = NULL, offset = rep(0, nobs), family = binomial(),
          control = glm.control(), control.brglm = brglm.control(),
          intercept = TRUE, pl = FALSE)

Arguments

formula

as in glm.

family

as in glm. brglm currently supports only the "binomial" family with links "logit", "probit", "cloglog", "cauchit".

data

as in glm.

weights

as in glm.

subset

as in glm.

na.action

as in glm.

start

as in glm.

etastart

as in glm.

mustart

as in glm.

offset

as in glm.

control.glm

control.glm replaces the control argument in glm but essentially does the same job. It is a list of parameters to control glm.fit.

control

same as in glm. Only available to brglm.fit.

intercept

as in glm.

model

as in glm.

method

the method to be used for fitting the model. The default method is "brglm.fit", which uses either the modified-scores approach to estimation or maximum penalized likelihood (see the pl argument below). The standard

pl

a logical value indicating whether the model should be fitted using maximum penalized likelihood, where the penalization is done using Jeffreys invariant prior, or using the bias-reducing modified scores. It is only used when method = "brglm

x

as in glm.

y

as in glm.

contrasts

as in glm.

control.brglm

a list of parameters for controlling the fitting process when method = "brglm.fit". See documentation of brglm.control for details.

...

further arguments passed to or from other methods.

Value

• brglm returns an object of class "brglm". A "brglm" object inherits first from "glm" and then from "lm" and is a list containing the following components:
• coefficients
  as in glm.
• residuals
  as in glm.
• fitted.values
  as in glm.
• effects
  as in glm.
• R
  as in glm.
• rank
  as in glm.
• qr
  as in glm.
• family
  as in glm.
• linear.predictors
  as in glm.
• deviance
  as in glm.
• aic
  as in glm (see Details).
• null.deviance
  as in glm.
• iter
  as in glm.
• weights
  as in glm.
• prior.weights
  as in glm.
• df.residual
  as in glm.
• df.null
  as in glm.
• y
  as in glm.
• converged
  as in glm.
• boundary
  as in glm.
• ModifiedScores
  the vector of the modified scores for the parameters at the final iteration. If pl = TRUE they are the derivatives of the penalized likelihood at the final iteration.
• FisherInfo
  the Fisher information matrix evaluated at the resultant estimates. Only available when method = "brglm.fit".
• hats
  the diagonal elements of the hat matrix. Only available when method = "brglm.fit"
• nIter
  the number of iterations that were required until convergence. Only available when method = "brglm.fit".
• cur.model
  a list with components ar and at which contains the values of the additive modifications to the responses (y) and to the binomial totals (prior.weights) at the resultant estimates (see modifications for more information). Only available when method = "brglm.fit".
• model
  as in glm.
• call
  as in glm.
• formula
  as in glm.
• terms
  as in glm.
• data
  as in glm.
• offset
  as in glm.
• control.glm
  as control in the result of glm.
• control.brglm
  the control.brglm argument that was passed to brglm. Only available when method = "brglm.fit".
• method
  the method used for fitting the model.
• contrasts
  as in glm.
• xlevels
  as in glm.
• pl
  logical having the same value with the pl argument passed to brglm. Only available when method = "brglm.fit".

code

modifications

eqn

$n$

pkg

• brlr
• logistf

Warnings

1. It is not advised to use methods associated with model comparison (add1, drop1, anova, etc.) on objects of class "brglm". Model comparison when estimation is performed using the modified scores or the penalized likelihood is an on-going research topic and will be implemented as soon as it is concluded.

2. The use of Akaike's information criterion (AIC) for model selection when method = "brglm.fit" is controversial. AIC was developed under the assumptions that (i) estimation is by maximum likelihood and (ii) that estimation is carried out in a parametric family of distributions that contains the true model. At least the first assumption is not valid when using method = "brglm.fit". However, since the MLE is asymptotically unbiased, asymptotically the modified-scores approach is equivalent to maximum likelihood. A more appropriate information criterion seems to be Konishi's generalized information criterion (see Konishi & Kitagawa, 1996, Sections 3.2 and 3.3), which will be implemented in a future version.

Details

brglm.fit is the workhorse function for fitting the model using either the bias-reduction method or maximum penalized likelihood. If method = "glm.fit", usual maximum likelihood is used via glm.fit.

The main iteration of brglm.fit consists of the following steps:

1. Calculate the diagonal components of the hat matrix (seegethatsandhatvalues).
2. Obtain the pseudo data representation at the current value of the parameters (seemodificationsfor more information).
3. Fit a local GLM, usingglm.fiton the pseudo data.
4. Adjust the quadratic weights to agree with the original binomial totals.

Iteration is repeated until either the iteration limit has been reached or the sum of the absolute values of the modified scores is less than some specified positive constant (see the br.maxit and br.epsilon arguments in brglm.control).

The default value (FALSE) of pl, when method = "brglm.fit", results in estimates that are free of any $\mathop{\rm O}(n^{-1})$ terms in the asymptotic expansion of their bias. When pl = TRUE bias-reduction is again achieved but generally not at such order of magnitude. In the case of logistic regression the value of pl is irrelevant since maximum penalized likelihood and the modified-scores approach coincide for natural exponential families (see Firth, 1993).

For other language related details see the details section in glm.

References

Bull, S. B., Lewinger, J. B. and Lee, S. S. F. (2007). Confidence intervals for multinomial logistic regression in sparse data. Statistics in Medicine 26, 903--918. Firth, D. (1992) Bias reduction, the Jeffreys prior and {GLIM}. In Advances in GLIM and statistical modelling: Proceedings of the GLIM 92 conference, Munich, Eds. L.~Fahrmeir, B.~Francis, R.~Gilchrist and G.Tutz, pp. 91--100. New York: Springer.

Firth, D. (1992) Generalized linear models and Jeffreys priors: An iterative generalized least-squares approach. In Computational Statistics I, Eds. Y. Dodge and J. Whittaker. Heidelberg: Physica-Verlag.

Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika 80, 27--38. Heinze, G. and Schemper, M. (2002). A solution to the problem of separation in logistic regression. Statistics in Medicine 21, 2409--2419. Konishi, S. and Kitagawa, G. (1996). Generalised information criteria in model selection. Biometrika 83, 875--890. Kosmidis, I. (2007). Bias reduction in exponential family nonlinear models. PhD Thesis, Department of Statistics, University of Warwick.

Kosmidis, I. and Firth, D. (2008). Bias reduction in exponential family nonlinear models. (submitted)

See Also

glm, glm.fit

Examples

## Begin Example
data(lizards)
# Fit the GLM using maximum likelihood
lizards.glm <- brglm(cbind(grahami, opalinus) ~ height + diameter +
                  light + time, family = binomial(logit), data=lizards,
                  method = "glm.fit")
# Now the bias-reduced fit:
lizards.brglm <- brglm(cbind(grahami, opalinus) ~ height + diameter +
                  light + time, family = binomial(logit), data=lizards,
                  method = "brglm.fit")
lizards.glm
lizards.brglm
# Other links
update(lizards.brglm, family = binomial(probit))
update(lizards.brglm, family = binomial(cloglog))
update(lizards.brglm, family = binomial(cauchit))
# Using penalized maximum likelihood
update(lizards.brglm, family = binomial(probit), pl = TRUE)
update(lizards.brglm, family = binomial(cloglog), pl = TRUE)
update(lizards.brglm, family = binomial(cauchit), pl = TRUE)