Mean and median bias reduction in generalized linear models

This paper presents an integrated framework for estimation and inference from generalized linear models using adjusted score equations that result in mean and median bias reduction. The framework unifies theoretical and methodological aspects of past research on mean bias reduction and accommodates, in a natural way, new advances on median bias reduction. General expressions for the adjusted score functions are derived in terms of quantities that are readily available in standard software for fitting generalized linear models. The resulting estimating equations are solved using a unifying quasi-Fisher scoring algorithm that is shown to be equivalent to iteratively reweighted least squares with appropriately adjusted working variates. Formal links between the iterations for mean and median bias reduction are established. Core model invariance properties are used to develop a novel mixed adjustment strategy when the estimation of a dispersion parameter is necessary. It is also shown how median bias reduction in multinomial logistic regression can be done using the equivalent Poisson log-linear model. The estimates coming out from mean and median bias reduction are found to overcome practical issues related to infinite estimates that can occur with positive probability in generalized linear models with multinomial or discrete responses, and can result in valid inferences even in the presence of a high-dimensional nuisance parameter.


Introduction
The flexibility of generalized linear models (McCullagh and Nelder 1989) in handling count, categorical, positive and real-valued responses under a common modelling framework has not only made them a typical choice in applications but also the focus of much methodological research on their estimation and use in inference.
Suppose that y 1 , . . . , y n are observations on independent random variables Y 1 , . . . , Y n , each with probability density or mass function of the exponential family form for some sufficiently smooth functions b(·), c 1 (·), a(·) and c 2 (·), and fixed observation weights m 1 , . . . , m n . The expected value and the variance of Y i are then E(Y i ) = μ i = b (θ i ) and var(Y i ) = φb (θ i )/m i = φV (μ i )/m i , respectively, where b (θ i ) and b (θ i ) are the first two derivatives of b(θ i ). Compared to the normal distribution, exponential family models are generally heteroscedastic because the response variance depends on the mean through the variance function V (μ i ), and the dispersion parameter φ allows shrinking or inflating that contribution of the mean. A generalized linear model (GLM) links the mean μ i to a linear predictor η i through a monotone, sufficiently smooth link function Estimated bias (B), root mean squared error (RMSE), percentage of underestimation (PU), mean absolute error (MAE) of maximum likelihood estimator and coverage of nominally 95% Wald-type confidence intervals (C), based on 10,000 samples under the ML fit. The summary B 2 /SD 2 is the relative increase in mean squared error from its absolute minimum due to bias. The results include the same summaries of the moment-based estimator of φ (row marked with * ). All reported figures are ×100 of their actual value and < 0.01 is used for a value that is less than 0.01 in absolute value g(μ i ) = η i with η i = p t=1 β t x it where x it is the (i, t)th component of a model matrix X , and β = (β 1 , . . . , β p ) . An intercept parameter is typically included in the linear predictor, in which case x i1 = 1 for all i ∈ {1, . . . , n}.
Estimation of the parameters of GLMs is commonly done using maximum likelihood (ML) because of the limiting guarantees that the ML estimator provides assuming that the model assumptions are adequate. Specifically, the ML estimator (β ,φ) is consistent, asymptotically unbiased and asymptotically efficient with a limiting normal distribution centred at the target parameter value and a variance-covariance matrix, given by the inverse of the Fisher information matrix, which is also the Cramér-Rao lower bound for the variance of unbiased estimators. These properties are used as reassurance that inferential procedures based on Wald, score or likelihood ratio statistics will perform well in large samples. Another reason that ML is the default estimation method for GLMs is that maximizing the likelihood can be conveniently performed by iteratively reweighted least squares (IWLS; Green 1984), requiring only standard algorithms for least squares and the evaluation of working weights and variates at each iteration.
Nevertheless, the properties of the ML estimator and of the associated inferential procedures that depend on its asymptotic normality may deteriorate for small or moderate sample sizes or, more generally, when the number of parameters is large relative to the number of observations. Example 1.1 To illustrate the differences between finite sample and limiting behaviour of the ML estimator and associate inferential procedures, consider the data in McCullagh and Nelder (1989, Sect. 8.4.2) of mean blood clotting times in seconds for nine percentage concentrations of normal plasma and two lots of clotting agent. The plasma concentrations are 5, 10, 15, 20, 30, 40, 60, 80, 100, with corresponding clotting times 118, 58, 42, 35, 27, 25, 21, 19, 18 for the first lot, and 69,35,26,21,18,16,13,12,12 for the second lot, respectively. We fit a gamma GLM with log μ i = 4 t=1 β t x it , where μ i is the expectation of the ith clotting time, x i1 = 1, x i2 is 1 for the second lot and 0 otherwise, x i3 is the corresponding (log) plasma concentration and x i4 = x i2 x i3 is an interaction term. The ML estimates areβ = (5.503, − 0.584, − 0.602, 0.034) andφ = 0.017. Table 1 shows the estimated bias, root mean squared error, percentage of underestimation and mean absolute error of the ML estimator from 10,000 simulated samples at the ML estimates, with covariates values fixed as in the original sample. The table also includes the same summaries of the moment-based estimator of φ (see, for example, McCullagh and Nelder 1989, Sect. 8.3, and the summary.glm function in R). The ML estimator of the regression parameters illustrates good bias properties, with distributions that have a mode around the parameter value used for simulation. On the other hand, the ML estimator of the dispersion parameter is subject to severe bias, which inflates the mean squared error by 54.13% from its absolute minimum, and has a severely right skewed distribution. Note here that the latter observation holds for any monotone transformation of the dispersion parameter. The moment-based estimator on the other hand has a much smaller bias, probability of underestimation closer to 0.5, and its use delivers a marked improvement to the coverage of standard confidence intervals for all model parameters.
Improvements in first-order inference based on ML can be achieved in several ways. For instance, bootstrap methods guarantee both correction of bias and higher-order accurate inference. Alternatively, analytical methods derived from higher-order asymptotic expansions based on the likeli-hood (see, for instance, Brazzale et al. 2007) have been found to result in accurate inference on model parameters. Nevertheless, bootstrap methods typically require intensive computation, and analytical methods, typically, require tedious, model-specific algebraic effort for their implementation. Furthermore, both bootstrap and analytical methods rely on the existence of the ML estimate, which is not always guaranteed. Such an example is GLMs with multinomial or discrete responses (Heinze and Schemper 2002;Kosmidis 2014b).
This paper presents a unified approach for mean and median bias reduction (BR) in GLMs using adjusted score functions (Firth 1993;Kosmidis and Firth 2009;and Kenne Pagui et al. 2017, respectively). Specifically, Firth (1993) and Kosmidis and Firth (2009) achieve higher-order BR of the ML estimator through the additive adjustment of the score equation. Kenne Pagui et al. (2017) use a similar approach in order to obtain component-wise higher-order median BR of the ML estimator, i.e. each component of the estimator has, to third order, the same probability of underestimating and overestimating the corresponding parameter component. We illustrate how those methods can be implemented without sacrificing the computational simplicity and the first-order inferential properties of the ML framework, and illustrate that they provide simple and practical solutions to the issue of boundary estimates in models with categorical responses.
Explicit, general formulae are derived for the adjusted score equations that produce higher-order mean and median unbiased estimators for GLMs. It is shown that, like ML, both mean and median BR can be conveniently performed by IWLS after the appropriate adjustment of the working variates for ML. Extensive empirical evidence illustrates that such an adjustment of IWLS leads to a stable estimation procedure even in case in which standard IWLS for ML estimation diverges.
Each method possesses invariance properties that can be more useful or less desirable depending on the GLM under consideration; the estimators resulting from mean BR (mean BR estimators, in short) are exactly invariant under linear transformations of the parameters in terms of the mean bias of the transformed estimators, which is useful, for example, when estimation and inference on arbitrary contrasts of the regression parameters is of interest. These invariance properties do not extend, though, to more general nonlinear transformations. On the other hand, median BR delivers estimators that are exactly invariant in terms of their improved median bias properties under general componentwise transformations of the parameters, which is useful, for example, when a dispersion parameter needs to be estimated from data. However, estimators from median BR are not invariant in terms of the median bias properties under more general transformations, for example, parameter contrasts. In order to combine the desirable invariance properties of each method when modelling with GLMs, we exploit the Fisher orthogonality (Cox and Reid 1987) of the mean and dispersion parameters to formally derive a novel mixed adjustment approach that delivers estimators of the regression parameters with improved mean bias and estimators for any unknown dispersion parameter with improved median bias.
Examples and simulation studies for various response distributions are used to demonstrate that both methods for BR are effective in achieving their respective goals and improve upon maximum likelihood, even in extreme settings characterized by high-dimensional nuisance parameters. Particular focus is given on special cases, like estimation of odds ratios from logistic regression models and estimation of log odds ratios from multinomial baseline category models.
All methods and algorithms discussed in this paper are implemented in the brglm2 R package (Kosmidis 2018), which has been used for all numerical computations and simulation experiments (see Supplementary Material).
The remaining of the paper is structured as follows. Section 2 gives a brief introduction to estimation using IWLS and shows how IWLS can be readily adjusted to perform mean or median BR. In particular, Sects. 2.1 and 2.2 review known results of ML estimation and explicit mean bias correction in generalized linear models. These subsections are useful to set up the notation and allow the introduction of mean and median bias-reducing adjusted score functions in Sects. 2.3 and 2.4, respectively. Inferential procedures based on the bias-reduced estimators are discussed in Sect. 3. Section 4 motivates the need for and introduces the mixed adjustment strategy for GLMs with a dispersion parameter. All methods are then assessed and compared through case studies and simulation experiments in Sects. 5 and 6. Section 6 also discusses how multinomial logistic regression models can be easily estimated with all methods using the equivalent Poisson log-linear model. Section 7 concludes the paper with a short discussion and possible extensions.

