---
abstract: |
    We develop a novel, general framework for the asymptotic reduction of the bias of $M$-estimators from unbiased estimating functions. The framework relies on additive, empirical adjustments to the estimating functions that depend only on the first two derivatives of the contributions to the estimating functions. The new estimation method has markedly broader applicability than previous bias-reduction methods by applying to models that are either partially-specified or that have a likelihood that is intractable or expensive to compute, and a surrogate objective is employed. The method also offers itself to easy, general implementations for arbitrary models by using automatic differentiation. This is in contrast to other popular bias-reduction methods that require either resampling or evaluation of expectations of products of log-likelihood derivatives. If $M$-estimation is by the maximization of an objective function, then, reduced-bias $M$-estimation can be achieved by maximizing an appropriately penalized objective. That penalized objective relates closely to information criteria based on the Kullback-Leibler divergence, establishing, for the first time, a strong link between reduction of estimation bias and model selection. The reduced-bias $M$-estimators are found to have the same asymptotic distribution, and, hence, the same asymptotic efficiency properties as the original $M$-estimators, and we discuss inference and model selection with reduced-bias $M$-estimates. The properties of reduced-bias $M$-estimation are illustrated in well-used, important modelling settings of varying complexity.\
    [Keywords: *automatic differentiation*; *composite likelihood*; *model selection*; *penalized likelihood*; *quasi likelihood*; *quasi Newton-Raphson*]{}
author:
- |
    Ioannis Kosmidis\
    `ioannis.kosmidis@warwick.ac.uk`\
    Department of Statistics, University of Warwick\
    Gibbet Hill Road, Coventry, CV4 7AL, UK\
    and\
    The Alan Turing Institute\
    British Library, 96 Euston Road, London, NW1 2DB, UK\
    and\
    Nicola Lunardon\
    `nicola.lunardon@unimib.it`\
    Department of Economics, Management and Statistics, University of Milano-Bicocca\
    Via Bicocca degli Arcimboldi 8, 20126 Milano, Italy
bibliography:
- 'robustbr.bib'
title: 'Empirical bias-reducing adjustments to estimating functions'
---

=1

Introduction
============

Reduction of estimation bias in statistical modelling is a task that has attracted immense research activity since the early days of the statistical literature. This ongoing activity resulted in an abundance of general bias-reduction methods of wide applicability. As is noted in @kosmidis:2014, the majority of those methods start from an estimator $\hat{{}\theta}$ and, directly or indirectly, attempt to produce an estimator $\tilde{{}\theta}$ of an unknown parameter ${{}\theta}$, which approximates the solution of the equation $$\label{bias_equation}
  \hat{{}\theta} - \tilde{{}\theta} = {{}B}_G(\bar{{}\theta}) \, ,$$ with respect to $\tilde{{}\theta}$. In the above equation, $G$ is the typically unknown, joint distribution function of the process that generated the data, ${{}B}_G({{}\theta}) = \expect_G(\hat{{}\theta} - {{}\theta})$ is the bias function, and $\bar{{}\theta}$ is the value that $\hat{{}\theta}$ is assumed to converge to in probability as information about ${{}\theta}$ increases, typically with the volume of the data. The need for approximating the solution of (\[bias\_equation\]) arises because, on one hand, $G$ is either unknown or the expectation with respect to $G$ is not available in closed form, and, on the other hand, the value of $\bar{{}\theta}$ is unknown.

Table \[br\_characteristics\] classifies prominent bias-reduction methods according to various criteria relating to their applicability and operation. Given the size of the literature on bias-reduction methods, we only cite key works that defined or greatly impacted the area.

Bias-reduction methods like the adjusted scores functions approach in @firth:1993, indirect inference in @gourieroux:1993, and iterated bootstrap in @kuk:1995 and @guerrier+dupuis-lozeron+etal:2019, assume that the model can be fully and correctly specified, in the sense that $G$ results from the assumed model for specific parameter values. That assumption allows to either have access to log-likelihood derivatives and expectations of products of those or to simulate from the model. In contrast, bias-reduction methods, like asymptotic bias correction [@efron:1975], bootstrap [@efron+tibshirani:1993; @hall+martin:1988], and jackknife [@quenouille:1956; @efron:1982] can also apply to at least partially-specified models. In this way, they provide the means for improving estimation in involved modelling settings, where researchers have, historically, resorted to surrogate inference functions in an attempt to either limit the number of hard-to-justify modelling assumptions or because the full likelihood function is impractical or cumbersome to compute or construct; see, for example, @wedderburn:1974 for quasi likelihoods, @liang+zeger:1986 for generalized estimating equations, and @lindsay:1988 and @varin+reid+firth:2011 for composite likelihood methods.

Another criterion for classifying bias-reduction methods comes from @kosmidis:2014. Therein, bias-reduction methods are classified according to whether they operate in an explicit manner by estimating ${{}B}_G(\bar{{}\theta})$ from the data and subtracting that from $\hat{{}\theta}$, or in an implicit manner by replacing ${{}B}_G(\bar{{}\theta})$ with $\hat{{{}B}}_G(\tilde{{}\theta})$ in (\[bias\_equation\]) for some estimator $\hat{{{}B}}_G$ of the bias function, and solving the resulting implicit equation.

Bias-reduction methods can also be classified according to whether the necessary approximation of the bias term in (\[bias\_equation\]) is performed analytically or through simulation. The vanilla implementations of asymptotic bias correction and adjusted score functions approximate ${{}B}_G({{}\theta})$ with a function ${{}b}({{}\theta})$ such that ${{}B}_G({{}\theta}) = {{}b}({{}\theta}) + O(n^{-3/2})$, where $n$ is a measure of how the information about ${{}\theta}$ accumulates. On the other hand, jackknife, bootstrap, iterated bootstrap, and indirect inference, generally, approximate the bias by simulating samples from the assumed model or an estimator of $G$, like the empirical distribution function. As a result, and depending on how demanding the computation of $\hat{{}\theta}$ is, simulation-based methods are typically more computationally intensive than analytical methods. Also, implicit, simulation-based methods require special care and ad-hoc considerations when approximating the solution of (\[bias\_equation\]), because the simulation-based estimator of ${{}B}_G({{}\theta})$ is not always differentiable with respect to ${{}\theta}$.

The requirement of differentiation of the log-likelihood or surrogate functions for some of the bias-reduction methods in Table \[br\_characteristics\] is also another area where considerable analytical effort has been devoted to (see, for example, @kosmidis+firth:2009 for multivariate generalized nonlinear models, and @grun+kosmidis+zeileis:2012 for Beta regression models). Nevertheless, differentiation is nowadays a task requiring increasingly less analytical effort because of the availability of comprehensive automatic differentiation routines [@griewank+walther:2008] in popular computing environments; such routines can be found, for example, in the FowardDiff Julia package [@revels+lubin+papamarkou:2016], and the CppAD package for C++ [@bell:2019] that enabled the development of software like the TMB package [@kristensen+nielsen+berg+etal:2016] for R [@rproject] which is a generic framework for fitting and inference from complex random effects models.

The vanilla versions of asymptotic bias correction and the adjusted score functions in @efron:1975 and @firth:1993, respectively, require the computation of expectations of products of log-likelihood derivatives under the model. Those expectations are intractable or expensive to compute for models with intractable or cumbersome likelihoods, and can be hard to derive even for relatively simple models [see, for example @grun+kosmidis+zeileis:2012 for the required expectations in Beta regression models].

  ----------------------------------------- --------------------- -------------------- ---------- ------------- ----------------- -------------------- --
                                                                                                                                                       
  Method                                     Model specification   Bias approximation     Type     Expectation   Differentiation   Original estimator  
  Asymptotic bias correction                       partial             analytical       Explicit       yes             yes                yes          
  Adjusted score functions                          full               analytical       Implicit       yes             yes                 no          
  Bootstrap                                        partial          simulation-based    Explicit       no              no                 yes          
  Jackknife                                        partial          simulation-based    Explicit       no              no                 yes          
  Iterated bootstrap / Indirect inference           full            simulation-based    Implicit       no              no                 yes          
  RB$M$-estimation                                 partial             analytical       Implicit       no              yes                 no          
  ----------------------------------------- --------------------- -------------------- ---------- ------------- ----------------- -------------------- --

\[br\_characteristics\]

Finally, except of the adjusted scores approach in @firth:1993, all the bias-reduction methods reviewed in Table \[br\_characteristics\] require the original estimator $\hat{{}\theta}$ and they cannot operate without it. For this reason, they directly inherit any of the instabilities that $\hat{{}\theta}$ may have. For example, in multinomial logistic regression, there is always a positive probability of data separation [@albert+anderson:1984] that results in infinite maximum likelihood estimates. Then, asymptotic bias correction, bootstrap, iterated bootstrap, and jackknife cannot be applied. The direct dependence on $\hat{{}\theta}$ may be more consequential for naive implementations of the latter three methods because they are simulation-based; even if data separation did not occur for the original sample, there is always positive probability that it will occur for at least one of the simulated samples. There is no easy way of knowing this before carrying out the simulation and such cases can only be handled in an ad-hoc way.

