Description Usage Arguments Details Value Warnings Note Author(s) References See Also Examples
Fits binomialresponse GLMs using the biasreduction method developed in Firth (1993) for the removal of the leading (O(n^{1})) term from the asymptotic expansion of the bias of the maximum likelihood estimator. Fitting is performed using pseudodata representations, as described in Kosmidis (2007, Chapter 5). For estimation in binomialresponse GLMs, the biasreduction method is an improvement over traditional maximum likelihood because:
the biasreduced estimator is secondorder unbiased and has smaller variance than the maximum likelihood estimator and
the resultant estimates and their corresponding standard errors are always finite while the maximum likelihood estimates can be infinite (in situations where complete or quasi separation occurs).
1 2 3 4 5 6 7 8 9 10  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)

formula 
as in 
family 
as in 
data 
as in 
weights 
as in 
subset 
as in 
na.action 
as in 
start 
as in 
etastart 
as in 
mustart 
as in 
offset 
as in 
control.glm 

control 
same as in 
intercept 
as in 
model 
as in 
method 
the method to be used for fitting the model. The default
method is 
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
biasreducing modified scores. It is only used when

x 
as in 
y 
as in 
contrasts 
as in 
control.brglm 
a list of parameters for controlling the fitting
process when 
... 
further arguments passed to or from other methods. 
brglm.fit
is the workhorse function for fitting the model using
either the biasreduction 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:
Calculate the diagonal components of the hat matrix (see
gethats
and hatvalues
).
Obtain the pseudodata representation at the current value of the
parameters (see modifications
for more information).
Fit a local GLM, using glm.fit
on the pseudo data.
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 O(n^{1}) terms in the asymptotic expansion of their bias. When
pl = TRUE
biasreduction 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
modifiedscores approach coincide for natural exponential families (see
Firth, 1993).
For other language related details see the details section in
glm
.
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 
residuals 
as in 
fitted.values 
as in 
effects 
as in 
R 
as in 
rank 
as in 
qr 
as in 
family 
as in 
linear.predictors 
as in 
deviance 
as in 
aic 
as in 
null.deviance 
as in 
iter 
as in 
weights 
as in 
prior.weights 
as in 
df.residual 
as in 
df.null 
as in 
y 
as in 
converged 
as in 
boundary 
as in 
ModifiedScores 
the vector of the modified scores for the
parameters at the final iteration. If 
FisherInfo 
the Fisher information matrix evaluated at the
resultant estimates. Only available when 
hats 
the diagonal elements of the hat matrix. Only available
when 
nIter 
the number of iterations that were required until
convergence. Only available when 
cur.model 
a list with components 
model 
as in 
call 
as in 
formula 
as in 
terms 
as in 
data 
as in 
offset 
as in 
control.glm 
as 
control.brglm 
the 
method 
the method used for fitting the model. 
contrasts 
as in 
xlevels 
as in 
pl 
logical having the same value with the 
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 ongoing
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 modifiedscores
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.
1. Supported methods for objects of class "brglm"
are:
print
through print.brglm
.
summary
through summary.brglm
.
coefficients
inherited from the
"glm"
class.
vcov
inherited from the"glm"
class.
predict
inherited from the"glm"
class.
residuals
inherited from the"glm"
class.
and other methods that apply to objects of class
"glm"
2. A similar implementation of the biasreduction method could be done for every GLM, following Kosmidis (2007) (see also Kosmidis and Firth, 2009). The full set of families and links will be available in a future version. However, biasreduction is not generally beneficial as it is in the binomial family and it could cause inflation of the variance (see Firth, 1993).
3. Basically, the differences between maximum likelihood, maximum penalized likelihood and the modified scores approach are more apparent in small sample sizes, in sparse data sets and in cases where complete or quasicomplete separation occurs. Asymptotically (as n goes to infinity), the three different approaches are equivalent to first order.
4. When an offset is not present in the model, the modifiedscores based estimates are usually smaller in magnitude than the corresponding maximum likelihood estimates, shrinking towards the origin of the scale imposed by the link function. Thus, the corresponding estimated asymptotic standard errors are also smaller.
The same is true for the maximum penalized likelihood estimates when for example, the logit (where the maximum penalized likelihood and modifiedscores approaches coincide) or the probit links are used. However, generally the maximum penalized likelihood estimates do not shrink towards the origin. In terms of meanvalue parameterization, in the case of maximum penalized likelihood the fitted probabilities would shrink towards the point where the Jeffreys prior is maximized or equivalently where the quadratic weights are simultaneously maximized (see Kosmidis, 2007).
5. Implementations of the biasreduction method for logistic
regressions can also be found in the logistf package. In
addition to the obvious advantage of brglm
in the range of
link functions that can be used ("logit"
, "probit"
,
"cloglog"
and "cauchit"
), brglm
is also more
efficient computationally. Furthermore, for any userspecified link
function (see the Example section of family
), the
user can specify the corresponding pseudodata representation to be
used within brglm
(see modifications
for
details).
Ioannis Kosmidis, [email protected]
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 leastsquares approach. In Computational Statistics I, Eds. Y. Dodge and J. Whittaker. Heidelberg: PhysicaVerlag.
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. (2009). Bias reduction in exponential family nonlinear models. Biometrika 96, 793–804.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20  ## 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 biasreduced 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)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.