Iteratively reweighted least squares
The log-likelihood function for a GLM is n is the inverse of the link function. Suppressing the dependence of the various quantities on the model parameters and the data, the derivatives of the loglikelihood function with respect to the components of β and φ are Table 2 Working variates for ML and additional quantities needed in mean and median BR, for the most popular combinations of distributions and link functions c 1 (y i )} and ρ i = m i a i are the ith deviance residual and its expectation, respectively, with a i = a (− m i /φ), where a (u) = da(u)/du.
The ML estimatorsβ of β andφ of φ can be found by solution of the score equations s β = 0 p and s φ = 0, where 0 p is a p-dimensional vector of zeros. Wedderburn (1976) derives necessary and sufficient conditions for the existence and uniqueness of the ML estimator ofβ. Given that the dispersion parameter φ appears in the expression for s β in (1) only multiplicatively, the ML estimate of β can be computed without knowledge of the value of φ. This fact is exploited in popular software like the glm.fit function in R (R Core Team 2018). The jth iteration of IWLS updates the current iterate β ( j) for β by solving the weighted least squares problem where the superscript ( j) indicates evaluation at β ( j) , and z = (z 1 , . . . , z n ) is the vector of "working" variates with (Green 1984). Table 2 reports the working variates for well-used combinations of exponential family models and link functions. The updated β from the weighted least squares problem in (2) is equal to the updated β from the Fisher scoring step where i ββ is the (β, β) block of the expected information matrix about β and φ with a i = a (−m i /φ), where a (u) = d 2 a(u)/du 2 . Efron (1975) has shown that under the usual regularity conditions, the asymptotic mean bias of the ML estimatorγ for a general parametric model M γ can be reduced by the explicit correction ofγ asγ