The current work develops a novel method for the reduction of the asymptotic bias of $M$-estimators from general, unbiased estimating functions. We call the new estimation method reduced-bias $M$-estimation, or RB$M$-estimation in short. Like the adjusted scores approach in @firth:1993, the new method relies on additive adjustments to the unbiased estimating functions that are bounded in probability, and results in estimators with bias of lower asymptotic order than the original $M$-estimators. The key difference is that the empirical adjustments introduced here depend only on the first two derivatives of the contributions to the estimating functions, and they require neither the computation of cumbersome expectations nor the, potentially expensive, calculation of $M$-estimates from simulated samples. Specifically, and as noted in the last row of Table \[br\_characteristics\], RB$M$-estimation i) applies to models that are at least partially-specified; ii) uses an analytical approximation to the bias function that relies only on derivatives of the contributions to the estimating functions; iii) does not depend on the original estimator; and iv) does not require the computation of any expectations. By relying only on derivatives of the contributions to the estimating functions, the new method is typically easier to implement for arbitrary models than other popular bias-reduction methods that require either resampling or the evaluation of moments of products of log-likelihood derivatives. In fact, automatic differentiation can be used to develop generic implementations of the new method whose only required input is an implementation of the contributions to the estimating functions. The GEEBRA Julia package (<https://github.com/ikosmidis/GEEBRA.jl>) by the Authors is a proof-of-concept of such an implementation.

If the estimating functions are the components of the gradient of an objective function, as is the case in maximum likelihood and maximum composite likelihood estimation, then, we show that bias reduction can always be achieved by the maximization of an appropriately penalized version of the objective. This is in contrast to the method in @firth:1993 which does not always have a penalized likelihood interpretation; see, for example, @kosmidis+firth:2009, who derive a necessary and sufficient condition for the existence of bias-reducing penalized likelihood for generalized linear models. Moreover, it is shown that the bias-reducing penalized objective closely relates to information criteria for model selection based on the Kullback-Leibler divergence. The functions of the parameters and the data that are used for bias reduction and model selection differ only by a known scalar constant. These observations, establish, for the first time, a strong link between reduction of bias in estimation and model selection. It is also shown that the RB$M$-estimators have the same asymptotic distribution, hence the same asymptotic efficiency properties, as the original $M$-estimators, and we discuss inference and model selection with RB$M$-estimates.

Section \[modelling\_setting\] introduces notation and sets up the general modelling setting we consider and the assumptions underpinning the theoretical developments. Section \[sec:theory\] and Section \[empirical\_adjustments\] are devoted on the derivation of the bias-reducing, empirical adjustments to the estimating functions and the specification of their components. Section \[sec:inference\] shows that the asymptotic distribution of the RB$M$ estimator is the same as that of the original $M$-estimator, and introduces Wald-type and generalized score approximate pivots that can be used for the construction of inferences. Section \[sec:penalties\] introduces and discusses bias-reducing penalized objectives and the links to model selection. A quasi Newton-Raphson iteration that has the BR$M$-estimates as its fixed point is introduced in Section \[sec:implementation\], along with discussion on its ingredients and on general implementations using automatic differentiation. The finite-sample properties of RB$M$-estimation are illustrated in well-used, important modelling settings of increasing complexity including, the estimation of the ratio of two means with minimal distributional assumptions (Section \[empirical\_adjustments\]), generalized linear models (Section \[sec:glms\]) and quasi likelihoods (Section \[sec:quasi\]), and composite likelihood methods for the estimation of gaussian max-stable processes (Section \[sec:maxstable\]), and for the estimation of multivariate probit models with correlated random effects for longitudinal ordinal responses (Section \[autoregressive\_probit\]). Section \[sec:discussion\] provides a discussion on the developments and possible extensions of the developments in this work, and introduces two alternative RB$M$-estimators.

Modelling setting and assumptions {#modelling_setting}
=================================

Estimating functions {#EF}
--------------------

Suppose that we observe the values ${{}y}_1, \ldots, {{}y}_k$ of a sequence of random vectors ${{}Y}_1, \ldots, {{}Y}_k$ with ${{}y}_i = (y_{i1}, \ldots, y_{ic_i})^\top \in \mathcal{Y} \subset
\Re^c$, possibly with a sequence of covariate vectors ${{}x}_1, \ldots, {{}x}_k$, with ${{}x}_i = (x_{i1}, \ldots, x_{iq_i})^\top \in \mathcal{X} \subset
\Re^{q_i}$. Any two distinct random vectors ${{}Y}_i$ and ${{}Y}_j$ may share random scalar components. We denote the distinct scalar random variables by $Z_1, \ldots, Z_m$, and by $G(z_1, \ldots, z_m)$ their typically unknown, underlying joint distribution function. In general $m \le \sum_{i = 1}^k c_i$, with equality only if ${{}Y}_1, \ldots, {{}Y}_k$ have distinct components. Let ${{}Y} = ({{}Y}_1^\top, \ldots, {{}Y}_k^\top)^\top$, and denote by $X$ the set of ${{}x}_1, \ldots, {{}x}_k$.

On of the common aims in statistical modelling is to estimate at least a sub-vector of an unknown parameter vector ${{}\theta} \in \Theta \subset \Re^p$ using data ${{}y}_1, \ldots, {{}y}_k$ and ${{}x}_1, \ldots, {{}x}_k$. This is most commonly achieved through a vector of $p$ estimating functions $\sum_{i = 1}^k{{}\psi}^i({{}\theta}) = (\sum_{i =
  1}^k\psi_1^i({{}\theta}), \ldots, \sum_{i = 1}^k\psi_{p}^i({{}\theta}))^\top$, where ${{}\psi}^i({{}\theta}) = {{}\psi}({{}\theta}, {{}Y}_i, {{}x}_i)$ and $\psi^i_r({{}\theta}) = \psi_r({{}\theta}, {{}Y}_i, {{}x}_i)$, $(r = 1,\ldots,p)$. In particular, ${{}\theta}$ is estimated by the $M$-estimator $\hat{{}\theta}$ [@vaart:1998 Chapter 5], which results by the solution of the system of estimating equations $$\label{estimating_eqautions}
  \sum_{i = 1}^k{{}\psi}^i({{}\theta}) = {{}0}_p \, ,$$ with respect to ${{}\theta}$, where ${{}0}_p$ is a $p$-vector of zeros. Prominent examples of estimation methods that fall within the above framework are estimation via quasi likelihoods [@wedderburn:1974] and generalized estimating equations [@liang+zeger:1986]. @stefanski+boos:2002 provide an accessible overview of estimating functions and demonstrate their generality and ability to tackle challenging estimation problems with fewer assumptions than likelihood-based estimation requires.

One way to derive estimating equations is through a, typically stronger, modelling assumption that ${{}Y}_i$ has a distribution function $F_i({{}y}_i | {{}x}_i, {{}\theta})$. The estimator $\hat{{}\theta}$ can then be taken to be the maximizer of the objective function $$\label{objective_function}
  l({{}\theta}) = \sum_{i =1}^k \log f_i({{}y}_i | {{}x}_i, {{}\theta}) \, ,$$ where $f_i({{}y}_i | {{}x}_i, {{}\theta})$ is the mixed joint density corresponding to $F_i({{}y}_i | {{}x}_i, {{}\theta})$. The word “mixed” is used here to allow for the fact that some of the components of ${{}y}_i$ may be continuous and some others may be discrete. If the objective function (\[objective\_function\]) is used, then the estimating functions in (\[estimating\_eqautions\]) have ${{}\psi}^i({{}\theta}) = \nabla \log f_i({{}y}_i | {{}x}_i,
{{}\theta})$, assuming that the gradient exists in $\Theta$. Prominent examples of estimation methods that involve an objective function of the form (\[objective\_function\]) are maximum likelihood and maximum composite likelihood [see, for example, @lindsay:1988; @varin+reid+firth:2011].

Assumptions
-----------

The assumptions we employ for the theoretical development in this work are listed below.

1.  \[consistency\] Consistency: $$\hat{{}\theta} \stackrel{p}{\longrightarrow} \bar{{{}\theta}} \, ,$$ where $\bar{{{}\theta}}$ is such that $\expect_G ({{}\psi}^i) = {{}0}_p$ for all $i \in \{1, \ldots, k\}$, with ${{}\psi}^i={{}\psi}^i(\bar{{}\theta})$ and $\expect_G(\cdot)$ denoting expectation with respect to the unknown joint distribution function $G$. In particular, we assume that $\hat{{}\theta} - \bar{{{}\theta}} = O_p(n^{-1/2})$, where $n \equiv n(k, m, q)$ is a measure of information about ${{}\theta}$.

2.  \[smoothness\] Smoothness: The derivatives of $\psi_r^i({{}\theta})$ exist up to the $5$-th order. In particular, $$l_{R_a}({{}\theta}) = \sum_{i=1}^k \frac{\partial^{a-1}
          \psi_{r_1}^i({{}\theta})}{\partial {{}\theta}^{r_2} \cdots \partial
          {{}\theta}^{r_a}} \, ,$$ exist for any set $R_a = \{r_1, \ldots, r_a\}$, with $r_j \in \{1, \ldots, p\}$ and $a \in \{1, \ldots, 5\}$, under the convention that $l_r({{}\theta}) = \sum_{i = 1}^k \psi_r^i({{}\theta})$ and that the components of ${}\theta$ are identified by superscripts.

3.  \[derivative\_orders\] Asymptotic orders of centred estimating function derivatives: $$H_{R_a} = l_{R_a} - \mu_{R_a} = O_p(n ^{1/2}) \, ,$$ where $\mu_{R_a} = \expect_G(l_{R_a})$, $l_{R_a} = l_{R_a}(\bar{{}\theta})$ and $\mu_{R_a} = \mu_{R_a}(\bar{{}\theta})$ exist for $a \in \{1, \ldots, 5\}$.

4.  \[moment\_orders\] Asymptotic orders of joint central moments of estimating functions and their derivatives: $$\nu_{R_{a_1}, S_{a_2}, \ldots, T_{a_b}} = \left\{
          \begin{array}{ll} O(n^{(b - 1)/2})\,, & \text{if $b$ is odd} \\
    O(n^{b/2})\,, & \text{if $b$ is even}
          \end{array} \right. \, ,$$ where $\nu_{R_{a_1}, S_{a_2}, \ldots, T_{a_b}} = \expect_G
    (H_{R_{a_1}} H_{S_{a_2}} \cdots H_{T_{a_b}})$ are joint central moments of estimating functions and their derivatives, with $R_{a_1}, S_{a_2}, \ldots, T_{a_b}$ being subsets of $a_1, a_2,
    \ldots, a_b > 0$ integers, respectively.

5.  \[invertibility\] The matrix with elements $\mu_{rs}$ $(r,s = 1,\ldots,p)$ is invertible.

Assumption \[consistency\] is a working assumption that we make about the unbiasedness of the estimating functions and the consistency of the $M$-estimators. Consistency can sometimes be shown to hold under weak assumptions about $G$ and the asymptotic unbiasedness of the estimating functions; see, for example, @vaart:1998 [Section 5.2] and @huber+ronchetti:2009 [Section 6.2], for theorems on the consistency of $M$-estimators. We assume that there is an index $n$, which is typically, but not necessarily, the number of observations, that measures the rate the information about the parameter ${{}\theta}$ accumulates, and that the difference $\hat{{}\theta} - \bar{{}\theta}$ is of order $n^{-1/2}$ in probability.

Assumption \[smoothness\] allows taking derivatives of the estimating functions across the parameter space when constructing the stochastic Taylor expansions required for the derivation of the empirical bias-reducing adjustments to the estimating functions in Section \[sec:theory\]. Such an assumption covers many well-used estimating functions, like the ones arising in quasi-likelihood estimation, estimation using generalized estimating equations, and maximum likelihood and maximum composite likelihood estimation for a wide range of models. The smoothness assumption does not cover, though, settings where the estimating function or one of its first few derivatives are non-differentiable at particular points in the parameter space. Examples of this kind are the estimating functions for quantile regression and robust regression with Huber loss; see @koenker:2005 and @huber+ronchetti:2009 for textbook-length expositions of topics in quantile and robust regression, respectively.

Assumptions \[derivative\_orders\] and \[moment\_orders\] ensure the existence of the expectations, under the underlying process $G$, of products of estimating functions and their derivatives, and that $\sqrt{n}$-asymptotic arguments are valid. Assumption \[invertibility\] is a technical assumption to ensure that the expectation of the jacobian of the estimating function is invertible, when inverting the stochastic Taylor expansions in Section \[sec:theory\], and is typically assumed for estimation using maximum likelihood and estimating equations [see, for example, @boos+stefanski:2013 Section 7.7].

Adjusted estimating equations for bias reduction {#sec:theory}
================================================

Asymptotic bias {#asymptotic_bias}
---------------

Under the assumptions of Section \[assumptions\] it can be shown that the bias of $\hat{{}\theta}$ is $$\label{bias}
  \expect_G(\hat{{}\theta} - \bar{{}\theta}) = O(n^{-1}) \, .$$ This provides some reassurance that, as the information about the parameter ${{}\theta}$ grows, the bias from $\bar{{}\theta}$ will converge to ${{}0}_p$. Nevertheless, the finite sample bias of $\hat{{}\theta}$ is typically not zero. Under the same assumptions, it is also possible to write down (\[bias\]) in the more explicit form $\expect_G(\hat{{}\theta} - \bar{{}\theta}) = {{}b}(\bar{{}\theta}) + O(n^{-3/2})$, where ${{}b}(\bar{{}\theta})$ depends on joint moments of estimating functions and their derivatives under $G$. It is tempting, then, to replace $\hat{{}\theta}$ with a new estimator $\hat{{}\theta} - {{}b}(\hat{{}\theta})$, on the basis that, under \[consistency\], ${{}b}(\hat{{}\theta}) \stackrel{p}{\longrightarrow} {{}b}(\bar{{}\theta})$, and hoping that the new estimator will have better bias properties.

The estimator $\hat{{}\theta} - {{}b}(\hat{{}\theta})$ has been shown to, indeed, have better bias properties when estimation is by maximum likelihood and the model is correctly specified [see, Section 10 @efron:1975 for a proof]. By the model being correctly specified we mean that the unknown joint distribution function $G(z_1, \ldots, z_m)$ is assumed to be a particular member of the family of distributions specified fully by $F_i({{}y}_i | {{}x}_i, {{}\theta})$ when forming (\[objective\_function\]). In that particular case, it is also possible to evaluate ${{}b}({{}\theta})$ and, hence, compute $\hat{{}\theta} - {{}b}(\hat{{}\theta})$ in light of data, because the expectations involved in the joint null moments are with respect to the modelling assumption. This is also the basis of more refined bias reduction methods, like the adjusted score function approach that has been derived in @firth:1993 and explored further in @kosmidis+firth:2009.

In the more general setting of Section \[modelling\_setting\], where the model can be only partially specified, naive evaluation of ${{}b}({{}\theta})$ using $F_i({{}y}_i | {{}x}_i, {{}\theta})$ not only does not lead to reduction of bias, in general, but it can also inflate the bias; see, for example, @lunardon+scharfstein:2017 [Section 2.1] who illustrate the impact of using an incorrect bias function in the estimation of log-odds in longitudinal settings. Furthermore, even in the case that the researcher is comfortable to assume that $F_i({{}y}_i | {{}x}_i, {{}\theta})$ is correctly specified, the applicability of standard bias-reduction methods is hampered whenever $F_i({{}y}_i | {{}x}_i, {{}\theta})$ is impossible or impractical to compute in closed form. In such cases, estimation and inference is based on pseudo likelihoods or estimating functions; an example of this kind is within the framework of max-stable processes for which the evaluation of the joint density becomes quickly infeasible when the number of site locations increases [@davison+gholamrezaee:2011].

In the remainder of this section, we show that the gap between improved estimation methods and applied modelling using general $M$-estimators can be bridged either through simple empirical bias-reducing adjustments to the estimating functions (\[estimating\_eqautions\]) that do not involve expectations, or empirical bias-reducing penalties to the objective function (\[objective\_function\]).

Family of bias-reducing adjustments to estimating functions
-----------------------------------------------------------

Suppose that the assumptions in Section \[assumptions\] hold and consider the estimator $\tilde{{}\theta}$ that results from the solution of the adjusted estimating equations $$\label{adjusted_estimating_equations}
  \sum_{i=1}^k {{}\psi}^i({{}\theta}) + {{}A}({{}\theta}) = {{}0}_p \, ,$$ where both ${{}A}({{}\theta}) = {{}A}({{}\theta}, {{}Y}, X)$ and its derivatives with respect to ${{}\theta}$ are $O_p(1)$ as $n$ grows.

Using assumptions \[consistency\]-\[derivative\_orders\] and index notation, with the indices taking values in the set $\{1,\ldots,p\}$, a calculation similar to that in @mccullagh:2018 [Section 7.3] can be used to show that the expansion of ${{}0}_p = \sum_{i=1}^k {{}\psi}^i(\tilde{{}\theta}) + {{}A}(\tilde{{}\theta})$ about $\bar{{}\theta}$ results in a stochastic Taylor expansion for $\tilde{{}\theta} - \bar{{}\theta}$ of the form $$\begin{aligned}
  \label{taylor}
  \tilde{\theta}^{r} - \bar{\theta}^{r} & =  H^{r} + H^{a} H_{a}^{r} + \frac{1}{2} H^{a} H^{b} \mu_{ab}^{r} + A^{r} +\\ \notag
                                   &+ H^{a} H_{a}^{b} H_{b}^{r} + \frac{1}{2} H^{a} H^{b} H_{c}^{r} \mu_{ab}^{c} + \frac{1}{2} H^{a} H^{b} H_{b}^{c} \mu_{ac}^{r}+ \\ \notag
                                  & + \frac{1}{2} H^{a} H^{b} H_{a}^{c} \mu_{cb}^{r} + \frac{1}{4} H^{a} H^{b} H^{c} \mu_{bc}^{d} \mu_{ad}^{r} + \frac{1}{4} H^{a} H^{b} H^{c} \mu_{ab}^{d} \mu_{dc}^{r} + \\ \notag
                                  & + \frac{1}{2} H^{a} H^{b} H_{ab}^{r} + \frac{1}{6} H^{a} H^{b} H^{c} \mu_{abc}^{r} +\\ \notag
                                  & + A^{a} H_{a}^{r} + \frac{1}{2} A^{a} H^{b} \mu_{ba}^{r} + \frac{1}{2} A^{a} H^{b} \mu_{ab}^{r} + A_{a}^{r} H^{a} + O_p(n^{-2}) \, ,\end{aligned}$$ where ${H^{r}_{r_1\ldots r_a}} = - {\mu^{rs}_{}} {H^{}_{sr_1\cdots r_a}}$, ${\mu^{r}_{r_1\cdots r_a}} = - {\mu^{rs}_{}} {\mu^{}_{sr_1\cdots r_a}}$, and ${A^{r}_{r_1\cdots r_a}} = - {\mu^{rs}_{}} {A^{}_{sr_1\cdots r_a}}$, with $\mu^{rs}$ denoting the matrix inverse of $\mu_{rs}$ (assumption \[invertibility\]) and $A_{r_1\cdots r_a} = \partial^{a-1} A_{r_1}({{}\theta}) / \partial
\theta^{r_2} \cdots \partial \theta^{r_a}$.

Taking expectations with respect to the underlying distribution $G$ on both sides of (\[taylor\]), assumption \[moment\_orders\] gives that the bias of $\tilde{{{}\theta}}$ is $$\expect_G{(\tilde{\theta}^{r} - \bar{\theta}^{r})} = -{\mu^{ra}_{}}
  \expect_G (A_a) + \frac{1}{2} {\mu^{ra}_{}} {\mu^{bc}_{}}
  \left({\nu^{}_{ab,c}} - {\mu^{de}_{}} {\nu^{}_{c,e}}
    {\mu^{}_{abd}} \right) + O(n^{-3/2})\, ,$$ where all terms in the right-hand side are understood as being evaluated at $\bar{{}\theta}$.

The above expansion for the bias implies that use of any adjustment ${{}A}$ with $$\label{adjustment}
  \expect_G(A_r) = \frac{1}{2} {\mu^{ab}_{}} \left( 2 {\nu^{}_{ra,b}} - {\mu^{cd}_{}} {\nu^{}_{b,d}} {\mu^{}_{rac}} \right) + O(n^{-1/2}) \, ,$$ in (\[adjusted\_estimating\_equations\]), will result in estimators with bias $\expect_G{(\tilde{\theta}^{r} - \bar{\theta}^{r})} = O(n^{-3/2})$, which is smaller, asymptotically, than the bias of $\hat{{}\theta}$ in (\[bias\]). Hence, expression (\[adjustment\]) defines a family of bias-reducing adjustments to estimating functions.

A clear candidate for ${{}A}$ has as $r$th component the first term in the right hand side of (\[adjustment\]). If we assume that the model $F_i({{}y}_i | {{}x}_i, {{}\theta})$ fully specifies $G$ and the estimation method is maximum likelihood, then the Bartlett relations $\mu_{cd} + \nu_{c, d} = 0$ hold [see, for example, @pace+salvan:1997 Section 9.2], and the adjustment (\[adjustment\]) becomes ${\mu^{ab}_{}} \left( 2 {\nu^{}_{ra,b}} + {\mu^{}_{rab}}
\right)/2$. The latter expression coincides with the bias-reducing score adjustment derived in @firth:1993.

For the more general setting of Section \[modelling\_setting\], however, the underlying distribution $G$ is at most only partially specified through $F_i({{}y}_i | {{}x}_i, {{}\theta})$, and the expectations involved in the right-hand side of (\[adjustment\]) cannot be computed. It turns out that, there is a family of empirical adjustments that can always be implemented and delivers bias reduction of $M$-estimators in that more general framework.

Empirical bias-reducing adjustments {#empirical_adjustments}
===================================

Derivation
----------

Consider data-dependent quantities ${l^{}_{r|s}}({{}\theta})$ and ${l^{}_{rs|t}}({{}\theta})$ such that $\expect_G({l^{}_{r|s}}) = {\nu^{}_{r, s}}$ and $\expect_G({l^{}_{rs|t}}) = {\nu^{}_{rs, t}}$, with ${H^{}_{r | s}} = {l^{}_{r | s}} - {\nu^{}_{r, s}} = O_p(n^{1/2})$ and ${H^{}_{rs | t}} = {l^{}_{rs | t}} - {\nu^{}_{rs, t}} =
O_p(n^{1/2})$. Then, a simple calculation under the assumptions of Section \[assumptions\] gives that equation (\[adjustment\]) is satisfied by an empirical adjustment of the form $$\label{general_empirical_adjustment}
  A_r({{}\theta}) = \frac{1}{2} {l^{ab}_{}}({{}\theta}) \left\{ 2 {l^{}_{ra|b}}({{}\theta}) - {l^{cd}_{}}({{}\theta}) {l^{}_{b|d}}({{}\theta}) {l^{}_{rac}}({{}\theta})\right\} \, ,$$ where ${l^{st}_{}}({{}\theta})$ is the matrix inverse of ${l^{}_{ts}}({{}\theta})$. The matrix form of expression (\[general\_empirical\_adjustment\]) sets the $r$th element of the vector of empirical bias-reducing adjustments to $$\label{general_empirical_adjustment_matrix}
  A_r({{}\theta}) = - \trace\left\{{{}j}({{}\theta})^{-1} {{}d}_r({{}\theta}) \right\} -
  \frac{1}{2}\trace\left[ {{}j}({{}\theta})^{-1} {{}e}({{}\theta}) \left\{{{}j}({{}\theta})^{-1}\right\}^\top {{}u}_r({{}\theta}) \right] \, ,$$ where ${{}u}_r({{}\theta})=\sum_{i=1}^k \nabla \nabla^\top \psi_r^i({{}\theta})$, and ${{}j}({{}\theta})$ is the matrix with $s$th row $-\sum_{i=1}^k \nabla \psi_s^i({{}\theta})$, assumed to be invertible but not necessarily symmetric. The matrix ${{}j}({{}\theta})$ coincides with the negative hessian matrix $\nabla \nabla^\top l({{}\theta})$ when estimation is through the maximisation of an objective function like (\[objective\_function\]), and is the observed information matrix when estimation is through maximum likelihood. The $p \times p$ matrices ${{}e}({{}\theta})$ and ${{}d}_r({{}\theta})$ correspond to the quantities ${l^{}_{r|s}}({{}\theta})$ and ${l^{}_{rs|t}}({{}\theta})$, respectively, and need to be exactly defined in order to use adjustment (\[general\_empirical\_adjustment\]) in estimation problems.

Specification of p and p under independence {#independence}
-------------------------------------------

If the assumption that ${{}Y}_1, \ldots, {{}Y}_k$ are independent under $G$ is plausible, then $$\nu_{r, s} = \sum_{i = 1}^k E_G( \psi_r^i \psi_s^i) + \sum_{i \ne j} E_G( \psi_r^i) E_G(\psi_s^j) \, .$$ Assumption \[consistency\] on the unbiasedness of the estimating function implies that the second term in the right-hand side of the above expression is zero. Hence, $l_{r|s} = \sum_{i = 1}^k \psi_r^i \psi_s^i$ is such that $\expect_G(l_{r|s}) = \nu_{r, s}$ and, under assumption \[derivative\_orders\], ${H^{}_{r | s}} = {l^{}_{r | s}} - {\nu^{}_{r, s}} = O_p(n^{1/2})$.

In matrix notation, ${{}e}({{}\theta})$ can be taken to be the $p \times p$ matrix ${{}\Psi}({{}\theta})^\top {{}\Psi}({{}\theta})$, where ${{}\Psi}({{}\theta})$ is the $n \times p$ matrix with $(i, s)$th element $\psi_s^i({{}\theta})$. In other words, the $(s, t)$th element of ${{}e}({{}\theta})$ is $$\label{e_mat}
  [{{}e}({{}\theta})]_{st} = \sum_{i = 1}^k \psi_s^i({{}\theta}) \psi_t^i({{}\theta}) \,.$$ A similar argument gives that ${{}d}_r({{}\theta})$ can be taken to be the $p \times p$ matrix $\tilde{{{}\Psi}}_r({{}\theta})^\top {{}\Psi}({{}\theta})$ where $\tilde{{{}\Psi}}_r({{}\theta})$ is the $n \times p$ matrix with $(i, s)$th element $\partial \psi_r^i / \partial \theta^s$. In other words, the $(s, t)$th element of ${{}d}_r({{}\theta})$ is $$\label{d_mat}
  [{{}d}_r({{}\theta})]_{st} = \sum_{i = 1}^k \left\{ \frac{\partial}{\partial \theta^s} \psi_r^i({{}\theta}) \right\} \psi_t^i({{}\theta}) \, .$$

The above expressions can directly be used in settings that involve observations of $k$ independent random vectors with dependent components. Examples of such settings are the generalized estimating equations in @liang+zeger:1986 for estimating marginal regression parameters for correlated responses, and the composite likelihood approach in @varin+czado:2010 for estimating multivariate probit models with random effects to account for serial dependence in longitudinal settings with binary and ordinal outcomes. Section \[autoregressive\_probit\] considers a case study using the latter models.

Reduced-bias $M$-estimation for partially-specified models
----------------------------------------------------------

Under the assumptions of Section \[assumptions\] and the assumptions about $G$ required for the specification of ${{}e}({{}\theta})$ and ${{}d}_r({{}\theta})$ in Section \[independence\], the solution of the adjusted estimating equations $\sum_{i=1}^k{{}\psi}^i({{}\theta}) + {{}A}({{}\theta}) =
0_p$, with ${{}A}({{}\theta})$ as in (\[general\_empirical\_adjustment\]), will result in RB$M$-estimators that have $O(n^{-3/2})$ bias. This is a major advantage over past bias-reduction methods [see, e.g. @cordeiro+mccullagh:1991; @firth:1993] whose applicability is limited to cases where $\sum_{i=1}^k{{}\psi}^i({{}\theta})$ is the gradient of the log-likelihood function of a correctly-specified model.

The example below demonstrates the benefits of using the empirical bias-reducing adjustment in the estimation of a ratio of means from realizations of independent pairs of random variables without any further assumptions on the joint distributions of the pairs.

\[ratio\_estimation\]

[**Ratio of two means**]{} Consider a setting where independent pairs of random variables $(X_1, Y_1), \ldots, (X_n, Y_n)$ are observed, and suppose that interest is in the ratio of the mean of $Y_i$ to the mean of $X_i$, that is $\theta = \mu_Y / \mu_X$, with $\mu_X = \expect_G(X_i)$ and $\mu_Y = \expect_G(Y_i) \ne 0$ $(i = 1, \ldots, n)$.

Assuming that sampling is from an infinite population, one way of estimating $\theta$ without any further assumptions about the joint distribution of $(X_i, Y_i)$ is to set up an unbiased estimating equation of the form (\[estimating\_eqautions\]), with $\psi^i(\theta) = Y_i - \theta X_i$. Then, the $M$-estimator is $$\label{M_ratio}
    \hat\theta = \arg \mathop{\rm solve}_{\theta \in \Re} \left\{\sum_{i = 1}^n \psi^i(\theta) = 0 \right\} = \frac{s_Y}{s_X} \, ,$$ where $s_X = \sum_{i = 1}^n X_i$ and $s_Y = \sum_{i = 1}^n Y_i$. The estimator $\hat\theta$ is generally biased, as can be shown, for example, by an application of the Jensen inequality assuming that $X_i$ is independent of $Y_i$, and efforts have been made in reducing its bias; see, for example, @durbin:1959, for an early work in that direction using the jackknife.

For the estimating function for the ratio of means, $j(\theta) = s_X$, $u(\theta) = 0$, and, using (\[e\_mat\]) and (\[d\_mat\]), $e(\theta) = s_{YY} + \theta^2 s_{XX} - 2\theta s_{XY}$ and $d(\theta) = -s_{XY} + \theta s_{XX}$, where $s_{XX} = \sum_{i = 1}^n X_i^2$, $s_{YY} = \sum_{i = 1}^n Y_i^2$ and $s_{XY} = \sum_{i = 1}^n X_i Y_i$. So, the empirical bias-reducing adjustment in (\[general\_empirical\_adjustment\]) is $s_{XY}/s_{X} - \theta s_{XX}/s_{X}$, and the solution of the adjusted estimating equations results in the RB$M$-estimator $$\label{BR_ratio}
    \tilde\theta = \arg \mathop{\rm solve}_{\theta \in \Re} \left\{\sum_{i = 1}^n \psi^i(\theta) + A(\theta) = 0 \right\} = \frac{s_Y + \frac{s_{XY}}{s_{X}}}{s_X + \frac{s_{XX}}{s_{X}}} \, .$$ In this case, the empirical bias-reducing adjustment has the side-effect of producing an estimator that is more robust to small values of $s_X$ than the standard $M$-estimator is. In particular, as $s_X$ becomes smaller in absolute value, $\hat\theta$ diverges, while $\tilde{\theta}$ converges to $s_{XY}/s_{XX}$, which is the slope of the regression line through the origin of $y$ on $x$. As a result, when $\mu_X$ is small in absolute value, $\tilde\theta$ has not only smaller bias, as granted by the developments in the current paper, but also smaller variance than $\hat\theta$, and, hence, smaller mean squared error.

To illustrate the performance of the RB$M$-estimator $\tilde\theta$ of the ration $\theta$ we assume that $(X_1, Y_1), \ldots, (X_n, Y_n)$ are independent random vectors from a bivariate distribution constructed through a gaussian copula with correlation $0.5$, to have an exponential marginal with rate $1/2$ for $X_i$ and a normal marginal for $Y_i$ with mean $10$ and variance $1$. For each $n \in \{10, 20, 40, 80, 160, 320\}$, we simulate $N_n = 250 \times 2^{16} / n$ samples. Calibrating the simulation size in this way guarantees a fixed simulation error for the simulation-based estimate of the bias $\expect_G(\hat\theta - \theta)$ of the $M$-estimator $\hat\theta$. For each sample we estimate $\theta$ using the $M$-estimator $\hat\theta$, the RB$M$-estimator $\tilde\theta$, and the jackknife estimator $$\label{Ja_ratio}
    \theta^* = n \hat\theta - \frac{n - 1}{n} \sum_{i = 1}^n \hat\theta_{(-i)} \, ,$$ where $\hat\theta_{(-i)} = \sum_{j \ne i} Y_i/\sum_{j \ne i}
  X_i$. Note here that use of the bias-reducing adjusted score function in @firth:1993 requires expectations of products of log-likelihood derivatives, which require fully specifying the bivariate distribution for $(X_i, Y_i)$. Despite the simplicity of the current setting, even if one could confidently specify that distribution, the required expectations involve non-trivial analytic calculations and may not even be available in closed-form.

  -- ---------------- ------- ------- ------ ------ ------ ------
                                                           
                                                           
     $\hat\theta$        0.53    0.25   0.12   0.06   0.03   0.01
     $\theta^*$         -0.07   -0.01   0.00   0.00   0.00   0.00
     $\tilde\theta$      0.17    0.05   0.01   0.00   0.00   0.00
     $\hat\theta$        3.76    1.47   0.65   0.31   0.15   0.07
     $\theta^*$          3.04    1.29   0.61   0.30   0.15   0.07
     $\tilde\theta$      3.12    1.31   0.61   0.30   0.15   0.07
     $\hat\theta$        1.37    0.91   0.62   0.44   0.30   0.21
     $\theta^*$          1.31    0.88   0.61   0.43   0.30   0.21
     $\tilde\theta$      1.29    0.88   0.61   0.43   0.30   0.21
     $\hat\theta$        0.46    0.47   0.48   0.48   0.49   0.49
     $\theta^*$          0.60    0.56   0.54   0.53   0.52   0.51
     $\tilde\theta$      0.54    0.54   0.53   0.52   0.52   0.51
  -- ---------------- ------- ------- ------ ------ ------ ------

  : Simulation-based estimates of the bias (B; $\expect_G\left(\theta_n - \theta\right)$), mean squared error (MSE; $\expect_G\left\{(\theta_n - \theta)^2\right\}$), mean absolute error (MAE; $\expect_G\left(\left|\theta_n - \theta\right|\right)$), and probability of underestimation (PU; $P_G(\theta_n < \theta)$) for the estimation of a ratio of means in the setting of Example \[ratio\_estimation\]. The estimator $\theta_n$ is either the $M$-estimator $\hat\theta$ in (\[M\_ratio\]), the Jackknife estimator $\theta^*$ in (\[Ja\_ratio\]) or the RB$M$-estimator $\tilde\theta$ in (\[BR\_ratio\]). The simulation-based estimates are computed from $N_n = 250 \times 2^{16} / n$ simulated samples with $n \in \{10, 20, 40, 80, 160, 320\}$ independent pairs each from a bivariate distribution constructed through a gaussian copula with correlation $0.5$, to have an exponential marginal with rate $1/2$ for $X_i$ and a normal marginal for $Y_i$ with mean $10$ and variance $1$ ($\theta = 5$). Figures are reported in 2 decimal places, and a figure of $0.00$ indicates an estimated bias between $-0.005$ and $0.005$. The simulation error for the estimates of the bias is between $0.001$ and $0.002$.

\[ratio\_simulation\]

Table \[ratio\_simulation\] shows the simulation-based estimates of the bias, mean-squared error, mean absolute deviation, and probability of underestimation for $\hat\theta$, $\tilde\theta$, and $\theta^*$. As with the jackknife, use of the empirical bias-reducing adjustment results in a marked reduction of the bias. The estimated slopes of the regression lines of the logarithm of the absolute value of the estimated biases for $\hat\theta$, $\theta^*$, and $\tilde\theta$ on the values of $n$ are $-1.046$, $-1.445$, and $-1.610$, respectively, which are in close agreement to the theoretical slopes of $-1$, $-3/2$, $-3/2$.

Reduction of the bias in this setting leads also in marked reduction in mean squared error and mean absolute deviation, with $\tilde\theta$ and $\theta^*$ performing similarly, and significantly better than $\hat\theta$. The RB$M$-estimator $\tilde{\theta}$ also has probability of underestimation closer to $0.5$, hence it is closer to being median unbiased than is $\theta^*$.

Asymptotic distribution of reduced-bias $M$-estimators and inference {#sec:inference}
====================================================================

Under the assumptions of Section \[assumptions\], the argument in @stefanski+boos:2002 [Section 2] applied to the first term of the stochastic Taylor expansion in (\[taylor\]) with ${{}A}({{}\theta}) = {{}0}_p$ gives that the $M$-estimator $\hat{{}\theta}$ is such that $${{}Q}(\bar{{}\theta})^{1/2} (\hat{{}\theta} - \bar{{}\theta}) \stackrel{d}{\longrightarrow} N_p({{}0}_p, {{}I}_p)$$ as $n$ increases, where ${{}Q}({{}\theta}) = {{{}V}({{}\theta})}^{-1}$ with ${{}V}({{}\theta}) = {{}B}({{}\theta})^{-1} {{}M}({{}\theta}) \{{{}B}({{}\theta})^{-1}\}^\top$. In the latter expression, ${{}M}({{}\theta}) = \expect_G[ \sum_{i=1}^k \sum_j^k {{}\psi^i}({{}\theta})\{ {{}\psi}^j({{}\theta})\}^\top ]$, and ${{}B}({{}\theta})$ is a $p \times p$ matrix with $r$th row $-\sum_{i=1}^k \expect_G\{\nabla \psi_r^i({{}\theta}) \}$. The notation ${{}D} \stackrel{d}{\longrightarrow} N_p({{}0}_p, {{}I}_p)$ indicates that the random $p$-vector ${{}D}$ converges in distribution to a multivariate Normal with mean ${}0_p$ and variance-covariance matrix the $p \times p$ identity matrix ${{}I}_p$.

From the stochastic Taylor expansion (\[taylor\]), it becomes directly apparent that $\tilde{{{}\theta}} - \bar{{{}\theta}}$ and $\hat{{{}\theta}} - \bar{{{}\theta}}$ have exactly the same $O_p(n^{-1/2})$ term in their expansions because ${{}A}({{}\theta}) = O_p(1)$. Hence, the RB$M$-estimator $\tilde{{}\theta}$ is also such that $$\label{asymptotic_normality}
  {{}Q}(\bar{{}\theta})^{1/2} (\tilde{{}\theta} - \bar{{}\theta}) \stackrel{d}{\longrightarrow} N_p({{}0}_p, {{}I}_p) \, .$$

An implication of (\[asymptotic\_normality\]) is that $\hat{{{}V}}({{}\theta}) = {{}j}({{}\theta})^{-1} {{}e}({{}\theta}) \left\{{{}j}({{}\theta})^{-1}\right\}^\top$ evaluated at ${{}\theta}=\tilde{{}\theta}$ is a consistent estimator of the variance-covariance matrix of $\tilde{{}\theta}$. This is exactly as in the case where $\hat{{{}V}}(\hat{{}\theta})$ is used as an estimator of the variance-covariance matrix of $\hat{{}\theta}$ in the framework of $M$-estimation [see, for example, @stefanski+boos:2002 Section 2]. It is important to note that the expression for $\hat{{{}V}}({{}\theta})$ appears unaltered in the second term of the right-hand-side of expression (\[general\_empirical\_adjustment\_matrix\]) for the empirical bias-reducing adjustment. As a result, the value of ${{}V}(\tilde{{}\theta})$ and, hence, that of estimated standard errors for the parameters, is, in general, readily available at the last step of the iterative process that is used to solve the adjusted estimating equations. Section \[sec:implementation\], provides more details on the implementation of solvers for the bias-reducing adjusted estimating equations.

In addition, if the model is correctly specified, and $l({{}\theta})$ in (\[objective\_function\]) is the log-likelihood, then the second Bartlett identity gives ${{}M}({{}\theta}) = {{}B}({{}\theta})$, which implies that ${{}Q}({{}\theta})$ is the expected information matrix. As a result, the RB$M$-estimator is asymptotically efficient, exactly as the maximum likelihood estimator and the reduced-bias estimator in @firth:1993 are.

  -------------------------------------- ------ ------ ------ ------ ------ ------
                                                                            
                                                                            
  $\var_G(\hat\theta)$                     3.49   1.41   0.64   0.30   0.15   0.07
  $\expect_G\{\hat{V}(\hat\theta)\}$       2.55   1.20   0.59   0.29   0.14   0.07
  $\var_G(\tilde\theta)$                   3.09   1.31   0.61   0.30   0.15   0.07
  $\expect_G\{\hat{V}(\tilde\theta)\}$     2.16   1.10   0.56   0.28   0.14   0.07
  -------------------------------------- ------ ------ ------ ------ ------ ------

  : Simulation-based estimates of the variances of the $M$-estimator $\hat\theta$ and the RB$M$-estimator $\tilde\theta$ ($\var_G(\hat\theta)$ and $\var_G(\tilde\theta)$, respectively) and of the mean of $\hat{V}(\hat\theta)$ and $\hat{V}(\tilde\theta)$ ($\expect_G\{\hat{V}(\hat\theta)\}$ and $\expect_G\{\hat{V}(\tilde\theta)\}$, respectively), from the simulation study of Example \[ratio\_estimation\]. See, also Table \[ratio\_simulation\].

\[ratio\_inference\]

[**Ratio of two means (continued)**]{} Table \[ratio\_inference\] shows the estimates of the actual variances of $\hat\theta$ and $\tilde\theta$ for $n \in \{10, 20, 40, 80, 160, 320\}$ and the estimate of the mean of $\hat{V}(\hat\theta)$ and $\hat{V}(\tilde\theta)$, respectively. As expected by the above discussion, the large sample approximations to the variance of the estimator converge to the actual variances as the sample size increases.

Another implication of (\[asymptotic\_normality\]) is that asymptotically valid inferential procedures, like hypothesis tests and confidence regions for the model parameters, can be constructed based on the Wald-type and generalized score approximate pivots of the form $$\begin{aligned}
  \label{pivots}
  W_{(e)}({{}\theta}) & = (\tilde{{}\theta} - {{}\theta})^\top \left\{\hat{{{}V}}(\tilde{{}\theta}) \right\}^{-1}(\tilde{{}\theta} - {{}\theta}) \, , \\ \notag
  W_{(s)}({{}\theta}) & = \left\{\sum_{i = 1}^k{{}\psi}^i({{}\theta}) + {{}A}({{}\theta})\right\}^\top  \left\{{{}e}(\tilde{{}\theta})\right\}^{-1}  \left\{\sum_{i = 1}^k{{}\psi}^i({{}\theta}) + {{}A}({{}\theta})\right\}\,,\end{aligned}$$ respectively, which, asymptotically, have a $\chi^2_p$ distribution. These pivots are direct extensions of the Wald-type and generalized score pivots, respectively, that are typically used in $M$-estimation [see, @boos:1992 for discussion about generalized score tests in $M$-estimation].

Empirical bias-reducing penalties to objective functions {#sec:penalties}
========================================================

Bias-reducing penalized objectives {#sec:penalty}
----------------------------------

When estimation is through the maximisation of (\[objective\_function\]), ${{}j}({{}\theta})$ is a symmetric matrix. Then, an extra condition that determines how $l_{r|s}$ behaves under differentiation guarantees that the bias-reducing adjustment (\[general\_empirical\_adjustment\]) always corresponds to an additive penalty to the objective function used for estimation. In particular, if $l_{r|s}$ admits a chain rule under differentiation, that is $$\label{diff_condition}
  \frac{\partial}{\partial {\theta^{a}_{}}} {l^{}_{r|s}} = {l^{}_{ar | s}} + {l^{}_{r | as}} \, ,$$ then bias reduction through empirical bias-reducing adjustments is formally equivalent to the maximization of a penalized objective function of the form $$\label{penalized_likelihood}
  l({{}\theta}) - \frac{1}{2}\trace \left\{{{}j}({{}\theta})^{-1} {{}e}({{}\theta}) \right\} \, ,$$ assuming that the maximum exists. This result greatly facilitates implementation of bias reduction for a much wider class of models and estimation methods than other bias-reduction methods via adjusted score functions [see, for example @firth:1993], where focus is on cases where $l({{}\theta})$ is the log-likelihood function of a correctly-specified model, and a bias-reducing penalized log-likelihood does not always exist.

The matrix ${{}e}({{}\theta})$ that has been derived under the assumption of independence in Section \[independence\] satisfies condition (\[diff\_condition\]).

Bias reduction and model selection based on Kullback-Leibler divergence {#sec:model_selection}
-----------------------------------------------------------------------

Suppose that $l({{}\theta})$ is the log-likelihood function based on an assumed parametric model $F$. @takeuchi:1976 showed that $$\label{tic}
  - 2 l(\hat{{}\theta}) + 2 \trace \left\{{{}j}(\hat{{}\theta})^{-1} {{}e}(\hat{{}\theta}) \right\}$$ is an estimator of the expected Kullabck-Leibler divergence of the underlying process $G$ to the assumed model $F$, where $\hat{{}\theta}$ is the maximum likelihood estimator. Expression (\[tic\]) is widely known as the Takeuchi information criterion (TIC), and, in contrast to the Akaike Information Criterion (AIC; @akaike:1974), is robust against deviations from the assumption that the model used is correct. @claeskens:2008 [Section 2.5] provide a more thorough discussion on TIC and its formal relationship to AIC.

Model selection from a set of parametric models then proceeds by computing $\hat{{}\theta}$ for each model and selecting the model with the smallest TIC value (\[tic\]), or equivalently, the model with the largest value for $$\label{tic2}
  l(\hat{{}\theta}) - \trace \left\{{{}j}(\hat{{}\theta})^{-1} {{}e}(\hat{{}\theta}) \right\} \, .$$ A direct comparison of expressions (\[tic2\]) and (\[penalized\_likelihood\]) reveals a previously unnoticed close connection between bias reduction in maximum likelihood estimation and model selection. Specifically, both bias reduction and TIC model selection rely on exactly the same penalty $\trace \left\{{{}j}({{}\theta})^{-1} {{}e}({{}\theta})
\right\}$, but differ in the strength of penalization; bias reduction is achieved by using half that penalty, while valid model selection requires stronger penalization by using one times the penalty.

As discussed in Section \[sec:inference\], the RB$M$-estimator $\tilde{{}\theta}$ has the same asymptotic distribution as $\hat{{}\theta}$. Then, the derivation of TIC [see, for example, @claeskens:2008 Section 2.3] works also with the RB$M$-estimator in place of the maximum likelihood estimator. As a result, TIC and, under more assumptions, AIC at the RB$M$-estimates are asymptotically equivalent to their typical versions at the maximum likelihood estimates. The same holds for reduced-bias estimators of @firth:1993. In other words, TIC model selection can proceed by selecting the model with the largest value of $$\label{tic3}
  l(\tilde{{}\theta}) - \trace \left\{{{}j}(\tilde{{}\theta})^{-1} {{}e}(\tilde{{}\theta}) \right\} \, ,$$ and AIC model selection using the largest value of $l(\tilde{{}\theta}) - p$. The quantity in (\[tic3\]) is readily available once (\[penalized\_likelihood\]) has been maximized to obtain the RB$M$-estimates; the only requirement for model selection is to adjust, from $1/2$ to $1$, the factor of the value of $\trace \left\{{{}j}({{}\theta})^{-1} {{}e}({{}\theta})
\right\}$ after maximization.

@varin+vidoni:2005 developed a model selection procedure when the objective $l({{}\theta})$ is a composite likelihood [see, @varin+reid+firth:2011 for a review of composite likelihood methods]. The composite likelihood information criterion (CLIC) derived in [@varin+vidoni:2005] has the same functional form as TIC in (\[tic\]). So, the link between model selection and bias reduction exists also when $l({{}\theta})$ is the logarithm of a composite likelihood.

Implementation {#sec:implementation}
==============

Apart from special cases, like the estimation of the ratio of two means in Example \[ratio\_estimation\], and as is the case for general estimating functions, the solution of the bias-reducing adjusted estimating equations (\[adjusted\_estimating\_equations\]) is, typically, not available in closed form, and iterative procedures are used to approximate that solution.

A general iterative procedure of this kind results from a modification of the Newton-Raphson iteration that in the $u$th iteration updates the current estimate ${{}\theta}^{(u)}$ to a new value ${{}\theta}^{(u + 1)}$ as $$\label{quasiNR}
  {{}\theta}^{(u + 1)} :=  {{}\theta}^{(u)} + a_u \left\{{{}j}\left({{}\theta}^{(u)}\right)\right\}^{-1} \left\{ \sum_{i=1}^k {{}\psi}^i\left({{}\theta}^{(u)}\right) + {{}A}\left({{}\theta}^{(u)}\right) \right\}\, ,$$ where $a_u$ is a deterministic sequence of positive constants that can be used to implement various schemes to further control the step size, like step-halving. Iteration (\[quasiNR\]) defines a quasi Newton-Raphson procedure with the correct fixed point, rather than a full Newton-Raphson iteration that would have $a_u = 1$ and the matrix of derivatives of $\sum_{i=1}^k {{}\psi}^i({{}\theta}) + {{}A}({{}\theta})$ in the place of ${{}j}({{}\theta})$. The $M$-estimates from the solution of $\sum_{i=1}^k {{}\psi}^i({{}\theta}) = {{}0}_p$ are obvious starting values for the quasi Newton-Raphson procedure. Candidate stopping criteria include $|{{}\theta}^{(u + 1)} - {{}\theta}^{(u)}|/a_u < \epsilon$ and $||\sum_{i=1}^k {{}\psi}^i\left({{}\theta}^{(u)}\right) + {{}A}\left({{}\theta}^{(u)}\right)||_1 < \epsilon$, for some $\epsilon > 0$, where $||\cdot||_1$ is the L1 norm.

Typically, quasi Newton-Raphson will have first-order convergence to the solution of the adjusted estimating equations, compared to the second-order convergence that full Newton-Raphson has. The advantage of using quasi Newton-Raphson instead of full Newton-Raphson is that all quantities required to implement (\[quasiNR\]) are readily available once an implementation of the empirical bias-reducing adjustments is done.

The fact that the empirical bias-reducing adjustment in (\[general\_empirical\_adjustment\]) depends only on derivatives of estimating functions, enables general implementations by deriving the derivatives $\partial \psi_r^i({{}\theta}) / \partial \theta_s$ and $\partial^2 \psi_r^i({{}\theta}) / \partial \theta_s \partial
\theta_t$ $(r, s, t = 1, \ldots, p)$, either analytically or by using automatic differentiation techniques [@griewank+walther:2008]. Those derivatives can be combined together to produce ${{}u}_r({{}\theta})$, ${{}j}({{}\theta})$, ${{}e}({{}\theta})$ and ${{}d}_r({{}\theta})$, and, then, matrix multiplication and a numerical routine for matrix inversion can be used for an easy, general implementation of (\[general\_empirical\_adjustment\_matrix\]).

For implementations using automatic differentiation, in particular, the only required input is an appropriate implementation of the contributions ${{}\psi}^i({{}\theta})$ to the estimating functions. The automatic differentiation routines will, then, produce implementations of the required first and second derivatives of the contributions. The GEEBRA Julia package (<https://github.com/ikosmidis/GEEBRA.jl>) provides a proof-of-concept of such an implementation.

When estimation is through the maximization of an objective function, the RB$M$-estimates can be computed using general numerical optimization procedures for the maximization of the penalized objective function (\[penalized\_likelihood\]), like those provided by the `optim` function in R or the Optim Julia package, that can operate by numerically approximating the gradient of (\[penalized\_likelihood\]). In such cases, bias reduction can be performed using only routines for matrix multiplication and inversion, and the contributions to the estimating function and their first derivatives for the implementation of ${{}j}({{}\theta})$ and ${{}e}({{}\theta})$. Those derivatives can again be obtained using an automatic differentiation library.

The discussion in Section \[sec:penalty\] implies that using the empirical bias-reducing adjustments in maximum likelihood and maximum composite likelihood estimation is readily available for all models for which there are implementations of TIC and CLIC, respectively [see, e.g., @padoan+bevilacqua:2015 for estimation of random fields based on composite likelihoods].

Generalized linear models {#sec:glms}
=========================

Bias-reducing penalized likelihood
----------------------------------

Consider a sequence of $n$ vectors $(y_1, {{}x}_1^\top)^\top, \ldots, (y_n, {{}x}_n^\top)^\top$, where $y_i \in \mathcal{Y} \subset \Re$, and ${{}x}_i \in \mathcal{X} \subset \Re^p$. Suppose that $y_1, \ldots ,y_n$ are realizations of random variables $Y_1, \ldots, Y_n$ which, conditionally on covariate vectors ${{}x}_1, \ldots, {{}x}_n$, are assumed to be independent and distributed according to a generalized linear model. The $i$th log-likelihood contribution is $$\label{loglikelihood_glm}
  \log f_i(y_i | {{}x}_i, {{}\beta}, \phi) = \frac{m_i}{\phi} \left\{ y_i \theta_i - \kappa(\theta_i) - c_1(y_i) \right\} - \frac{1}{2} a\left(-\frac{m_i}{\phi}\right) \, ,$$ for sufficiently smooth functions $\kappa(\cdot)$, $c_1(\cdot)$, and $a(\cdot)$, where the conditional mean and variance of $Y_i$ are associated to ${}\beta$ and ${{}x}_i$ as $\mu_i = d \kappa(\theta_i) / d\theta_i = h(\eta_i)$ and $\var(Y_i | {{}x}_i) = \phi v(\mu_i) / m_i$, respectively, with $\eta_i = {{}x}_i^\top {{}\beta}$, and a sufficiently smooth function $h(\cdot)$ $(i = 1, \ldots, n)$. The parameter $\phi$ is a known or unknown dispersion parameter, $m_1, \ldots, m_n$ are known observation weights, and $v(\mu_i) = d^2 \kappa(\theta_i) / d \theta_i^2$ is the variance function. A few prominent generalized linear models are binomial logistic regression models, Poisson log-linear models, gamma regression models, and normal linear regression models.

Bias reduction of the maximum likelihood estimator $\hat{{}\beta}$ and $\hat\phi$ for fully-specified generalized linear models has been studied extensively. Landmark studies in this direction include: @cordeiro+mccullagh:1991, who derive a closed-form estimator of the bias of $\hat{{}\beta}$ and $\hat\phi$ and subtract that from the estimates; @kosmidis+firth:2009, who show that iterating @cordeiro+mccullagh:1991’s reweighted least squares reduced-bias estimator results in the solution of the bias-reducing adjusted score equations in @firth:1993, and derive a necessary and sufficient condition for the existence of a bias-reducing penalized likelihood for generalized linear models; and @kosmidis+kennepagui+sartori:2019, who derive the connections between mean and median bias reduction for generalized linear models and propose a unifying implementation through a quasi Fisher scoring procedure.

In contrast to the bias-reduction methods proposed for generalized linear models in @kosmidis+firth:2009, the empirical bias-reducing adjustment to the score function always corresponds to a penalty to the log-likelihood function about ${{}\beta}$ and $\phi$. According to Section \[sec:penalty\], if $\phi$ is unknown, the only ingredients required in the penalty are the observed information matrix about ${{}\beta}$ and $\phi$, ${{}j}^{(G)}({{}\beta}, \phi)$, and the sum of the outer products of the gradient of (\[loglikelihood\_glm\]) across observations, ${{}e}^{(G)}({{}\beta}, \phi)$. The bias-reducing penalized log-likelihood is, then $$\sum_{i = 1}^n \log f_i(y_i | {{}x}_i, {{}\beta}, \phi) -
  \frac{1}{2} \trace \left[ \left\{ {{}j}^{(G)}({{}\beta},
      \phi)\right\}^{-1} {{}e}^{(G)}({{}\beta}, \phi) \right] .$$ The closed-form expressions for ${{}j}^{(G)}({{}\beta}, \phi)$ and ${{}e}^{(G)}({{}\beta}, \phi)$ are given in Appendix \[glm\_appendix\]. If $\phi$ is fixed, as is, for example, for the binomial and Poisson distributions, then the penalty involves only the $({{}\beta}, {{}\beta})$ blocks of ${{}j}^{(G)}({{}\beta}, \phi)$ and ${{}e}^{(G)}({{}\beta}, \phi)$. Some algebra, after plugging the expressions in Appendix \[glm\_appendix\] in (\[penalized\_likelihood\]), shows that the RB$M$-estimator of ${{}\beta}$ results as the maximizer of the bias-reducing penalized log-likelihood $$\label{pen_loglik_glm}
  \sum_{i = 1}^n m_i \left\{ y_i \theta_i - \kappa(\theta_i) - \frac{1}{2} s_i \frac{d_i}{v_i}(y_i - \mu_i) \right\}  \, ,$$ where $d_i = d h(\eta_i) / d\eta_i$ and $s_i$ is the $i$th diagonal element of the matrix ${{}X}({{}X}^\top {{}Q} {{}X})^{-1} {{}X}^\top \tilde{{{}W}}$ with ${{}X}$ the $n \times p$ matrix with rows ${{}x}_1, \ldots, {{}x}_n$, and with the diagonal matrices $\tilde{{{}W}}$ and ${{}Q}$ as defined in Appendix \[glm\_appendix\].

  -- -- ----------- -- -------- -------- -------- --------
                                                  
         Parameter                                
         $\beta_1$       -0.047   -0.018   -0.011   -0.008
         $\beta_4$        0.019    0.007    0.004    0.005
         $\beta_5$        0.058    0.027    0.016    0.007
         $\beta_1$       -0.001    0.003    0.002   -0.002
         $\beta_4$       -0.001   -0.001   -0.001    0.002
         $\beta_5$        0.004   -0.001    0.001    0.000
         $\beta_1$       -0.013   -0.000    0.000   -0.002
         $\beta_4$        0.002   -0.001   -0.001    0.002
         $\beta_5$        0.020    0.005    0.002    0.001
         $\beta_1$        0.140    0.064    0.033    0.017
         $\beta_4$        0.133    0.062    0.034    0.016
         $\beta_5$        0.059    0.024    0.011    0.006
         $\beta_1$        0.121    0.060    0.032    0.016
         $\beta_4$        0.118    0.059    0.033    0.016
         $\beta_5$        0.046    0.021    0.011    0.005
         $\beta_1$        0.124    0.060    0.032    0.016
         $\beta_4$        0.120    0.059    0.033    0.016
         $\beta_5$        0.050    0.022    0.012    0.005
  -- -- ----------- -- -------- -------- -------- --------

  : Bias and mean squared error (MSE) of the maximum likelihood (ML) estimator, RB$M$-estimator, and the adjusted scores (AS) estimator in @firth:1993 of $\beta_1$, $\beta_4$ and $\beta_5$ in the probit regression model with $\mu_i = \Phi(\beta_1 + \beta_4 x_{i4} + \beta_4 x_{i5})$ $(i = 1, \ldots, n)$, $n \in \{75, 150, 300, 600\}$. The figures are based on $10\, 000$ samples simulated from the model with $\mu_i = \Phi(-0.5 + 0.5 x_{i4} + 0.5 x_{i5})$. For each value of $n$, the covariate values $x_{i4}$ and $x_{i5}$ $(i = 1, \ldots, n)$ are generated independently and independent of each other from a Bernoulli distribution with probability $3/4$, and an exponential distribution with rate $1$, respectively, and are held fixed across samples. No boundary estimates were encountered, though the figures should be interpreted as simulation-based estimates of the bias and mean-squared error conditional on all estimates being finite.

\[glm\_bias\_mse\]

Probit regression
-----------------

The performance of the RB$M$-estimator is compared here to the performance of the maximum likelihood estimator and the reduced-bias estimator of @firth:1993 in a probit regression model with $\mu_i = \Phi(\beta_1 + \sum_{t = 2}^5 \beta_t x_{it})$ $(i = 1, \ldots, n)$, where $\Phi(\cdot)$ is the cumulative distribution function of a standard Normal distribution. The covariate values $x_{i2}, x_{i3}, x_{i4}, x_{i5}$ $(i = 1, \ldots, n)$ are generated independently and independent of each other from a standard normal distribution, Bernoulli distributions with probabilities $1/4$ and $3/4$, and an exponential distribution with rate $1$, respectively. Note that the necessary and sufficient condition of @kosmidis+firth:2009 [Theorem 1] is not satisfied for this model, and hence, there is no bias-reducing penalized log-likelihood that corresponds to the adjusted score functions of @firth:1993. Nevertheless, the bias-reducing penalized log-likelihood (\[pen\_loglik\_glm\]) is well-defined.

There are $16$ possible nested probit regression models with an intercept $\beta_1$, depending on which of $\beta_2, \ldots, \beta_5$ are zero or non-zero. For $n \in \{75, 150, 300, 600\}$, we simulate $n$ covariate values, as detailed in the previous paragraph, and conditional on those we simulate $10\, 000$ samples of $n$ response values from a probit regression model with ${{}\beta} = (-0.5, 0, 0, 0.5, 0.5)^\top$. For each sample, we estimate all $16$ possible models using maximum likelihood, the adjusted score functions approach in @firth:1993, and maximum bias-reducing penalized likelihood. Maximum likelihood estimates are computed using the `glm` function in R [@rproject], and the @firth:1993 adjusted scores estimates are computed using the `brglm_fit` method from the `brglm2` R package [@brglm2]. The RB$M$-estimates resulting from the numerical maximization of are computed using the `nlm` R function (see the scripts for probit regression supplied in the Supporting Materials).

There were $7$ samples for which separation occurred for at least one of the 16 models, and no separated samples were observed under the data generating model that involves only the intercept, $x_{i4}$ and $x_{i5}$. The detection of separated data sets was done prior to fitting each model using the linear programming algorithms in @konis:2007, as implemented in the `detect_separation` method of the `brglm2` R package [@brglm2]. In those cases, maximum likelihood and maximum penalized likelihood result in estimates on the boundary of the parameter space. In contrast, as has also been observed in the 2007 PhD thesis by Ioannis Kosmidis [@kosmidis:2007], the adjusted score approach of @firth:1993 always resulted in finite estimates.

Table \[glm\_bias\_mse\] shows the simulation-based estimates of the bias and mean squared error of the three estimators for the true model that has $\beta_2 = \beta_3 = 0$, conditionally on the maximum likelihood estimate not being on the boundary of the parameter space. The conditioning is carried out by ignoring the separated samples when taking averages. As is immediately apparent, both the adjusted score and maximum penalized likelihood approaches to bias reduction result in estimators with smaller conditional bias and mean squared errors than the maximum likelihood estimator. The mean squared error of the RB$M$-estimator tends to be slightly larger than that of the adjusted score estimator, but the differences diminish fast as the sample size increases.

Figure \[model\_selection\] shows the selection proportion among the 16 models based on AIC and TIC for each of the three estimators. As expected by the discussion in Section \[sec:model\_selection\], the probability of selecting the model with $\beta_2 = \beta_3 = 0$ increases with the sample size for both information criteria and for all estimation methods. There are only small discrepancies on the selection proportions between estimation methods, which tend to disappear as the sample size increases. Finally, it is worth noting that AIC model selection tends to be more confident on what the true model is than TIC, illustrating less variability in the selected proportions. In the current study, for the larger samples sizes this results in selecting the correct model more often than TIC does. However, in smaller samples sizes, AIC selects the wrong model with $\beta_2 = \beta_3 = \beta_4 = 0$ notably more times than the correct model.

Quasi likelihoods {#sec:quasi}
=================

Preamble
--------

Compared to generalized linear models, much less work has been carried out for reducing the bias of quasi-likelihood estimators. At the time of writing this paper, apart from bias reduction based on resampling schemes, like the bootstrap and the jackknife [see, for example @wu:1986 Section 9], we were unable to find an analytical approach for the reduction of the bias for quasi-likelihood estimators. @heyde:1997 [Section 4.4] notes that an adjusted score approach to bias reduction is possible for quasi likelihoods, but provides no further development in that direction. Also, @paul+zhang:2014 develop a bias-reduction method for generalized estimating equations though their development seems to focus on cases where the distribution of the longitudinal outcomes is correctly specified [see, also @lunardon+scharfstein:2017 who show that @paul+zhang:2014 development is not delivering bias reduction under mispecification of the working intra-subject covariance matrix].

Consider again a sequence of $n$ vectors $(y_1, {{}x}_1^\top)^\top, \ldots, (y_n, {{}x}_n^\top)^\top$, where $y_i \in \mathcal{Y} \subset \Re$, and ${{}x}_i \in \mathcal{X} \subset \Re^p$. Suppose that $y_1, \ldots ,y_n$ are realizations of random variables $Y_1, \ldots, Y_n$ which are assumed independent conditionally on covariate vectors ${{}x}_1, \ldots, {{}x}_n$, and that interest is in the estimation of a parameter ${{}\beta}$ in the assumed mean relation $\mu_i = \expect(Y_i | {{}x}_i) = h(\eta_i)$ with $\eta_i = {{}x}_i^\top {{}\beta}$, for some known function $h(\cdot)$ $(i = 1, \ldots, n)$. Interest is also in allowing for the variance to vary with the mean as $\var(Y_i | {{}x}_i) = \phi v(\mu_i) / m_i$, where $\phi$ is a dispersion parameter, $m_1, \ldots, m_n$ are known observation weights, and $v(\mu_i)$ is a variance function.

Under the assumption that the assumed relationship between the mean and the regression parameters ${{}\beta}$ is correctly specified, a consistent estimator $\hat{{}\beta}$ for ${{}\beta}$ can be obtained by solving the quasi-likelihood equations $\sum_{i= 1}^n \psi^i_s({{}\beta}, \phi) = 0$ $(s = 1, \ldots, p)$ with $$\label{quasilikelihood}
  \psi^i_s({{}\beta}, \phi) = \frac{m_i d_i}{\phi v_i} (Y_i - \mu_i) x_{is} \,,$$ where $d_i = \partial \mu_i / \partial \eta_i$ and $v_i = v(\mu_i)$ [see, @wedderburn:1974; @mccullagh:1983 for introduction and study of quasi-likelihood methods]. @mccullagh:1983 [Section 4] has shown that the quasi-likelihood estimator generalizes the Gauss-Markov optimality of least squares estimators by having the minimum asymptotic variance amongst all estimators resulting from unbiased estimating equations that are linear in $y_i$.

By the discussion in Section \[sec:inference\], if the expression for $\var(Y_i | {{}x}_i)$ is also correctly specified, then the variance-covariance matrix for $\hat{{}\beta}$ can be estimated as $\hat\phi \left\{{{}X}^\top {{}W}(\hat{{{}\beta}})^{-1} {{}X}\right\}^{-1}$, for some estimator $\hat\phi$ of the dispersion parameter, where ${{}W}({{}\beta})$ is a diagonal matrix with $i$th diagonal element $m_i d_i^2/v_i$ . Otherwise, the variance-covariance matrix for $\hat{{}\beta}$ can be estimated by $\hat{{{}V}}(\hat{{}\beta}, \hat\phi)$ as defined in Section \[sec:inference\] or variants of it; see, for example, @firth:1993b [Section 2] for a concise discussion on the estimation of the variance-covariance matrix of $\hat{{}\beta}$.

@mccullagh+nelder:1989 recommend to estimate $\phi = m_i \var(Y_i | {{}x}_i) / v(\mu_i)$ using the moment estimator $$\label{dispersion}
  \hat\phi_R = \frac{1}{R} \sum_{i = 1}^n \frac{m_i}{v(\hat\mu_i)} (Y_i - \hat\mu_i)^2 \, ,$$ with $\hat\mu_i = h(\hat{{}\beta}^\top {{}x}_i)$ and $R = n - p$ to account for the estimation of the $p$-vector of regression parameters ${{}\beta}$. The estimator $\hat\phi_n$ is consistent under the same assumptions that $\hat\phi_{n - p}$ is consistent, but $\hat\phi_{n - p}$ is considered to be less biased than $\hat\phi_n$.

It is well-known that quasi-likelihood methods about ${{}\beta}$ are closely related to maximum likelihood methods for generalized linear models. If the distribution of $Y_i | {{}x}_i$ is fully-specified to be an exponential family with mean $\mu_i$ and variance $\phi v_i /m_i$, then it is straightforward to see that (\[quasilikelihood\]) is the partial derivative of the $i$th log-likelihood contribution (\[loglikelihood\_glm\]) with respect to $\beta_s$. So, the discussion and results below also cover estimation of the parameters of generalized linear models with unknown dispersion estimated by $\hat\phi_R$.

Bias reduction in the estimation of p and $\phi$
------------------------------------------------

There exist two immediate possibilities for the use of the empirical bias-reducing adjustments (\[general\_empirical\_adjustment\_matrix\]) in quasi-likelihood estimation. The first is to adjust the quasi-likelihood equations for ${{}\beta}$, and after solving those, to use $\hat\phi_{n-p}$ in (\[dispersion\]) for the estimation of $\phi$ with $\hat\mu_i$ replaced by $\tilde\mu_i = h({{}x}_i^\top\tilde{{}\beta})$, where $\tilde{{}\beta}$ is the RB$M$-estimator.

The second possibility is to adjust the quasi-likelihood estimating functions for ${{}\beta}$ and the estimating function $\sum_{i = 1}^n\psi_{p + 1}^i({{}\beta}, \phi)$ with $$\label{phi_ee}
  \psi_{p + 1}^i({{}\beta}, \phi) = \frac{m_i}{v(\mu_i)} (Y_i - \mu_i)^2 - \frac{R}{n} \phi \, .$$ Under the assumption that $\expect(Y_i | {{}x}_i)$ and $\var(Y_i | {{}x}_i)$ are correctly specified, the estimating function $\sum_{i = 1}^n\psi_{p + 1}^i({{}\beta}, \phi)$ has expectation $(n - R)\phi$. So, assumption \[consistency\] on the unbiasedness of the estimating equations is satisfied for $R = n$. In other words, the estimating equations that result in $\hat\phi_n$ in (\[dispersion\]), rather than those for $\hat\phi_{n - p}$, are the ones that should be adjusted to produce the RB$M$-estimator of $\phi$.

According to Section \[sec:implementation\], the necessary quantities for the implementation of the empirical bias-reducing adjustment to the quasi-likelihood equations are the first and second derivatives of (\[quasilikelihood\]) with respect to ${{}\beta}$. If reduced-bias estimation of $\phi$ is also required, then the first and second derivatives of (\[quasilikelihood\]) with respect to $\phi$ and those of $m_i (Y_i - \mu_i)^2/v(\mu_i) - \phi$ with respect to ${{}\beta}$ and $\phi$ are also required. The closed-form expressions for ${{}j}^{(Q)}({{}\beta}, \phi)$, ${{}e}^{(Q)}({{}\beta}, \phi)$, ${{}d}_r^{(Q)}({{}\beta}, \phi)$ and $u_r^{(Q)}({{}\beta}, \phi)$ $(r = 1, \ldots, p + 1)$ are given in Appendix \[quasi\_appendix\], and are ready to be used for the implementation of the quasi Newton-Raphson iteration in Section \[sec:implementation\] for bias reduction in the estimation of both ${{}\beta}$ and $\phi$.

For bias reduction in the estimation of only ${{}\beta}$, the empirical bias-reducing adjustment is composed using only the $({{}\beta}, {{}\beta})$ blocks of ${{}j}^{(Q)}({{}\beta}, {{}\phi})$, ${{}e}^{(Q)}({{}\beta}, \phi)$, ${{}d}_r^{(Q)}({{}\beta}, \phi)$ and ${{}u}_r^{(Q)}({{}\beta}, \phi)$ $(r = 1, \ldots, p)$ in Appendix \[quasi\_appendix\]. Then, a direct computation shows that the empirical bias-reducing adjustment is inversely proportional to $\phi$, exactly as (\[quasilikelihood\]) is. Hence, the value of the dispersion parameter is not needed for getting the RB$M$-estimate of ${{}\beta}$ with the empirical bias-reducing adjustment in quasi likelihoods and generalized linear models. In fact, similarly to (\[pen\_loglik\_glm\]), the RB$M$-estimates of ${{}\beta}$ can be obtained by maximizing the bias-reducing penalized quasi log-likelihood $$\sum_{i = 1}^n m_i \left\{ \int_{y_i}^{\mu_i} \frac{y_i - t}{v(t)} d
    t - \frac{1}{2} s_i \frac{d_i}{v_i}(y_i - \mu_i) \right\} \, ,$$ with respect to ${{}\beta}$ which does not involve the dispersion parameter $\phi$. In contrast, @kosmidis+kennepagui+sartori:2019 show that the bias-reduction method in @firth:1993 requires the joint estimation of ${{}\beta}$ and $\phi$ in generalized linear models with an unknown dispersion parameter.

Over-dispersion in count responses {#sec:over_simu}
----------------------------------

In order to illustrate the impact of the empirical bias-reducing adjustment in quasi-likelihood estimation we assume that the responses $Y_1, \ldots, Y_n$ are generated independently from a negative binomial regression model with mean $\mu_i = \expect(Y_i | {{}x}_i) = \exp\{\beta_0 + \beta_1 x_{i1}+
\beta2 x_{i2}\}$ and variance $\var(Y_i | {{}x}_i) = \phi \mu_i$ $(i = 1, \ldots, n)$. This model is discussed in @mccullagh+nelder:1989 [Section 6.2.3] as a way to model count responses with dispersion that is greater than what is implied by a Poisson log-linear model. We generate two sets $\{x_{11}, \ldots, x_{n_01}\}$ and $\{x_{12}, \ldots, x_{n_02}\}$ of covariate values, independently and independent of each other from a Bernoulli distribution with probability $0.5$ and an exponential distribution with rate $2$, respectively. Conditionally on those covariate values, we simulate $N_0 = 4000$ samples of $n_0 = 20$ response values at $(\beta_0, \beta_1, \beta_2, \phi)^\top = (2, 1, -1, 6)$. Such a parameter setting implies that the conditional variances of the observed counts are 6 times greater than what would have been prescribed by a Poisson log-linear model. We also simulate $N_r= 2^r N_0$ samples of $n_r = 2^r n_0$ response values, after replicating each of the $n_0$ covariate settings $2^r$ times for $r \in \{1, \ldots, 7 \}$. This tuning of the simulation and sample sizes with $r$ guarantees that the simulation standard error when estimating the scaled bias ($n_r$ times the bias) of an estimator is $O(\sqrt{n_0} / \sqrt{N_0})$ and, hence, asymptotically bounded for any $r$. According to the theory in Section \[sec:theory\], the scaled bias for the quasi-likelihood estimator is $O(1)$ and for the reduced-bias estimator is $O(n^{-1/2})$. As a result, the bias of the reduced-bias estimator should converge to zero as $r$ increases, while that of the quasi-likelihood estimator should stabilize at a value that is not necessarily zero. Figure \[negbin\], shows the estimated scaled biases for various estimators for $\beta_0, \beta_1, \beta_2, \phi$ along with approximate $99\%$ Wald-type intervals for the scaled biases. The estimators for ${{}\beta}$ examined are the maximum likelihood estimator, the quasi-likelihood estimator, and the RB$M$-estimators in the estimation of only ${{}\beta}$, and of both ${{}\beta}$ and $\phi$. The quasi-likelihood and RB$M$-estimators of ${{}\beta}$ are examined not only for the correctly-specified variance function $V(\mu) = \mu$, but also for the incorrectly-specified variance function $V(\mu) = \mu^2$. As shown in the bottom right of Figure \[negbin\], there were very few cases where iteration (\[quasiNR\]), as we implemented it, did not converge (see the scripts for quasi likelihoods and negative binomial regression supplied in the Supporting Materials). These have been excluded from the computation of the scaled biases. The scaled biases for the estimators of $\phi$ are only provided for the correctly-specified model.

As expected from theory, the scaled biases of the maximum likelihood and quasi-likelihood estimators of ${{}\beta}$ converge to non-zero values. In contrast, the scaled biases of both versions of the RB$M$-estimator converge to zero, demonstrating that the empirical bias-reducing adjustment delivers estimators with bias of smaller order. When the variance function is correctly specified, the empirical bias-reducing adjustment delivers an estimator of $\phi$ with second-order bias. In contrast, the bias of $\hat\phi_{n - p}$ is of the same asymptotic order to that of the maximum likelihood estimator, and, hence, $\hat\phi_{n - p}$ does not have bias of smaller asymptotic order for arbitrary quasi likelihoods.

Figure \[negbin\_variance\] shows the variances of the estimators examined in Figure \[negbin\]. For correctly-specified variance function, the variances of the RB$M$-estimators are almost the same to the variances of the quasi-likelihood estimators for the sample sizes considered. Also, the variance of the RB$M$-estimator of $\phi$ is almost the same to that of $\hat\phi_{n - p}$. When the variance function is misspecified, the RB$M$-estimators do considerably better in terms of finite sample variance than the quasi-likelihood estimator for $n_0 = 20$. As expected from Section \[sec:inference\], the differences in the variances of all estimators vanish as the sample size increases.

Gaussian max-stable processes {#sec:maxstable}
=============================

Extreme climate events are becoming more frequent all over the world, with strong impact on the built environment. As a result, their statistical modelling and prediction has attracted a lot of research. Vanilla likelihood and Bayesian approaches for spatial extreme processes face challenges, because the direct generalization of the classical multivariate extreme value distributions to the spatial case is a max-stable process for which the evaluation of the likelihood becomes increasingly more intractable as the number of site locations increases [@davison+gholamrezaee:2011].

Several works have proposed the use of computationally appealing surrogates to the likelihood, like composite likelihoods, which are formed by specifying marginal or conditional densities for subsets of site locations [see @padoan+ribadet+sisson:2010; @genton+ma+sang:2011; @davison+gholamrezaee:2011; @huser+davison:2013 among others]. Nevertheless, standard bias-reduction methods, like the one in @firth:1993 and @kuk:1995 in Table \[br\_characteristics\], are either infeasible or computationally expensive because the calculation of the bias function involves either integrals with respect to the true underlying joint density or requires resampling and refitting.

An example of such a surrogate to the likelihood function is the pairwise likelihood introduced in @padoan+ribadet+sisson:2010 under a block maxima approach to the modelling of extremes. Suppose that $y_1({}s),\ldots,y_k({}s)$ with ${}s \in \{{}s_1,\ldots,{}s_L\}$, ${}s_j \in \Re^2$, are $k$ independent observations at $L$ site locations. The pairwise log-likelihood formed from the collection of $L(L-1)/2$ distinct pairs of locations is $$\label{pwl_mev}
l({{}\theta})=\sum_{i=1}^k\sum_{l>m} \log f(y_i({}s_l),y_i({}s_{m}) | {{}\theta}) \,,$$ where $$\begin{aligned}
  \label{BivDensMev}
  f(y_i({}s_l),y_i({}s_{m})|{{}\theta})  = \exp\left\{-\frac{\Phi(w_{lm})}{y_i({}s_l)}-\frac{\Phi(v_{lm})}{y_i({}s_{m})}\right\} \left[ \frac{v_{lm}\phi(w_{lm})}{a^2_{lm}y^2_i({}s_l)y_i({}s_{m})} + \frac{w_{lm}\phi(v_{lm})}{a^2_{lm}y^2_i({}s_{m})y_i({}s_l)} \right. & \\ \notag
  + \left. \left\{ \frac{\Phi(w_{lm})}{y^2_i({}s_l)} + \frac{\phi(w_{lm})}{a_{lm}y^2_i({}s_l)} - \frac{\phi(v_{lm})}{a_{lm}y_i({}s_l)y_i({}s_{m})} \right\}
  \left\{ \frac{\Phi(v_{lm})}{y^2_i({}s_{m})} + \frac{\phi(v_{lm})}{a_{lm}y^2_i({}s_{m})} - \frac{\phi(w_{lm})}{a_{lm}y_i({}s_l)y_i({}s_{m})} \right\} \right],\end{aligned}$$ is the joint density of $Y_i({}s_l)$ and $Y_i({}s_{m})$ $(l, m = 1, \ldots, L; l \ne m)$, with $\Phi(\cdot)$ and $\phi(\cdot)$ the distribution and density function of the standard normal distribution, respectively. In the above expression, $a_{lm} \equiv a_{lm}({{}\theta}) = \{({}s_l-{}s_m)^\top{{}\Sigma}^{-1}({{}\theta})({}s_l-{}s_m)\}^{1/2}$, $w_{lm}\equiv w_{lm}({{}\theta}) = a_{lm}/2+\log\{y_i({}s_l)/y_i({}s_{m})\}/a_{lm}$, and $v_{lm} = a_{lm}-w_{lm}$. The parameter vector is ${{}\theta}=(\sigma^2_{1}, \sigma^2_{2}, \sigma^2_{12})^\top$, which collects the distinct elements of the $2\times 2$ covariance matrix ${{}\Sigma}({{}\theta})$ that governs the spatial dependence, and has diagonal elements $\sigma^2_1$ and $\sigma^2_2$, with $\sigma_{12}^2$ in the off-diagonals.

The maximizer of is the maximum pairwise likelihood estimator $\hat{{}\theta}$, while the RB$M$-estimator $\tilde{{}\theta}$ results by maximizing the penalized version of $\l({}\theta)$ in (\[penalized\_likelihood\]). Closed-form expressions for ${{}j}({{}\theta})$ and ${{}e}({{}\theta})$ are given in Section S2 of the Supporting Materials.

Simulations are run by generating independent observations $y_1({}s),\ldots,y_k({}s)$ from a Gaussian max-stable process observed at $L=50$ site locations. The locations are generated uniformly on a $[0, 40] \times [0, 40]$ region. We consider sample sizes $k \in \{10, 20, 40, 80, 160\}$ with corresponding number of simulations equal to $4000 k$, and true parameter settings ${{}\theta} = (2000, 3000, 1500)^\top$ and ${{}\theta} = (20, 30, 15)^\top$, imposing strong and weak spatial dependence, respectively. These parameter values correspond to ${{}\Sigma}_4$ and ${{}\Sigma}_5$ in Table 1 of @padoan+ribadet+sisson:2010. In Figure \[MEV\_plots\] (top panel), we show the simulation-based estimates of $\log |\expect_G\{\hat{{}\theta}-\bar{{}\theta}\}|$ and $\log |\expect_G\{\tilde{{}\theta}-\bar{{}\theta}\}|$ as functions of $\log n$, where $n = k$. These curves have roughly slopes $-1$, and between $-3/2$ and $-2$, respectively as expected by the asymptotic theory in Section \[sec:theory\], demonstrating the reduction of the bias that $\tilde{{}\theta}$ delivers. Furthermore, $\tilde{{}\theta}$ appears to have smaller finite-sample mean squared error, and hence smaller variance than the maximum pairwise likelihood estimator $\hat{{}\theta}$; see bottom panel in Figure \[MEV\_plots\]. The simulation results provide evidence for the superiority of the RB$M$-estimator. The differences between $\hat{{}\theta}$ and $\tilde{{}\theta}$ will tend to be significant in the small to moderate samples that are typically observed in settings involving block-maxima.

[![Results of the simulation study for Gaussian max-stable process: right and left panels refer to $\Sigma_4$ and $\Sigma_5$, respectively. Summaries for $\hat{{}\theta}$ and $\tilde{{}\theta}$ are depicted in black and gray, respectively. The top panels show simulation-based estimates of the logarithm of the absolute bias of $\hat{{}\theta}$ and $\tilde{{}\theta}$ (solid). The dashed line is a line with slope $-1$ corresponding to the theoretical rate of the bias of $\hat{{}\theta}$. The shadowed region is defined by lines with slopes $-3/2$ and $-2$, which correspond to the theoretical rate of the bias of $\tilde{{}\theta}$. The bottom panels show the simulation-based estimates of the root mean squared error of $\hat{{}\theta}$ and $\tilde{{}\theta}$.[]{data-label="MEV_plots"}](Figures/RATE_Sigma4 "fig:")]{} [![Results of the simulation study for Gaussian max-stable process: right and left panels refer to $\Sigma_4$ and $\Sigma_5$, respectively. Summaries for $\hat{{}\theta}$ and $\tilde{{}\theta}$ are depicted in black and gray, respectively. The top panels show simulation-based estimates of the logarithm of the absolute bias of $\hat{{}\theta}$ and $\tilde{{}\theta}$ (solid). The dashed line is a line with slope $-1$ corresponding to the theoretical rate of the bias of $\hat{{}\theta}$. The shadowed region is defined by lines with slopes $-3/2$ and $-2$, which correspond to the theoretical rate of the bias of $\tilde{{}\theta}$. The bottom panels show the simulation-based estimates of the root mean squared error of $\hat{{}\theta}$ and $\tilde{{}\theta}$.[]{data-label="MEV_plots"}](Figures/RATE_Sigma5 "fig:")]{} [![Results of the simulation study for Gaussian max-stable process: right and left panels refer to $\Sigma_4$ and $\Sigma_5$, respectively. Summaries for $\hat{{}\theta}$ and $\tilde{{}\theta}$ are depicted in black and gray, respectively. The top panels show simulation-based estimates of the logarithm of the absolute bias of $\hat{{}\theta}$ and $\tilde{{}\theta}$ (solid). The dashed line is a line with slope $-1$ corresponding to the theoretical rate of the bias of $\hat{{}\theta}$. The shadowed region is defined by lines with slopes $-3/2$ and $-2$, which correspond to the theoretical rate of the bias of $\tilde{{}\theta}$. The bottom panels show the simulation-based estimates of the root mean squared error of $\hat{{}\theta}$ and $\tilde{{}\theta}$.[]{data-label="MEV_plots"}](Figures/VAR_Sigma4 "fig:")]{} [![Results of the simulation study for Gaussian max-stable process: right and left panels refer to $\Sigma_4$ and $\Sigma_5$, respectively. Summaries for $\hat{{}\theta}$ and $\tilde{{}\theta}$ are depicted in black and gray, respectively. The top panels show simulation-based estimates of the logarithm of the absolute bias of $\hat{{}\theta}$ and $\tilde{{}\theta}$ (solid). The dashed line is a line with slope $-1$ corresponding to the theoretical rate of the bias of $\hat{{}\theta}$. The shadowed region is defined by lines with slopes $-3/2$ and $-2$, which correspond to the theoretical rate of the bias of $\tilde{{}\theta}$. The bottom panels show the simulation-based estimates of the root mean squared error of $\hat{{}\theta}$ and $\tilde{{}\theta}$.[]{data-label="MEV_plots"}](Figures/VAR_Sigma5 "fig:")]{}

Migraine severity study {#autoregressive_probit}
=======================

We consider the data from a longitudinal study on the determinants of migraine severity, as has been analyzed in @varin+czado:2010. The study involved $n = 133$ Canadian patients that were asked to rate headache pain levels on an ordinal scale and record them at four pre-specified daily occasions. Potential explanatory variables, such as subject specific information and meteorological data, were collected with the pain ratings.

Suitable regression models for longitudinal ordinal responses that account for subject-specific variability are probit models with random effects. Let $Y_{ij}$ be the pain rating by the $i$th individual at time $t_{ij}$, and let ${{}x}_{ij}$ be the $r$-dimensional vector containing explanatory variables with $x_{i1} = 1$ $(i = 1, \ldots, n; j = 1, \ldots, n_i)$. Suppose that $Y_{ij}$ takes values from $\{1, \ldots, h\}$ and that it results from the discretization of a continuous random variable $Z_{ij}$ according to the rule $Y_{ij} = y_{ij} \iff \alpha_{y_{ij}-1} < Z_{ij} \leq
\alpha_{y_{ij}}$, where ${{}\alpha}=(\alpha_0,\alpha_1,\ldots,\alpha_h)^\top$ is a vector of cut-points with $\alpha_0<\alpha_1 \ldots <\alpha_h$, $\alpha_0 = -\infty$, $\alpha_1 = 0$, and $\alpha_h = \infty$, where $Z_{ij} = {{}x}^\top_{ij}{{}\beta} + U_i + \epsilon_{ij}$. The $r$-dimensional vector ${{}\beta}$ collects the regression parameters, $U_i$ are random effects, and $\epsilon_{ij}$ are random errors. Under normality assumptions on both the random effects and errors, the joint distribution of $(Z_{i1}, \ldots, Z_{in_i})$ is multivariate normal. The computation of the log-likelihood function about ${{}\alpha}$, ${{}\beta}$ and the parameters of the distribution of the random effects and errors is typically impractical because it requires the approximation of $n$ integrals, with the $i$th integral being of dimension $n_i$.

@varin+czado:2010 have overcome those computational challenges by resorting to a pairwise likelihood whose specification requires the joint distribution for the pairs $(Y_{ij}, Y_{ik})$ $(j, k = 1, \ldots, n_i; j\neq k)$ and, in turn, the numerical approximation of bivariate integrals only.

Under the assumption that $U_i \sim N(0, \sigma^2)$ and $\epsilon_{ij} \sim N(0, 1)$, with $\cov(U_i, U_j)=0$, $\cov(\epsilon_{ik}, \epsilon_{il}) = 0$, and $\cov(U_{i}, \epsilon_{ik}) = 0$, the pairwise log-likelihood about ${{}\theta}=({{}\alpha}^\top, {{}\beta}^\top, {{}\omega}^\top,
\sigma^2)^\top$ that is introduced in @varin+czado:2010 is $$\label{pllik_varin_czado} 
  l({{}\theta}) = \sum_{i=1}^n\sum_{j>k}^{n_i}\log P(Y_{ij}=y_{ij}, Y_{ik}=y_{ik} ; {{}\theta}) \mathbb{I}_{[-q,q]}(t_{ij} - t_{ik}),$$ where $\mathbb{I}_{[-q,q]}(a)$ takes value $1$ if $a \in [-q, q]$ and $0$, otherwise, and $$\label{int_varin_czado}
  P(Y_{ij}=y_{ij}, Y_{ik}=y_{ik}; {{}\theta})=\int^{\tilde\alpha_{y_{ij}}}_{\tilde\alpha_{y_{ij}-1}}\int^{\tilde\alpha_{y_{ik}}}_{\tilde\alpha_{y_{ik}-1}}\phi_2(u,v; \rho_{ijk})du dv \, ,$$ with $\tilde\alpha_{y_{ij}} = (\alpha_{y_{ij}} - {{}x}^\top_{ij}
{{}\beta})(1 + \sigma^2)^{-1/2}$ and $\phi_2(\cdot; \rho_{ijk})$ denoting the bivariate normal density function with zero means and correlation $\rho_{ijk} = (1 + \sigma^2)^{-1}\omega_i^{|t_{ij}-t_{ik}|}$, allowing for subject heterogeneity in the correlation functions. Computational speed and statistical efficiency gains [see, for example @joe+lee:2009; @bevilacqua+gaetan+mateu+etal:2012] can be achieved by discarding too distant pairs through a suitable choice of $q$.

@varin+czado:2010 analyzed the migraine data by using and derived the associated inferential tools. We repeat their analysis using the bias-reducing penalized version of and compare the results; the derivation and closed-form expressions for the matrices ${}e({}\theta)$ and ${}j({}\theta)$ in (\[penalized\_likelihood\]) are given in Section S3 of the Supporting Materials.

The five explanatory variables we consider are university degree status, analgesics intake, change of atmospheric pressure between consecutive days, humidity, and windchill. A brief description of those variables and the response is in Table \[table:migraine\]. We set a base model that has the university degree status and analgesics intake as explanatory variables, and use the CLIC at the maximum pairwise likelihood and the RB$M$-estimates to select the best out of the $8$ possible models that result by including or not the possible combinations of the remaining $3$ explanatory variables. As in @varin+czado:2010, we use $q = 12$ and set the autocorrelation structure in all models to be $\rho_{ijk} = (1 + \sigma^2)^{-1}\omega_u^{|t_{ij}-t_{ik}|}$, with $u = 0$ if the $i$th patient does not make use of analgesics, and $u = 1$, otherwise.

We define the bias-reducing penalized pairwise log-likelihood for ${{}\xi}={{}\xi({}\theta)} = ({{}\gamma}^\top, {{}\beta}^\top,
{{}\tau}^\top, \lambda)^\top$, where $\tau_u = \log\{\omega_u/(1 - \omega_u)\}$ $(u = 0, 1)$, $\lambda = \log\sigma^2$, $\gamma_0=0$, and $\gamma_m = \log(\alpha_{m+1}-\alpha_{m})$ $(m = 1, \ldots, h-1)$. In this way, all parameters take values on the real line avoiding numerical issues during maximization. By the equivariance properties of the maximum pairwise likelihood estimator, $\hat{{}\xi}={}\xi(\hat{{}\theta})$, where $\hat{{}\theta}$ is the maximizer of . Note here that back-transforming the RB$M$-estimators ${{}\gamma}$, ${{}\tau}$ and $\lambda$ does not result in RB$M$ estimators for ${{}\alpha}$, ${{}\omega}$, $\sigma^2$ because the bias-reducing penalty is not parametrization invariant under non-linear transformations of the parameters. Nevertheless, the RB$M$-estimators of the parameters of interest ${{}\beta}$ will have improved bias in either parameterization.

Table \[table:migraine clic\] reports the CLIC information criteria for models estimated using and its penalized version. Table \[table:migraine clic\] also reports the CLIC weights, defined as $w_k = \exp(-\Delta_k)/\sum_{i=1}^8 \exp(-\Delta_i)$, $\Delta_k=(\text{CLIC}_k-\min_c \text{CLIC}_c)/2$. According to the CLIC weights for models estimated by maximum pairwise likelihood, the best and second best models are the base model with change of atmospheric pressure, and the base model, respectively. We arrive at the same conclusion by using the CLIC weights at the RB$M$-estimates, but there is stronger evidence in favour of the base model with the CLIC weight increasing from $0.28$ to $0.38$.

  ------------------------------------------ ----------------------------------------------------------------------------------------------------------------------------------------------------------
  Variable                                                                                                           Categories
  \[.5ex\] Headache pain rating (response)                                                            no/mild/moderate/painful/severe/intense
  University degree                                                                                                    yes/no
  Analgesics intake                                                                                                    yes/no
  Change of atmospheric pressure                                                                         high to low/low to high/unchanged
  Humidity                                                                                                     $<$60%/60%-80%/$>$80%
  Windchill                                   $(-50^\circ\text{C},-10^\circ\text{C})$/$(-10^\circ\text{C},0^\circ\text{C})$/$(0^\circ\text{C},10^\circ\text{C})$/$(10^\circ\text{C},30^\circ\text{C})$
  ------------------------------------------ ----------------------------------------------------------------------------------------------------------------------------------------------------------

  : Variables in the migraine data. The response variable was recorded on an ordinal scale with six categories whose full length description may be found in @varin+czado:2010 [Table 1]. The correspondence between the levels of the variable referring to the change of atmospheric pressure and the actual change in pressure are from high ($>$1013 hPa) to low pressure ($\leq$ 1013 hPa), from low to high pressure, and unchanged level of pressure (either from low to low or from high to high).

\[table:migraine\]

  -- -------- ---------- ----------- --------- ------- --------- -------
                                                                 
                                                                 
      Change   Humidity   Windchill    CLIC     $w_k$    CLIC     $w_k$
        -         -           -       5919.17   0.28    5940.74   0.38
        \*        -           -       5918.26   0.44    5940.40   0.45
        -         \*          -       5921.11   0.11    5943.80   0.08
        -         -          \*       5925.00   0.02    5949.91   0.00
        \*        \*          -       5920.60   0.14    5943.87   0.08
        \*        -          \*       5924.45   0.02    5950.01   0.00
        -         \*         \*       5927.40   0.00    5953.45   0.00
        \*        \*         \*       5927.09   0.01    5953.80   0.00
  -- -------- ---------- ----------- --------- ------- --------- -------

  : CLIC statistics and CLIC weights for the models considered for the migraine data set, with $q=12$. CLIC for pairwise likelihood refers to , whereas CLIC for the bias-reducing penalized pairwise likelihood refers to once multiplied by $-2$. The symbol \* indicates that a variable is included in the model and $w_k$ refers to the CLIC weight of each model. CLIC values for pairwise likelihood displayed here differ slightly from the ones in @varin+czado:2010 [Table 4] because in our analysis the penalty term in the CLIC has been evaluated using the actual hessian of rather than the approximation used in @varin+czado:2010 [equation 3.1].

\[table:migraine clic\]

Estimates from the best two models according to CLIC weights are displayed in Table \[table:migraine est\]. While the maximum pairwise likelihood estimator is equivariant, the corresponding estimated standard errors are not, so that hypotheses tests based on Wald-type statistics for the main effects refer to the ${{}\xi}$ parametrization. However, the differences between the estimated standard errors for ${{}\beta}$ in the two parametrizations are small and the conclusions are the same as the ones stated in @varin+czado:2010. Furthermore, the difference between the maximum pairwise likelihood and RB$M$-estimates of ${{}\beta}$ are small, indicating that estimation bias has not been a major concern in that case study. In particular, analgesics intake is related to higher pain level ratings, the university degree acts in the opposite direction, and the decrease of atmospheric pressure is associated with worse pain ratings. Overall, the estimated standard errors for the RB$M$-estimator $\tilde{{}\xi}$ are only slightly smaller than those for $\hat{{}\xi}$.

  ---------------------- ---------- ------- ---------- ------- -- ---------- ------- ---------- -------
  $\gamma_2$               $-$0.530   0.088   $-$0.529   0.088      $-$0.534   0.088   $-$0.533   0.088
  $\gamma_3$               $-$0.602   0.069   $-$0.601   0.069      $-$0.602   0.069   $-$0.602   0.069
  $\gamma_4$               $-$0.430   0.061   $-$0.429   0.061      $-$0.428   0.060   $-$0.427   0.061
  $\gamma_5$               $-$0.330   0.092   $-$0.330   0.092      $-$0.329   0.092   $-$0.329   0.093
  Intercept                $-$0.474   0.230   $-$0.521   0.230      $-$0.461   0.230   $-$0.510   0.231
  University               $-$0.523   0.156   $-$0.524   0.158      $-$0.537   0.156   $-$0.537   0.158
  Analgesics                  0.558   0.221      0.560   0.223         0.544   0.220      0.546   0.222
  Change - unchanged                             0.031   0.046                            0.033   0.046
  Change - high to low                           0.164   0.051                            0.164   0.052
  $\tau_0$                    0.886   0.322      0.891   0.329         0.915   0.306      0.921   0.312
  $\tau_1$                    1.253   0.110      1.252   0.110         1.250   0.106      1.249   0.106
  $\tau_1-\tau_0$             0.368   0.347      0.361   0.353         0.335   0.334      0.327   0.340
  $\lambda$                $-$0.567   0.196   $-$0.566   0.198      $-$0.541   0.189   $-$0.541   0.191
  ---------------------- ---------- ------- ---------- ------- -- ---------- ------- ---------- -------

  : Estimates and estimated standard errors (s.e.) from maximum pairwise likelihood ($\hat{{}\xi}$) and from maximum bias-reducing penalized pairwise likelihood ($\tilde{{}\xi}$) for the best two models, according to the CLIC values in Table \[table:migraine clic\]. Columns 1-2 and 5-6 refer to base models, whereas 3-4 and 7-8 to the models including the change of atmospheric pressure between consecutive days. The baseline categories for university degree status and analgesics intake is “no”, and for the variable change of atmospheric pressure between consecutive days is “low to high”.

\[table:migraine est\]

Discussion and further work {#sec:discussion}
===========================

We have developed a novel, general framework for the asymptotic reduction of the bias of general $M$-estimators from sufficiently smooth, unbiased estimating functions. Bias reduction is achieved by the adjustment of the estimating functions by quantities that are bounded in probability, and depend only on the first and second derivatives of contributions to the estimating functions. The RB$M$-estimates can be computed using the general quasi Newton-Raphson iteration in Section \[sec:implementation\] that again requires only the contributions to the estimating functions and the first two derivatives of those.

The resulting estimators have the same asymptotic distribution, and, hence, they are asymptotically as efficient as the original, unadjusted $M$-estimators. As detailed in Section \[sec:inference\], uncertainty quantification can be carried out using the empirical estimate $\hat{{{}V}}({{}\theta})$ of the variance-covariance matrix of that asymptotic distribution. The expression for $\hat{{{}V}}({{}\theta})$ appears in expression (\[general\_empirical\_adjustment\_matrix\]) for the empirical adjustment, and, hence, is readily available at the last iteration of the quasi Newton-Raphson fitting procedure. Inferences can be constructed using the Wald-type and generalized score pivots in expression (\[pivots\]).

An alternative estimator with $o(n^{-1})$ bias in general $M$-estimation problems under the assumptions of Section \[assumptions\] is $$\label{bias_correction}
  {{}\theta}^\dagger = \hat{{}\theta} + {{}j}(\hat{{}\theta})^{-1} {{}A}(\hat{{}\theta}) \, ,$$ where the $r$th component of ${{}A}({{}\theta})$ is as in expression (\[general\_empirical\_adjustment\_matrix\]) for the empirical bias-reducing adjustments. According to the classification of bias-reduction methods in @kosmidis:2014, the estimator (\[bias\_correction\]) defines an explicit bias-reduction method. The value of ${{}\theta}^\dagger$ in (\[bias\_correction\]) can be computed with a single step of the quasi Newton-Raphson procedure of Section \[sec:implementation\] with $a_1 = 1$, starting at the $M$-estimator $\hat{{}\theta}$. The choice between using ${{}\theta}^\dagger$ and the fully iterated version $\tilde{{}\theta}$ that has been introduced earlier in this paper is application-dependent. For example, for the estimation of a ratio of two means in Example \[ratio\_estimation\], $\theta^\dagger = \hat\theta (1 - s_{XX} / s_{X}^2) + s_{XY}/ s_X^2$, which, like $\hat\theta = s_Y/s_X$ and unlike the fully-iterated $\tilde\theta$, is not robust to small values of $s_X$. Furthermore, if $M$-estimation is through the maximization of an objective function, then $\tilde{{}\theta}$ is the maximizer of the penalized objective (\[penalized\_likelihood\]). As discussed in Section \[sec:implementation\], the maximization of (\[penalized\_likelihood\]) can be performed using a general numerical optimization routine that operates by numerically approximating the gradient of (\[penalized\_likelihood\]) and requires only the estimating function contributions and their first derivatives. In contrast, by its definition, ${{}\theta}^\dagger$ requires also the second derivatives of the estimating function contributions.

A derivation similar to that in Section \[derivation\] and Section \[sec:penalty\] can be used to show that another estimator with $o(n^{-1})$ bias results by the maximization of the bias-reducing penalized objective $$\label{penalized_likelihood2}
  l({{}\theta}) + \frac{1}{2} \log \det \left\{{{}j}({{}\theta}) \right\} - \frac{1}{2} \log \det\left\{ {{}e}({{}\theta}) \right\} \, .$$ The resulting estimator has, again, the same asymptotic distribution to the original $M$-estimators (see Section \[sec:inference\]), however, it is, typically, less attractive than $\tilde{{}\theta}$ from a computational point of view, because of the need to compute the logarithms of two determinants.

We also showed that the bias-reducing penalized objective (\[penalized\_likelihood\]) closely relates to model selection procedures using the Kullback-Leibler divergence. Such model selection procedures include the use of the TIC for maximum likelihood estimation, and the CLIC for maximum composite likelihood estimation. The functions of the parameters and the data that are used for bias reduction and model selection differ only by a known scalar constant. These results establish, for the first time, a close relation between bias reduction in estimation and model selection.

Furthermore, TIC and CLIC are still consistent information criteria when evaluated at the RB$M$-estimates, and, hence, their value is readily available once the bias-reducing penalized likelihood has been maximized. The only difference with the standard versions of TIC and CLIC is that when evaluated at the RB$M$-estimates these criteria are only asymptotically invariant to non-linear transformations of the parameters. The same justification we provided in Section \[sec:model\_selection\] for the use of information criteria at the reduced-bias estimates can be used to justify the use of information criteria at estimates arising from the additive adjustment of estimating functions by alternative $O_p(1)$ quantities, like the median reduced-bias estimates discussed in @kenne+salvan+sartori:2017 and @kosmidis+kennepagui+sartori:2019, and the mean reduced-bias estimators in @firth:1993.

If there is only $k = 1$ multivariate observation from the underlying process, application-dependent conditions need to be used for the appropriate definition of ${{}e}({{}\theta})$ in the penalized objective function in (\[penalized\_likelihood\]) or of ${{}e}({{}\theta})$ and ${{}d}_r({{}\theta})$ in the empirical bias-reducing adjustment (\[general\_empirical\_adjustment\_matrix\]). For example, in the context of time series and spatial data that do not seriously depart from the condition of stationarity, one can consider window sub-sampling for the definition of ${{}e}({{}\theta})$ and ${{}d}_r({{}\theta})$ [see, for example @carlstein:1986; @heagerty+lumley:2000 for definitions and guidance on the choice of the window size].

@lunardon:2018 showed that bias reduction in maximum likelihood estimation using the adjustments in @firth:1993 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. Current research investigates the performance of the reduced-bias estimator from empirically adjusted estimating functions for general $M$-estimation in stratified settings, extending the arguments and optimality results in @lunardon:2018 when maximum composite likelihood and other $M$-estimators are employed.

Supporting materials
====================

The supporting materials include computer code to fully reproduce all numerical results and figures in the paper. The organization of the computer code is detailed in Section S1 of the supporting materials document. Section S2 and Section S3 of that document derive the mathematical expressions that are used in Section \[sec:maxstable\] and Section \[autoregressive\_probit\], respectively.

Acknowledgements
================

Ioannis Kosmidis is supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1.

Appendix
========

For notational simplicity, the dependence of the various quantities below on ${{}\beta}$ and/or $\phi$ is suppressed.

Expressions for the bias-reducing penalty for generalized linear models {#glm_appendix}
-----------------------------------------------------------------------

$${{}j}^{(G)} = \left[
              \begin{array}{cc}
                {{}j}_{{{}\beta}{{}\beta}} & {{}j}_{{{}\beta}\phi} \\
                {{}j}_{{{}\beta}\phi}^\top &  {{}j}_{\phi\phi}^{(G)}  \\
              \end{array}
              \right] \quad \text{and} \quad
  {{}e}^{(G)} = \left[ 
              \begin{array}{cc}
                {{}e}_{{{}\beta}{{}\beta}} & {{}e}_{{{}\beta}\phi}^{(\text{G})} \\
                \left\{{{}e}_{{{}\beta}\phi}^{(\text{G})}\right\}^\top &  {{}e}_{\phi\phi}^{(G)}  \\
              \end{array}
              \right] \, ,$$

where $$\begin{aligned}
  {{}j}_{{{}\beta}{{}\beta}} & = \frac{1}{\phi} {{}X}^\top {{}Q} {{}X} \, , &
  {{}j}_{\phi\phi}^{(G)} & = \frac{1}{\phi^3} {{}1}_n^\top ({{}R} - {{}A}') {{}1}_n + \frac{1}{2\phi^4} {{}1}_n^\top {{}A}'' {{}1}_n \, , &
  {{}j}_{{{}\beta}\phi} & = \frac{1}{\phi^2} {{}X}^\top \tilde{{{}W}} {{}1}_n \, , \\
  {{}e}_{{{}\beta}{{}\beta}} & = \frac{1}{\phi^2} {{}X}^\top \tilde{{{}W}}^2 {{}X} \, , &
  {{}e}_{\phi\phi}^{(G)} & = \frac{1}{4 \phi^4} {{}1}_n^\top ({{}R} - {{}A}')^2 {{}1}_n \, , &
  {{}e}_{{{}\beta}\phi}^{(G)} & = \frac{1}{2\phi^3} {{}X}^\top \tilde{{{}W}} ({{}R} - {{}A}') {{}1}_n \, .\end{aligned}$$ In the above expressions, ${{}1}_n$ is an $n$-vector of ones. The $n \times n$ diagonal matrices ${{}Q}$, $\tilde{{{}W}}$, ${{}R}$, ${{}A}'$, ${{}A}''$ have $i$th diagonal element $q_i = b_i d_i - b_i' (y_i - \mu_i)$, $\tilde{g}_i = b_i (y_i - \mu_i)$, $r_i = -2 m_i (y_i \theta_i - b(\theta_i) - c_1(y_i))$ (deviance residual), $a_i' = m_i a'(-m_1/\phi)$, $a_i'' = m_i^2 a''(-m_1/\phi)$, respectively, where $b_i = m_i d_i / v_i$, $b_i' = m_i (d_i'/v_i - d_i^2 v_i'/v_i^2)$, $d_i' = d^2 \mu_i / d\eta_i^2$, $v_i' = d v_i /d \mu_i$ $(i = 1, \ldots, n)$, and $a'(u) = d a(u) / d u$, $a''(u) = d^2 a(u) / d u^2$.

Expressions for the empirical bias-reducing adjustment for quasi likelihoods {#quasi_appendix}
----------------------------------------------------------------------------

$${{}j} = \left[
    \begin{array}{cc}
      {{}j}_{{{}\beta}{{}\beta}} & {{}j}_{{{}\beta}\phi} \\
      {{}j}_{\phi{{}\beta}}^{(\text{Q})} &  {{}j}_{\phi\phi}^{(Q)}  \\
    \end{array}
              \right] \quad \text{and} \quad
  {{}e} = \left[ 
              \begin{array}{cc}
                {{}e}_{{{}\beta}{{}\beta}} & {{}e}_{{{}\beta}\phi}^{(\text{Q})} \\
                \left\{{{}e}_{{{}\beta}\phi}^{(\text{Q})}\right\}^\top &  {{}e}_{\phi\phi}^{(Q)}  \\
              \end{array}
              \right] \, ,$$

where $$\begin{aligned}
  {{}j}_{\phi\phi}^{(Q)} & = n \, , &
  {{}j}_{\phi{{}\beta}}^{(Q)} & = {{}1}_n^\top {{}F} {{}X} \, , \\
  {{}e}_{\phi\phi}^{(Q)} & = {{}1}_n^\top ({{}K} - \phi {{}I}_n)^2 {{}1}_n \, , &
  {{}e}_{{{}\beta}\phi}^{(Q)} & = \frac{1}{\phi} {{}X}^\top \tilde{{{}W}} ({{}K} - \phi {{}I}_n) {{}1}_n \, .\end{aligned}$$ In the above expressions, ${{}I}_n$ is the $n \times n$ identity matrix. The $n \times n$ diagonal matrices ${{}F}$ and ${{}K}$ have $i$th diagonal element $f_i = 2 c_i d_i (y_i - \mu_i) - c_i' (y_i - \mu_i)^2$ and $k_i = c_i (y_i - \mu_i)^2$, respectively, with $c_i = m_i / v_i$ and $c_i' = - m_i v_i' d_i / v_i^2$ $(i = 1, \ldots, n)$.

In addition, $${{}u}_r({{}\beta}, \phi) = \left[
    \begin{array}{cc}
      {{}u}_{r, {{}\beta}{{}\beta}} & {{}u}_{r, {{}\beta}\phi} \\
      {{}u}_{r, {{}\beta}\phi}^\top &  {{}u}_{r, \phi\phi}  \\
    \end{array}
  \right] \quad \text{and} \quad
  {{}u}_{p+1}({{}\beta}, \phi) = \left[ 
    \begin{array}{cc}
      {{}u}_{p+1, {{}\beta}{{}\beta}} & {{}u}_{p+1, {{}\beta}\phi} \\
      {{}u}_{p+1, {{}\beta}\phi}^\top &  {{}u}_{p+1, \phi\phi}  \\      
    \end{array}
  \right] \quad \quad (r = 1, \ldots, p) \, ,$$ where $$\begin{aligned}
  {{}u}_{r, {{}\beta}{{}\beta}} & = -\frac{1}{\phi} {{}X}^\top {{}Q}' {{}T}_r {{}X} \, &
  {{}u}_{r, \phi\phi} & = \frac{2}{\phi^3} {{}1}_n^\top {{}T}_r \tilde{{{}W}} {{}1}_n \, , &
  {{}u}_{r, {{}\beta}\phi} & = \frac{1}{\phi^2} {{}X}^\top {{}Q} {{}T}_r {{}1}_n \, , \\
  {{}u}_{p + 1, {{}\beta}{{}\beta}} & = {{}X}^\top {{}S} {{}X} \,, &
  {{}u}_{p + 1, \phi\phi} & = 0 \, , &
  {{}u}_{p + 1, {{}\beta}\phi} & = {{}0}_p \, . \\  \end{aligned}$$ The $n \times n$ diagonal matrices ${{}Q}'$, ${{}S}$, and ${{}T}_r$ have $i$th diagonal element $q_i' = 2 b_i' d_i + b_i d_i' - b_i'' (y_i - \mu_i)$, $s_i = c_i'' (y_i - \mu_i)^2 - 4 c_i' d_i (y_i - \mu_i) - 2 c_i d_i'
(y_i - \mu_i) + 2 c_i d_i^2$, and $x_{ir}$, respectively $(i = 1, \ldots, n)$, with $c_i'' = - m_i \left\{ v_i'' d_i^2 / v_i^2 - v_i' d_i'/v_i^2 - 2 d_i^2
  (v_i')^2 / v_i^3 \right\}$ and $b_i'' = m_i (d_i'' v_i - 3 d_i d_i' v_i' - d_i^3 v_i'')/v_i^2 + m_i
d_i^3 (v_i')^2/v_i^3$.

Finally, $${{}d}_r({{}\beta}, \phi) = \left[
    \begin{array}{cc}
      {{}d}_{r, {{}\beta}{{}\beta}} & {{}d}_{r, {{}\beta}\phi} \\
      {{}d}_{r, \phi{{}\beta}} &  {{}d}_{r, \phi\phi}  \\
    \end{array}
  \right] \quad \text{and} \quad
  {{}d}_{p+1}({{}\beta}, \phi) = \left[ 
    \begin{array}{cc}
      {{}d}_{p+1, {{}\beta}{{}\beta}} & {{}d}_{p+1, {{}\beta}\phi} \\
      {{}d}_{p+1, \phi{{}\beta}}^\top &  {{}d}_{p+1, \phi\phi}  \\      
    \end{array}
  \right] \quad \quad (r = 1, \ldots, p) \, ,$$ where $$\begin{aligned}
  {{}d}_{r, {{}\beta}{{}\beta}} & = \frac{1}{\phi^2} {{}X}^\top {{}Q} {{}T}_r \tilde{{{}W}} {{}X} \, ,&
  {{}d}_{r, \phi\phi} & = \frac{1}{\phi^2} {{}1}_n^\top \tilde{{{}W}} {{}T}_r ({{}K} - \phi {{}I}_n) {{}1}_n \, , \\
  {{}d}_{r, {{}\beta}\phi} & = \frac{1}{\phi} {{}X}^\top {{}Q} {{}T}_r ({{}K} - \phi {{}I}_n) {{}1}_n \, , &
  {{}d}_{r, \phi{{}\beta}} & = \frac{1}{\phi^3} {{}1}_n ^\top \tilde{{{}W}}^2 {{}T}_r {{}X} \, ,  \\
  {{}d}_{p + 1, {{}\beta}{{}\beta}} & = \frac{1}{\phi} {{}X}^\top {{}F} \tilde{{{}W}} {{}X} \,, &
  {{}d}_{p + 1, \phi\phi} & = {{}1}_n^\top ({{}K} - \phi {{}I}_n) {{}1}_n \, , \\
  {{}d}_{p + 1, {{}\beta}\phi} & = {{}X}^\top {{}F} ( {{}K} - \phi {{}I}_n) {{}1}_n\, , &
  {{}d}_{p + 1, \phi{{}\beta}} & = \frac{1}{\phi} {{}1}_n^\top \tilde{{{}W}} {{}X} \, . \end{aligned}$$