Explicit mean bias reduction
is the first term in the expansion of the mean bias ofγ . Kosmidis (2014a) provides a review of explicit and implicit methods for mean BR. The general form of b γ is given in Cox and Snell (1968) in index notation and in Kosmidis and Firth (2010, Section 2) in matrix notation. For where ξ = (ξ 1 , . . . , ξ n ) T with ξ i = h i d i /(2d i w i ) and d i = d 2 μ i /dη 2 i , h i is the "hat" value for the ith observation, obtained as the ith diagonal element of the matrix H = X (X W X) −1 X W , and a i = a (− m i /φ), with a (u) = d 3 a(u)/du 3 . The derivation of b φ above is done using Kosmidis and Firth (2010, expressions (4.8) in Remark 3) to write b φ in terms of the first term in the expansion of the bias of 1/φ, which is given in Cordeiro and McCullagh (1991).
Note here that neither i φφ nor A * φ depend on β, and hence, the bias-reduced estimator for φ can be computed by knowledge ofφ only aŝ . Some algebra gives that the biasreduced estimator for β is whereB denotes evaluation of B at the ML estimator. Equivalently, and as also noted in Cordeiro and McCullagh (1991), the explicit correctionβ − b β (β,φ) can be performed by IWLS as in (2) up to convergence, and then making one extra step, where the working variate z is replaced by its adjusted version z + φξ . Table 2 gives the quantity φξ for some well-used GLMs. Firth (1993) shows that the solution of the adjusted score equations

Mean bias-reducing adjusted score functions
with A * β and A * φ as in (4) result in estimators β * and φ * with mean bias of smaller asymptotic order than the ML estimator.
A natural way to solve the adjusted score equations is through quasi-Fisher scoring (see Kosmidis and Firth 2010, for the corresponding quasi-Newton-Raphson iteration), where at the jth step the values for β and φ are updated as The term "quasi-" here reflects the fact that the expectation of the negative second derivatives of the scores, instead of the adjusted scores, is used for the calculation of the step size. Setting φ ( j+1) − φ ( j) = 0 in the above iteration shows that it has the required stationary point. Furthermore, if the starting values β (0) and φ (0) for iteration (7) are the ML estimates, then β (1) and φ (1) are the estimates from explicit BR, because s (0) β = 0 p and s (0) φ = 0. Figure 1 illustrates the quasi-Fisher scoring iterations for an one-parameter problem, starting from the ML estimate.
A similar calculation to that in Sect. 2.2 can be used to show that (7) can be written in terms of an IWLS step for β and an appropriate update for φ. In particular, Expression 8 makes apparent that, in contrast to ML, solving the mean bias-reducing adjusted score functions in GLMs with unknown dispersion parameter involves updating β and φ simultaneously. This is because b β generally depends on φ.
Despite that the stationary point of the iterative scheme (8) is the mean BR estimates, there is no theoretical guarantee for its convergence for general GLMs. However, substantial empirical studies have shown no evidence of divergence, even in cases in which standard IWLS (2) fails to converge. Some of those empirical studies are presented in Sects. 4, 5 and 6 of the present paper.

Median bias-reducing adjusted score functions
Kenne Pagui et al. (2017) introduce a family of adjusted score functions whose solution has smaller median bias than the ML estimator. Specifically, the solution γ † of s γ + A † γ = 0 is such that each of its components has probability 1/2 of underestimating the corresponding component of the parameter γ with an error of order O(n −3/2 ), as opposed to the error of order O(n −1/2 ) forγ . A useful property of the method is that it is invariant under component-wise monotone reparameterizations in terms of the improved median bias properties of the resulting estimators.
Some tedious but straightforward algebra starting from Kenne Pagui et al. (2017, expression (10)) gives that the median bias-reducing adjustments A † β and A † φ for GLMs have the form Fig. 1 Illustration of the quasi-Fisher scoring iterations for a model with a scalar parameter β, starting at the maximum likelihood estimateβ. One step gives the explicit mean reduced-bias estimator β − b β (β) of Sect. 2.2, and iterating until convergence results in the solution β * of the mean bias-reducing adjusted score equation where u = (u 1 , . . . , u p ) with In the above expressions, Similarly to the case of mean BR, the median biasreducing adjusted score equations can be solved using quasi-Fisher scoring or equivalently IWLS, where at the jth iteration Note here that the working variate for median BR is the one for mean BR plus the extra term φ Xu. Equivalently, and since the extra term is in the column space of X , the median BR IWLS update for β consists of a mean BR update for β as in (8), and a translation of the result by φu. Figure 2 illustrates that procedure. The core quantities in the definition of u are (10), and Table 2 includes their expressions for some well-used GLMs.
Similarly to (8), there is no theoretical guarantee for the convergence of the iterative scheme (11) for general GLMs. However, even in this case, our extensive empirical studies have produced no evidence of divergence.
Illustration of the IWLS update for computing the iterates of β for a given φ when performing mean BR and median BR . All quantities in the figure should be understood as being pre-multiplied by W 1/2 . The left figure shows the addition of φξ to the maximum likelihood working variates z and the subsequent projection onto C (the column space of W 1/2 X ) that gives the updated value for the mean BR estimates β * . The right figure illustrates the addition of φu on β * to give the updated value for the median BR estimates β †

Wald-type inference by plug-in
According to the results in Firth (1993) and Kenne Pagui et al. (2017), both θ * and θ † have the same asymptotic distribution as the ML estimator and hence are all asymptotically unbiased and efficient. Hence, the distribution of those estimators for finite samples can be approximated by a normal with mean θ and variance-covariance matrix {i(θ )} −1 , where i(θ ) is given in (3). The derivation of this result relies on the fact that both A * θ and A † θ are of order O(1) and hence dominated by the score function as information increases.
The implication of the above results is that standard errors for the components of θ * and θ † can be computed as for the ML estimator, using the square roots of the diagonal elements of {i(β * , φ * )} −1 and {i(β † , φ † )} −1 , respectively. As a result, first-order inference like standard Wald tests and Wald-type confidence intervals and regions are constructed in a plug-in fashion, by replacing the ML estimates with the mean BR or median BR estimates in the usual procedures in standard software.
Of course, for finite samples, Wald-type procedures based on the use of ML, mean and median bias reduction will yield different results. Such differences will disappear as the samples size increases. Sect. 3.2 explores those differences in normal linear regression models.

Normal linear regression models
Consider a normal regression model with y 1 , . . . , y n real- The adjustment terms A * β and A † β are zero for this model. As a result, the ML, mean BR and median BR estimators of β coincide with the least squares estimator On the other hand, the ML, mean BR and median BR estimators for The estimator φ * is mean unbiased for φ, and for this reason, it is the default choice for estimating the precision parameter in normal linear regression models. On the other hand, and as shown by Theorem 3.1, the use of φ † for Waldtype inference about β j based on asymptotic normality leads to inferences that are closer to the exact ones, based on the Student t n− p distribution, than when φ * is used, for all practically relevant values of n − p and α.
be the Wald-type confidence intervals for β j of nominal level 1 − α, based on the asymptotic normal distribution ofβ, β * and β † , respectively, where z α is the quantile of level α of the standard normal and where t n− p;α is the quantile of level α of the Student t distribution with n − p degrees of freedom, and define Len(I ) to be the length of interval I .
The proof of Theorem 3.1 is given in Appendix. Exact inferential solutions are not generally available for other GLMs with unknown dispersion parameter. It is therefore of interest to investigate whether the desirable behaviour of inference based on the median BR estimator, as demonstrated in Theorem 3.1 for the normal linear regression model, is preserved, at least approximately, in other models. Section 5.2 considers an example with gamma regression.

Mixed adjustments for dispersion models
In contrast to ML, mean BR is inherently not invariant to general transformations of the model parameters, in terms of its smaller asymptotic mean bias properties. This imposes a level of arbitrariness when carrying out inference on β in GLMs with unknown dispersion parameters, mainly because φ appears as a factor on the variance-covariance matrix {i(β, φ)} −1 of the estimators. For example, standard errors for β * will be different if the bias is reduced for φ or 1/φ. The mean BR estimates are exactly invariant under general affine transformations, which is useful in regressions that involve categorical covariates where invariance under parameter contrasts is, typically, required. On the other hand, median BR is invariant, in terms of smaller asymptotic median bias, under component-wise monotone transformations of the parameters, but it is not invariant under more general parameter transformations, like parameter contrasts.
In order to best exploit the invariance properties of each method, we propose the default use of a mixed adjustment that combines the mean bias-reducing adjusted score for β with the median bias-reducing adjusted score for φ by jointly solving with A * β and A † φ as in expressions (4) and (9), respectively. For GLMs with known φ, like Poisson or Binomial models, the mixed adjustment results in mean BR. On the contrary, for the normal linear models of Sect. 3.2 the mixed adjustment results in median BR because For general GLMs with unknown φ, the mixed adjustment provides the estimators β ‡ and φ ‡ , which are asymptotically equivalent to third order to β * and φ † , respectively. The proof of this result is a direct consequence of the orthogonality (Cox and Reid 1987) between β and φ and makes use of the expansions in Appendix of Kenne Pagui et al. (2017). Specifically, parameter orthogonality implies that terms up to order O(n −1 ) in the expansion of β ‡ − β are not affected by terms of order O(1) in s φ + A † φ . As a result, and up to order O(n −1 ), the expansion of β ‡ −β is the same as that of β * −β. The same reasoning applies if we switch the roles of β and φ, i.e. the expansion of φ ‡ − φ is the same to the expansion of φ † − φ, up to order O(n −1 ). Hence, β ‡ has the same mean bias properties as β * and φ ‡ has the same median bias properties as φ † . For this reason, we use the term mixed BR to refer to the solution of adjusted score functions resulting from the mixed adjustment.
In order to illustrate the stated invariance properties of the estimators coming from the mixed adjustment, we consider a gamma regression model with independent response random variables Y 1 , . . . , Y 12 , where, conditionally on covariates s i and t i , each Y i has a gamma distribution with mean μ i = exp(η i ) and variance φμ 2 i . The predictor η i is a function of regression parameters and the covariates, s i is a categorical covariate with values L1, L2 and L3, and t 1 , . . . , t 12 are generated from an exponential distribution with rate 1. Consider the three alternative parameterizations in Table 3. The identities β 1 = γ 1 , β 2 = γ 1 + γ 2 and β 3 = γ 1 + γ 3 follow directly.
We simulate 1000 independent response vectors from the parameter value (β 1 , β 2 , β 3 , β 4 , φ) = (− 1, − 0.5, 3, 0.2, 0.5) and estimate the three parameter vectors in Table 3 for each sample using the ML estimator, and the estimators resulting from the mean, median and mixed bias-reducing adjusted scores. The estimates for parameterizations I and III are used to estimate the probability P(|β 2 −γ 1 −γ 2 | > 1 ), and those for parameterizations I and II are used to estimate the probability P(|φ − exp(ζ )| > 2 ) for various values of Table 4 Probability P(|β 2 −γ 1 −γ 2 | > 1 ) for parameterizations I and III, and P(|φ − exp(ζ )| > 1 ) for parameterizations I and II for various values of 1 P(|β 2 −γ 1 −γ 2 | > 1 ) The ML estimator, the estimators from the mean, median and mixed bias-reducing adjusted scores are used in place of the tilded quantities. The figures are based on 1000 simulated response vectors from the gamma regression model of Table 3 with (β 1 , β 2 , β 3 , β 4 , φ) = (− 1, − 0.5, 3, 0.2, 0.5) not for mean BR. In contrast, the mixed adjustment strategy inherits the relevant properties of mean and median BR, and delivers estimators that are numerically invariant under linear contrasts of the mean regression parameters, and monotone transformations of the dispersion parameter. Section 5.2 further evaluates the use of the mixed adjustment in the estimation of gamma regression models.

Case studies and simulation experiments
In this section, we present results from case studies and confirmatory simulation studies that provide empirical support to the ability of mean and median BR to achieve their corresponding goals, i.e. mean and median bias reduction, respectively. In particular, in Sect. 5.2 we consider gamma regression, in which we also evaluate the mixed adjustment strategy of Sect. 4, while in Sect. 5.3 we consider logistic regression, showing how both mean and median BR provide a practical solution to the occurrence of infinite ML estimates. Finally, Sect. 5.4 evaluates the performance of mean and median BR in a logistic regression setting characterized by the presence of many nuisance parameters. In this case, ML estimation and inference are known to be unreliable, while both mean and median BR practically reproduce the behaviour of estimation and inference based on the conditional likelihood, which, in this particular case, is the gold standard.
All numerical computations are performed in R using the brglm2 R package (Kosmidis 2018). The brglm2 R package provides the brglmFit method for the glm R function that implements mean and median BR for any GLM using the quasi-Fisher scoring iteration introduced in Sect. 2.

Gamma regression model for blood clotting times
The regression model for the clotting data in Example 1.1 is fitted, here, using the mean, median and mixed bias-  Table 5. The estimates of regression parameters are practically the same for all methods. More marked differences between ML and the three adjusted score methods are noted in the estimates of the dispersion parameter. In particular, the estimates from the adjusted score methods result in notable inflation of the estimated standard errors for the regression parameters, with the median and mixed bias-reducing adjustments resulting in the largest inflation. In order to assess the quality of the estimates in Table 5, the simulated data sets in Example 1.1 are used to estimate the bias, the root mean squared error, the percentage of underestimation, and the mean absolute error of the various estimators, and the coverage of nominally 95% Wald-type confidence intervals. Table 6 reports the results. A comparison with the results of ML in Table 1 shows that the ML, mean BR, median BR and mixed BR estimators of β 1 , . . . , β 4 have similar bias and variance properties. On the other hand, the mean BR estimator of the dispersion parameter almost fully compensates for the mean bias of the ML estimator, while median BR and mixed BR give almost exactly 50% probability of underestimation. Furthermore, all BR methods deliver marked improvements in terms of empirical coverage

Logistic regression for infant birthweights
We consider a study of low birthweight using the data given in Hosmer and Lemeshow (2000, Table 2.1), which are also publicly available in the MASS R package. The focus here is on the 100 births for which the mother required no physician visits during the first trimester. The outcome of interest is a proxy of infant birthweight (1 if ≥ 2500g and 0 otherwise), whose expected value μ i is modelled in terms of explanatory variables using a logistic regression model with log{μ i /(1 − μ i )} = 7 t=1 β t x it , where x i1 = 1, x i2 and x i3 are the age and race (1 if white, 0 otherwise) of the mother, respectively, x i4 is the mother's smoking status during pregnancy (1 if yes, 0 if no), x i5 is a proxy of the history of premature labour (1 if any, 0 if none), x i6 is history of hypertension (1 if yes, 0 if no) and x i7 is the logarithm of the mother's weight at her last menstrual period. Table 7 gives the parameter estimates from ML, mean BR and median BR. Both mean BR and median BR deliver estimates that are shrunken versions of the corresponding ML estimates, with mean BR delivering the most shrinkage. This shrinkage translates to smaller estimated standard errors for the regression parameters. Kosmidis and Firth (2018) provide geometric insights for the shrinkage induced by mean BR in binary regression and prove that the mean BR estimates are always finite for full rank X .
The frequency properties of the resulting estimators are assessed by simulating 10,000 samples at the ML estimates in Table 7, with covariates fixed as in the observed sample, and re-estimating the model from each simulated sample. A total of 103 out of the 10,000 samples result in ML estimates with one or more infinite components due to data separation (Albert and Anderson 1984). The detection of infinite estimates was done prior to fitting the model using the linear programming algorithms in Konis (2007), as implemented in the detect_separation method of the brglm2 R package (Kosmidis 2018). The separated data sets were excluded when estimating the bias and coverage of Waldtype confidence intervals for the ML estimator. In contrast, the estimates from mean and median BR estimates were finite in all cases. For this reason, the corresponding summaries are based on all 10,000 samples. Table 8 shows the results. Both mean BR and median BR have excellent performance in terms of mean bias and probability of underestimation, respectively. Table 8 also includes summaries for the estimatorsψ t = eβ t , ψ * t = e β * t , ψ † t = e β † t of the odds ratios ψ t = e β t . Estimators of ψ t with improved bias properties have also been recently investigated in Lyles et al. (2012). The invariance properties of ML and median BR guarantee thatψ and ψ † are the ML and median BR All reported summaries, described in the caption of Table 1, for ML are conditional to the finiteness of the estimates. B ψ is the estimated bias in the ψ parameterization, and < 0.01 is used for a value that is less than 0.01 in absolute value estimators of ψ, respectively. As a result, ψ † t preserves its improved median bias properties. On the other hand, ψ * t is not, formally, the mean BR estimator of ψ. Nevertheless, it behaves best in terms of bias. The improved estimation and inference provided by mean and median BR become even more evident in more extreme modelling settings, as shown by the example in the next section.

Logistic regression for the link between sterility and abortion
We consider data from a retrospective, matched case-control study on the role of induced and spontaneous abortions in the aetiology of secondary sterility (Trichopoulos et al. 1976). The data are available in the infert data frame from the datasets R package. The two healthy control subjects from the same hospital were matched to each of 83 patients according to their age, parity and level of education. One of the cases could be matched with only one control; thus, there are a total of 248 records. Each record also provides the number of induced and spontaneous abortions, taking values 0, 1 and 2 or more. As is meaningful for retrospective case-control studies (see, for example, McCullagh and Nelder 1989, Sect. 4.3.3), we consider a logistic regression model with one fixed effect for each matched combination of cases and controls, and the number of induced and spontaneous abortions as the two categorical covariates of interest. In particular, the log odds of secondary sterility for the jth individual in the ith casecontrol combination are assumed to be where n i ∈ {2, 3}, x i j , x i j are indicator variables of 1 and 2 or more spontaneous abortions, respectively, and z i j and z i j are indicator variables of 1 and 2 or more induced abortions, respectively. The parameters λ 1 , . . . , λ 83 are the fixed effects for each matched combination of cases and controls, and the parameters of interest are β 1 , . . . , β 4 . Due to the many nuisance parameters, the maximum likelihood estimators of β 1 , . . . , β 4 are highly biased leading to misleading inference. A solution that is specific to logistic regression is to eliminate the fixed effects by conditioning on their sufficient statistics and maximize the conditional likelihood (CL). This can be done, for example, using the clogit function in the survival R package. As shown in Table 9, both mean and median BR give estimates that are close to the maximum CL estimates, practically removing all the bias from the ML estimates, and resulting also in a correction for the estimated standard errors.
This desirable behaviour of mean BR and median BR is in line with published theoretical results in stratified settings with nuisance parameters. In particular, Lunardon (2018) has recently shown that inferences based on mean BR in stratified settings with strata-specific nuisance parameters are valid under the same conditions for the validity of inference (Sartori 2003) based on modified profile likelihoods (see, for example, Barndorff-Nielsen 1983;Cox and Reid 1987;McCullagh and Tibshirani 1990;Severini 1998). The same equivalence is shown for median BR in Kenne Pagui et al. (2017).
The advantage of mean and median BR over maximum CL is their generality of application. As is shown in Table 2 mean and median BR can be used in models where a sufficient statistic does not exist, and hence, direct elimination of the nuisance parameters is not possible. One such example is probit regression, which is typically the default choice in many econometric applications stemming out from prospective studies. The further algorithmic simplicity for mean and median BR makes them also competitive to the various modified profile likelihoods.

The Poisson trick
Suppose that y 1 , . . . , y n are k-vectors of counts with k j=1 y i j = m i and that x 1 , . . . , x n are corresponding p-vectors of explanatory variables. The multinomial logistic regression model assumes that conditionally on x 1 , . . . , x n the vectors of counts y 1 , . . . , y n are realizations of independent multinomial vectors, with y i = (y i1 , . . . , y ik ), where the probabilities for the ith multinomial vector satisfy with k j=1 π i j = 1. Typically, x i1 = 1 for every i ∈ {1, . . . , n}. The above model is also known as the baseline category logit (see, for example, Agresti 2002, Sect. 7.1) because it uses one of the multinomial categories as a baseline for the definition of the log odds. Expression (13) has the kth category as baseline, but this is without loss of generality since any other log odds can be computed using simple contrasts of the parameter vectors γ 1 , . . . , γ k−1 .
Maximum likelihood estimation can be done either by direct maximization of the multinomial log-likelihood for (13) or using maximum likelihood for an equivalent Poisson log-linear model. Specifically, if y 11 , . . . , y nk are realizations of independent Poisson random variables with means μ 11 , . . . , μ nk , where then the score equations for λ i are m i = k j=1 μ i j , forcing the Poisson means to add up to the multinomial totals and the maximum likelihood estimates for γ 1 , . . . , γ k−1 to be exactly those that result from maximizing the multinomial likelihood for model (13) directly. Kosmidis and Firth (2011) proved that the equivalence of the multinomial logistic regression model (13) and the Poisson log-linear model (14) extends to the mean BR estimates of γ 1 , . . . , γ k−1 , if at each step of the iterative procedure for solving the adjusted score equations, the current values of the Poisson expectations μ i1 , . . . , μ ik are rescaled to sum up to the corresponding multinomial totals. Specifically, the results in Kosmidis and Firth (2011) suggest to prefix the IWLS update in (8) for the Poisson log-linear model (14) with the extra step . . . , n; s = 1, . . . , k) that rescales the Poisson means to sum to the multinomial totals. Then, W and the ML and mean BR quantities in the last row of Table 2 are computed usingμ The same argument applies the case of median BR. Given that the extra term in the IWLS update for median bias reduction in (11) depends on the parameters only through the response means, the same extra step of rescaling the Poisson means before the IWLS update of the parameters will result in an iteration that delivers the median BR estimates of the multinomial logistic regression model using the equivalent Poisson log-linear model.

Invariance properties
The mean BR estimator is invariant under general affine transformations of the parameters, and hence, direct contrasts result in mean BR estimators for any other baseline category for the response and any reference category in the covariates, without refitting the model. This is a particularly useful guarantee when modelling with baseline category models. In contrast, a direct transformation of the median BR estimates with baseline category k or a specific set of contrasts for the covariates is not guaranteed to result in median BR estimates for other baseline categories or contrasts in general.

Primary food choices of alligators
In order to investigate the extent that non-invariance impacts estimation and inference, we consider the data on food choice of alligators analysed in Agresti (2002, Sect. 7.1.2). The data come from a study of factors influencing the primary food choice of alligators. The observations are 219 alligators captured in four lakes in Florida. The nominal response variable is the primary food type, in volume, found in an alligator's stomach, which has five categories (fish, invertebrate, reptile, bird and other). The data set classifies the primary food choice according to the lake of capture (Hancock, Oklawaha, Trafford, George), gender (male and female) and size of the alligator (≤ 2.3 m long, > 2.3 m long).
Let s = 1 for alligator size > 2.3 metres and 0 otherwise, and let z H , z O , z T , z G be indicator variables for the lakes; for instance, z H = 1 for alligators on the lake Hancock and 0 otherwise. A possible model for the probabilities of food choice is where π ic is the probability for category c, with values corresponding to fish (c = 1), invertebrate (c = 2), reptile (c = 3), bird (c = 4) and other (c = 5). Model (15) is based on the choice of contrasts that would be selected by default in R. In order to investigate the effects of lack of invariance of median bias reduction, the set of contrasts used in Agresti (2002, Section 7.1.2) is considered where George is the reference lake and > 2.3 is the reference alligator size. These choices result in writing the food choice log odds as 3,4,5), (16) where s = 1 for alligator size ≤ 2.3 metres and 0 otherwise. The coefficients in the linear predictors of (15) and (16) are related as γ c1 = γ c1 +γ c2 +γ c3 , γ c2 = −γ c2 , γ c3 = γ c4 −γ c3 , γ c4 = γ c5 − γ c3 and γ c5 = −γ c3 . Table 10 gives the ML, mean BR and median BR estimates, along with the corresponding estimated standard errors of the coefficients of model (15). Table 10 shows also results of median BR γ , which correspond to the median BR estimates of γ transformed in the γ parameterization. As in logistic regression, the mean and median BR estimates are shrunken relative to the maximum likelihood ones with a corresponding shrinkage effect on the estimated standard errors.
The median BR and median BR γ estimates are almost the same, indicating that median BR, in this particular setting, is not affected by its lack of invariance under linear contrasts. The differences between the three methods are more notable when the observed counts are divided by two, as given in Table 11. In this case, data separation results in two of the ML estimates being infinite. This can generally happen with positive probability when data are sparse or when there are large covariate effects (Albert and Anderson 1984). As is the case for logistic regression (see Sect. 5.3), both mean and median BR deliver finite estimates for all parameters. The finiteness of the mean BR estimates has also been observed in Bull et al. (2002).
In order to better assess the properties of the estimators considered in Tables 10 and 11, we designed a simulation study where the multinomial totals for each covariate setting in the alligator food choice data set are progressively increased as a fraction of their observed values. Specifically, we consider the sets of multinomial totals {rm 1 , . . . , rm n } for r ∈ {0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 3, 4, 5}, where m i (i = 1, . . . , n) is the observed multinomial total for the ith combination of covariate values. For each value of r , we simulate 10,000 data sets from the ML fit of model (15) given in Table 10 and then compare the mean BR, median BR and median BR γ estimators in terms of relative bias and percentage of underestimation. The ML estimator is not considered in the comparison because the probability of infinite estimates is very high, ranging from 1.3% for r = 5 up to 76.4% for r = 0.5. In contrast, mean BR and median BR produced finite estimates for all data sets and r values considered. Figures 3 and 4 show the relative bias and the percentage of underestimation, respectively, for each parameter as a function of r . Overall, mean BR is preferable in terms of mean bias, while median BR achieves better median cen- tring for all the parameters. We note that even median BR γ has bias and probabilities of underestimation very close to those obtained directly under the γ parameterization. This confirms the indications from the observed data that, even if not granted by the theory, median BR is close to invariant under contrasts in the current model setting. As expected, the frequency properties of the three estimators converge to what we expect from standard ML asymptotics as r increases.
In particular, the bias converges to 0 and the percentage of underestimation to 50%.

Discussion
Fisher orthogonality (Cox and Reid 1987) of the mean and dispersion parameters dictates that the mixed approach to bias reduction is valid also for generalized linear models with dispersion covariates in Smyth (1989), and that estimation can be done by direct generalization of the IWLS iterations in (5) and (11), for mean and median bias reduction, respectively. Inference and model comparison has been based on Waldtype statistics. For special models, it is possible to form penalized likelihood ratio statistics based on the penalized  Table 10, for each r ∈ {0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 3, 4, 5}. The curves correspond to the mean BR (solid), median BR (dashed) and median BR γ (long dashed) estimators. The grey horizontal line is at zero log-likelihood that corresponds to the adjusted scores. A prominent example is logistic regression where the mean bias-reducing adjusted score is the gradient of the loglikelihood penalized by the logarithm of the Jeffreys' prior (see Heinze and Schemper 2002, where the profiles of the penalized log-likelihood are used for inference). In that case, the estimator from mean BR coincides with the mode of the posterior distribution obtained using the Jeffreys' prior (see also Ibrahim and Laud 1991). The same happens for Poisson log-linear models and for multinomial baseline category models. Even when a penalized log-likelihood corresponding to adjusted scores is not available (see Theorem 1 in Kosmidis and Firth 2009, for necessary and sufficient conditions for the existence of mean bias-reducing penalized likelihoods for generalized linear models), the adjustments to the score can, however, be seen as model-based penalties to the inferential quantities for maximum likelihood. In this sense, the adjustments introduce some implicit regularization to the estimation problem, which is just enough to achieve mean or median BR. In this framework, a general alternative to Wald-type statistics is score-type statistics with known asymptotic distributions, which can be readily defined as in Lindsay and Qu (2003). Let (β , φ) = (ψ , λ ) , with dim(ψ) = p 1 and dim(λ) = p − p 1 , i ψψ (ψ, λ) be a p 1 × p 1 matrix collecting the rows and columns of {i(ψ, λ)} −1 corresponding to ψ, and λ * ψ the estimator of λ resulting from the solution of the mean bias-reducing adjusted score equations on λ for fixed  Table 10, for each r ∈ {0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 3, 4, 5}. The curves corre-spond to the mean BR (solid), median BR (dashed) and median BR γ (long dashed) estimators. The grey horizontal line is at 50 ψ. Since the scores have an asymptotic normal distribution with mean zero and variance-covariance matrix i(ψ, λ) and the mean bias-reducing adjustment is of order O(1), has an asymptotic null χ 2 p 1 distribution. The same result holds for median BR, by replacing λ * ψ and A * ψ with λ † ψ and A † ψ . The adjusted score statistic can then be used for constructing confidence intervals and regions and testing hypotheses on any set of parameters of the generalized linear models, including constructing tables similar to analysis of deviance tables for maximum likelihood.
Finally, as is illustrated in the example of Sect. 5.4 and shown in Lunardon (2018) and Kenne Pagui et al. (2017), mean BR and median BR can be particularly effective for inference about a low-dimensional parameter of interest in the presence of high-dimensional nuisance parameters, while providing, at the same time, improved estimates of the nuisance parameters.

Supplementary material
The supplementary material includes R code and a report to fully reproduce all numerical results and figures in the paper.