# Mean and median bias reduction: A concise review and application to
adjacent-categories logit models

Ioannis Kosmidis
<EMAIL_ADDRESS>

Department of Statistics, University of Warwick
Coventry CV4 7AL, UK

###### Abstract

The estimation of categorical response models using bias-reducing adjusted
score equations has seen extensive theoretical research and applied use. The
resulting estimates have been found to have superior frequentist properties to
what maximum likelihood generally delivers and to be finite, even in cases
where the maximum likelihood estimates are infinite. We briefly review mean
and median bias reduction of maximum likelihood estimates via adjusted score
equations in an illustration-driven way, and discuss their particular
equivariance properties under parameter transformations. We then apply mean
and median bias reduction to adjacent-categories logit models for ordinal
responses. We show how ready bias reduction procedures for Poisson log-linear
models can be used for mean and median bias reduction in adjacent-categories
logit models with proportional odds and mean bias-reduced estimation in models
with non-proportional odds. As in binomial logistic regression, the reduced-
bias estimates are found to be finite even in cases where the maximum
likelihood estimates are infinite. We also use the approximation of the bias
of transformations of mean bias-reduced estimators to correct for the mean
bias of model-based ordinal superiority measures. All developments are
motivated and illustrated using real-data case studies and simulations.
Keywords: infinite estimates; bias reduction; adjusted score equations; data
separation

##  1 Overview

The first part of this chapter provides an example-driven, concise review of
the developments in a fast growing body of literature about mean and median
bias reduction (BR) in parametric estimation via adjusted score equations; see
Firth (1993) for mean BR (mBR) and Kenne Pagui et al. (2017) for median BR
(mdBR). Particular focus is placed on how these methods can be used as a
remedy for the numerical and inferential consequences of boundary maximum
likelihood (ML) estimates in categorical response models, which are
illustrated in Section 2. Sections 3 and 4 describe how the mean and median
bias of the ML estimator can be reduced in general parametric models through
the appropriate adjustment of the gradient of the log-likelihood. Section 5
discusses the validity of inference when the ML estimates are replaced by mBR
or mdBR estimates in standard first-order procedures. Section 6 takes a close
look at the equivariance properties of mBR and mdBR estimators under
transformation of the model parameters. We also present an approximation of
the bias of general transformations of mBR estimators, which can be used to
correct for the bias of transformations of the model parameters using only the
mBR estimates, the second derivatives of the transformation, and the expected
information matrix. The bias approximation is used to get mBR estimates of
odds ratios from mBR estimates of regression coefficients in logistic
regression models.

The second part of this chapter uses the results from the first to develop,
for the first time, mBR and mdBR procedures for adjacent-categories logit
(ACL) models for ordinal responses (see, for example, Agresti, 2010, Chapter 4
for an introduction). Section 7 reviews the proportional odds (PO) and non-
proportional odds (NPO) versions of the ACL models, and their key properties,
including their equivalence to baseline-category logit (BCL) models, and
discusses how that equivalence can be exploited for ML estimation. A real-data
case study is used to illustrate that boundary estimates can also cause
numerical and inferential issues for ACL models. Section 8 then details how
and when the equivariance properties of mBR and mdBR, and implementations of
the latter for BCL models, can be used for mBR and mdBR for the PO and NPO
versions of ACL models. Finally, Section 9 details how the mBR estimates can
be used for the explicit correction of the estimates of ordinal superiority
summaries.

##  2 Boundary estimates in categorical response models

It is well known that ML estimation of regression models with categorical
responses may result in estimates on the boundary of the parameter space. The
data patterns that result in boundary estimates in general multinomial
logistic regression models (also known as baseline category models models; see
Agresti, 2002, Section 7.1) have been studied extensively and are completely
characterized. For a range of binomial regression models, Silvapulle (1981)
proves that a certain degree of “overlap” on the data is a necessary and
sufficient condition for the ML estimates to have finite values. Albert and
Anderson (1984) enrich the arguments in Silvapulle (1981) generalizing the
results in the case of baseline-category logit (BCL) models for nominal
responses. In particular, Albert and Anderson (1984) categorize the possible
configurations for the sample points into complete separation, quasi-complete
separation, and overlap, and then show that separation is necessary and
sufficient for the ML estimate to have at least one infinite-valued component.
Geometric representations of (quasi-)complete separation for binomial logistic
regression — when the ACL and BCL models reduce to exactly the same form — are
given in Albert and Anderson (1984, Figure 1), and for multinomial responses
in Lesaffre and Albert (1989, Figure 1).

######  Example 1:

Separation in logistic regression A simple illustration of a completely
separated data set is shown in Figure 1. The data consists of $100$
realizations of two continuous covariates $x_{2}$ and $x_{3}$, and a response
$y$ that ends up being $0$ whenever $x_{2}+2x_{3}>0$. ML estimation of the
logistic regression model with
$\log\\{\pi/(1-\pi)\\}=\beta_{1}+\beta_{2}x_{2}+\beta_{3}x_{3}$, where $\pi$
is the probability of observing $y=1$ given $x$, results in the estimated
logistic discriminant line in Figure 1, with the log-likelihood attaining its
global maximum value of $0$, and the fitted value $0$ being assigned to all
observations with $y=0$, and $1$ to the rest. The detectseparation R package
(Kosmidis and Schumacher, 2021) that implements the methods in the unpublished
PhD thesis by Konis (2007) can be used to show that the ML estimates of
$\beta_{1}$, $\beta_{2}$ and $\beta_{3}$ are $-\infty$, $+\infty$ and
$+\infty$ respectively.

Table 1: ML, mBR and mdBR estimates for the logistic regression model in Example 1. The estimated standard errors (S.E.) are based on the expected information matrix at the estimates. The $z$-statistic is computed as estimate over estimated S.E., and the $p$-value is computed as $2\min(\Phi($z$),1-\Phi($z$))$, where $\Phi(\cdot)$ is the cumulative distribution function of the standard Normal distribution. Parameter | Estimate | Estimated S.E. |  $z$-statistic |  $p$-value
---|---|---|---|---
Maximum likelihood
$\beta_{1}$ | -22.397 | 13879.616 | -0.002 | 0.999
$\beta_{2}$ | 62.578 | 20968.761 | 0.003 | 0.998
$\beta_{3}$ | 132.228 | 44964.541 | 0.003 | 0.998
Mean bias reduction
$\beta_{1}$ | -2.001 | 1.552 | -1.289 | 0.197
$\beta_{2}$ | 5.266 | 1.997 | 2.637 | 0.008
$\beta_{3}$ | 11.166 | 3.984 | 2.803 | 0.005
Median bias reduction
$\beta_{1}$ | -2.583 | 1.984 | -1.302 | 0.193
$\beta_{2}$ | 6.325 | 2.628 | 2.406 | 0.016
$\beta_{3}$ | 13.321 | 5.293 | 2.517 | 0.012
Figure 1: The data described in Example 1. The dashed line is the line
$0=\hat{\beta}_{2}x_{2}+\hat{\beta}_{3}x_{3}$, where the fitted probabilities
are all $0.5$.

While there is no ambiguity in reporting infinite estimates, estimates on the
boundary of the parameter space can i) cause numerical instabilities to
fitting procedures, ii) lead to misleading output when estimation is based on
iterative procedures with a stopping criterion, and more importantly, iii)
cause havoc to asymptotic inferential procedures, and especially to the ones
that depend on estimates of the standard error of the estimators (for example,
Wald tests and related confidence intervals), oftentimes leading to wrong
inferences. For example, the ML estimates in Table 1 have been obtained using
the glm() function in R (R Core Team, 2021). Despite the fact that the ML
estimates for $\beta_{1}$, $\beta_{2}$ and $\beta_{3}$ are in reality
infinite, the stopping criteria of the fitting procedure that glm() implements
are met for finite values of the parameters, which are returned. The reported
estimated standard errors are also finite and substantially larger than the
estimates. This results in small, in absolute value, $z$-statistics, and hence
no evidence against the individual hypotheses $\beta_{2}=0$ and $\beta_{3}=0$;
one would expect at least some evidence against the hypotheses given that the
value of the response has been fully determined from the values of $x_{2}$ and
$x_{3}$.

One way to circumvent the numerical and inferential issues associated with
boundary ML estimates is to replace ML with an alternative estimation method
that i) has comparable or sometimes better asymptotic properties than the ML
estimator generally does, and ii) tends to result or results in estimates away
from the boundary of the parameter space. Popular examples of such alternative
estimation methods are the mean bias-reducing adjusted score functions
approach in Firth (1993), and the median bias-reducing adjusted score
functions approach in Kenne Pagui et al. (2017), which we briefly review in
Section 3, Section 4, and Section 6.

##  3 Mean bias reduction

Let $\ell(\theta)$ be the log-likelihood about a parameter vector $\theta$
with $\theta\in\Re^{v}$. Assuming that the model at hand is appropriate, then
under fairly general regularity conditions about the model, the ML estimator
$\hat{\theta}=\arg\max\ell(\theta)$ has mean bias ${\mathop{\rm
E}}_{\theta}(\hat{\theta}-\theta)=O(N^{-1})$, where $N$ is a measure of
information about $\theta$, usually — but not necessarily — the sample size.

If $S(\theta)=\nabla\ell(\theta)$, Firth (1993) shows that we can define an
alternative estimator $\theta^{*}$ with mean bias ${\mathop{\rm
E}}_{\theta}(\theta^{*}-\theta)=O(N^{-2})$, which is asymptotically smaller
than the bias of $\hat{\theta}$, as the solution of

$S(\theta)+A(\theta)=0_{v}\,,$ (1)

where

$A_{t}(\theta)=\frac{1}{2}{\mathop{\rm
trace}}\left[i(\theta)^{-1}\left\\{P_{t}(\theta)+Q_{t}(\theta)\right\\}\right]\quad(t=1,\ldots,v)\,.$

In the above expression, $P_{t}(\theta)={\mathop{\rm
E}}_{\theta}(S(\theta)S(\theta)^{\top}S_{t}(\theta))$ and
$Q_{t}(\theta)=-{\mathop{\rm E}}_{\theta}(j(\theta)S_{t}(\theta))$, and
$j(\theta)=-\nabla\nabla^{\top}\ell(\theta)$ and $i(\theta)={\mathop{\rm
E}}_{\theta}(S(\theta)S(\theta)^{\top})$ are the observed and expected
information matrix about $\theta$, respectively, with all expectations taken
with respect to the model.

Mean bias reduction has been found to result in estimates away from the
boundary of the parameter space in a range of categorical data models; see,
for example, Firth (1993) and Heinze and Schemper (2002) for binomial logistic
regression; Mehrabi and Matthews (1995) for the estimation of simple
complementary log–log-models; Kosmidis and Firth (2009, Section 6) for row-
column association models; Bull et al. (2002), Kosmidis and Firth (2011), and
Kosmidis et al. (2020, Section 6) for BCL models; and Kosmidis (2014) for
cumulative link models.

If $\theta$ is the canonical parameter of a full exponential family (see Pace
and Salvan, 1997, Chapter 5), like in binomial and multinomial logistic
regression, then $j(\theta)=i(\theta)$ and $j(\theta)$ does not depend on the
stochastic part of the model. Hence, $Q_{t}(\theta)=0_{v\times v}$, where
$0_{v\times v}$ is a $v\times v$ matrix of zeros, and some algebra (see Firth,
1993, Section 3) gives that the solution of the mean bias-reducing adjusted
score equations (1) is equivalent to the maximization of the penalized log-
likelihood

$\ell(\theta)+\frac{1}{2}\log{\mathop{\rm det}}\\{i(\theta)\\}\,,$ (2)

where the penalty is the logarithm of the Jeffreys prior. Recent work by
Kosmidis and Firth (2021) considers the impact of penalized likelihoods like
(2) in the estimation of many well-used binomial-response generalized linear
models, including logistic, probit, complementary log-log, and cauchit
regression. Among other results, Kosmidis and Firth (2021) prove that
maximizing the likelihood after penalizing it by arbitrary positive powers of
the Jeffreys prior always results in finite estimates, and derive the
shrinkage directions implied by the penalty.

##  4 Median bias reduction

The median bias-reducing adjusted score functions of Kenne Pagui et al. (2017)
is another method that has been found to result in finite estimates in
extensive simulation studies with logistic regression and BCL models (see
Kosmidis et al., 2020, Section 6)) and with cumulative link models (Gioia et
al., 2021).

The ML estimator generally has median bias
$P(\hat{\theta}_{t}\leq\theta_{t})=1/2+O(N^{-1/2})$. Kenne Pagui et al. (2017)
show that we can define an alternative estimator $\theta^{\dagger}$ with
$P(\theta_{t}^{\dagger}\leq\theta_{t})=1/2+O(N^{-3/2})$, which is
asymptotically closer to $1/2$ than the median bias of $\hat{\theta}$, as the
solution of

$S(\theta)+A(\theta)-i(\theta)F(\theta)=0_{v}\,.$ (3)

In the above expression,
$F_{t}(\theta)=[i(\theta)^{-1}]_{t}^{\top}\tilde{F}_{t}(\theta)$, with

$\tilde{F}_{tu}(\theta)={\mathop{\rm
trace}}\left[\tilde{i}_{u}(\theta)\left\\{\frac{1}{3}P_{t}(\theta)+\frac{1}{2}Q_{t}(\theta)\right\\}\right]\quad(t=1,\ldots,g)\,,$

and
$\tilde{i}_{u}(\theta)=[i(\theta)^{-1}]_{u}[i(\theta)^{-1}]_{u}^{\top}/[i(\theta)^{-1}]_{uu}$
$(u=1,\ldots,v)$, where $A_{u}$ and $A_{tu}$ denote the $u$th column and
$(t,u)$th element of a matrix $A$.

When $j(\theta)=i(\theta)$, expression (3) simplifies in a similar manner as
expression (1) does. In fact, for one-parameter models ($v=1$) that are
exponential families in canonical parameterization, it can be shown that mdBR
is formally equivalent to the maximization of $\ell(\theta)+\log{\mathop{\rm
det}}\\{i(\theta)\\}/6$ (see Kenne Pagui et al., 2017, Section 2.1). However,
mdBR has no penalized likelihood interpretation for $v>1$.

##  5 Inference with mean and median bias reduction

According to the results in Firth (1993) and Kenne Pagui et al. (2017), both
$\theta^{*}$ and $\theta^{\dagger}$ have the same asymptotic distribution as
the ML estimator generally does, and hence are asymptotically efficient.
Therefore, the distribution of those estimators for finite samples can be
approximated by a Normal with mean $\theta$ and variance-covariance matrix
$\\{i(\theta)\\}^{-1}$. The derivation of this result relies on the fact that
both the adjustments $A(\theta)$ and $A(\theta)-i(\theta)F(\theta)$ to the
score functions for mBR and mdBR in (1) and (3), respectively, are of order
$O(1)$ as $N\to\infty$. Hence, the score function $S(\theta)$, which is
$O_{p}(\sqrt{N})$, dominates the adjustments as information increases. The
implication is that standard errors for the components of $\theta^{*}$ and
$\theta^{\dagger}$ can be computed exactly as for the ML estimator, using the
square roots of the diagonal elements of $\\{i(\theta)\\}^{-1}$ of
$\\{j(\theta)\\}^{-1}$ at the estimates. Furthermore, first-order inferences,
like standard Wald tests and Wald-type confidence intervals and regions are
constructed in a plugin fashion, by replacing the ML estimates with the mBR or
mdBR estimates in the usual procedures in standard software.

######  Example 2:

Separation in logistic regression (continued) Continuing from Example 1, Table
1 provides the estimates of $\beta_{1}$, $\beta_{2}$ and $\beta_{3}$ from mBR
and mdBR. The estimates have been computed using the default arguments of the
brglm_fit() method of the brglm2 R package (Kosmidis, 2021). brglm_fit()
implements a variant of the quasi-Fisher scoring procedure

$\theta^{(k+1)}=\theta^{(k)}+\\{i(\theta^{(k)})\\}^{-1}U(\theta^{(k)})\,,$ (4)

where $U(\theta):=S(\theta)+A(\theta)$ if the intention is to compute the mean
BR estimates, and $U(\theta):=S(\theta)+A(\theta)-i(\theta)F(\theta)$ if the
intention is to compute the mdBR estimates; see Kosmidis et al. (2020) for
details on the quasi-Fisher iterations and the form of the adjusted scores for
mBR and mdBR in generalized linear models. Convergence has been rapid and
brglm_fit() reported no issues for either mBR or mdBR. Furthermore, the
estimates and estimated standard errors appear to be finite. Note that the
estimates and estimated standard errors from mBR are typically closer in
absolute value to zero than those from mdBR. Importantly, the $z$-statistics
for $\beta_{2}$ and $\beta_{3}$ are all away from zero, and, in contrast to
ML, both mBR and mdBR suggest at least some evidence against the individual
hypothesis $\beta_{2}=0$ and $\beta_{3}=0$, which agrees with the fact that
the value of the response has been fully determined from the values of $x_{2}$
and $x_{3}$.

##  6 Bias reduction and parameter transformation

###  6.1 Maximum likelihood estimation and general parameter transformations

The ML estimator is equivariant in the sense that the ML estimator of
$g(\theta)$ is exactly $g(\hat{\theta})$ for any one-to-one transformation
$g(\cdot)$. Hence, there is no need to maximize the log-likelihood about
$g(\theta)$ if the ML estimator of $\theta$ has already been computed. In
contrast, the mBR and mdBR estimators are equivariant only for specific
transformations $g(\cdot)$.

###  6.2 Mean bias reduction and linear parameter transformations

The mBR estimator is equivariant under linear transformations for the
parameters, in the sense that the mBR estimator of $C\theta$ for a known
matrix $C$ is exactly $C\theta^{*}$. The same is not true for the mdBR
estimator.

For example, using Table 1, the mBR estimate of $\beta_{2}-\beta_{3}$ in
Example 2 is simply $5.266-11.166=-5.9$. The mdBR estimate, however, is not
$6.325-13.321=6.996$, but rather $-7.227$, which is obtained by
reparameterizing the model in terms of $\beta_{2}-\beta_{3}$ and computing the
mdBR estimate by solving (3) in the new parameterization.

###  6.3 Median bias reduction and component-wise parameter transformations

On the other hand, the mdBR estimator of
$(g_{1}(\theta_{1}),\ldots,g_{v}(\theta_{v}))^{\top}$ is
$(g_{1}(\theta_{1}^{\dagger}),\ldots,g_{v}(\theta_{v}^{\dagger}))^{\top}$ for
any set of one-to-one functions $g_{1}(\cdot),\ldots,g_{v}(\cdot)$. In other
words, the mdBR estimator is equivariant under component-wise transformations.
The same is not true for the mBR estimator. For example, the mdBR estimate of
the odds-ratio $\exp(\beta_{2})$ in Example 2 is exactly $\exp(6.325)$, but
$\exp(5.266)$ is not an mBR estimate of $\exp(\beta_{2})$.

###  6.4 Mean bias reduction and general parameter transformations

Di Caterina and Kosmidis (2019) show that there is a simple way to derive the
mean bias of $h(\theta^{*})$ for any three-times differentiable function
$h:C\to D$, with $C\subset\Re^{p}$ and $D\subset\Re$, where $\theta^{*}$ is an
mBR estimator of $\theta$ with $O(N^{-2})$ bias. In particular, Di Caterina
and Kosmidis (2019) show that the estimator $h(\theta^{*})$ of
$\zeta=h(\theta)$ has mean bias

$E(h(\theta^{*})-h(\theta))=\frac{1}{2}{\mathop{\rm
trace}}\left\\{i(\theta)^{-1}\nabla\nabla^{\top}h(\theta)\right\\}+O(N^{-2})\,,$
(5)

where $\nabla\nabla^{\top}h(\theta)$ is the hessian of $h(\cdot)$ at $\theta$.
Note that for linear transformations, $\nabla\nabla^{\top}h(\theta)=0_{v\times
v}$, and hence $E(h(\theta^{*})-h(\theta))=O(N^{-2})$, which confirms the
discussion in Section 6.2 that the mBR estimator is exactly equivariant for
linear transformations of the parameters. The first term in the right-hand
side of (5) can be evaluated at $\theta^{*}$ and be used to derive mean BR
estimators of $h(\theta)$, based only on $\hat{\theta}^{*}$,
$i(\hat{\theta}^{*})$, and $\nabla\nabla^{\top}h(\theta^{*})$. An obvious mean
BR estimator resulting from (5) is $h(\theta^{*})-{\mathop{\rm
trace}}\left\\{i(\theta^{*})^{-1}\nabla\nabla^{\top}h(\theta^{*})\right\\}/2$.

For example, consider the special case of estimation of the odds-ratio
$\exp(\beta_{j})$ in Example 1, which was estimated using the equivariance
properties of mdBR in Section 6.3. Expression (5) gives that the odds-ratio at
the mBR estimator has

$E(\exp(\beta_{j}^{*}))=\exp(\beta_{j})\left[1+\frac{1}{2}v_{jj}(\theta)\right]+O(N^{-2})\,,$
(6)

where $v_{jj}(\theta)=[i(\theta)^{-1}]_{jj}$. Hence, two mean BR estimators of
$\zeta_{j}=\exp(\beta_{j})$ with $O(N^{-2})$ bias are

$\zeta_{j}^{*}=\exp(\beta_{j}^{*})\left[1-\frac{1}{2}v_{jj}(\theta^{*})\right]\quad\text{and}\quad\zeta_{j}^{**}=\frac{\exp(\beta_{j}^{*})}{1+v_{jj}(\theta^{*})/2}\,,$

arising from subtracting an estimate of the bias at $\theta:=\theta^{*}$ from
$\exp(\beta_{j}^{*})$, and dividing $\exp(\beta_{j}^{*})$ by the correction
factor $1+v_{jj}(\theta^{*})/2$ from the right-hand side of (6), respectively.
The estimator $\zeta_{j}^{**}$ for the odds-ratio $\zeta_{j}$ has the
advantage of being always positive, while $\zeta_{j}^{*}$ takes negative
values if $v_{jj}(\theta^{*})>2$. For example, to the accuracy reported in
Table 1, $\zeta_{2}^{*}=\exp(5.266)(1-1.997^{2}/2)=-192.48$, which is clearly
nonsensical as an odds-ratio estimate. In contrast,
$\zeta_{2}^{*}=\exp(5.266)/(1+1.997^{2}/2)=64.66$. The approximation
$\exp\\{v_{jj}(\theta)/2\\}\approx 1+v_{jj}(\theta)/2$ for small
$v_{jj}(\theta)$ can be used to show that the mean BR estimator
$\zeta_{j}^{**}$ closely relates to the mean BR estimator
$\zeta_{j}^{***}=\exp\\{\beta_{j}^{*}-v_{jj}(\theta^{*})/2\\}$ derived in
Lyles et al. (2012).

The discussion in Section 5 implies that estimated standard errors for mBR
estimators of transformed parameters constructed on the basis of (6) can be
computed using the delta method, as for the ML estimator.

##  7 Adjacent-categories logit models

###  7.1 Proportional and non-proportional odds models

We now turn our attention in applying mBR and mdBR from Section 3 and Section
4 to ACL models.

Adjacent-categories logit models (see, for example, Agresti, 2010, Chapter 4
for an introduction) are a prominent family of regression models for ordinal
responses, where the local odds ratios of consecutive categories of an ordinal
response variable are linked with linear combinations of parameters and
explanatory variables. Suppose that we observe realizations of $n$ independent
random vectors of frequencies $Y_{1},\ldots,Y_{n}$, where
$Y_{i}=(Y_{i1},\ldots,Y_{ik})^{\top}$ has a $k$-category multinomial
distribution with ordered categories $1<2<\ldots<k$, total
$m_{i}=\sum_{j=1}^{k}Y_{ij}$ and probability vector
$(\pi_{1}(x_{i}),\ldots,\pi_{k}(x_{i}))^{\top}$ with
$\sum_{j=1}^{k}\pi_{j}(x_{i})=1$, where $x_{i}=(x_{i1},\ldots,x_{ip})^{\top}$
is a $p$-vector of covariate values. An ACL model has

$\log\frac{\pi_{j}(x)}{\pi_{j+1}(x)}=\eta_{j}(x)\quad(j=1,\ldots,k-1)\,,$ (7)

where $\eta_{j}(x)$ is typically a linear combination of unknown model
parameters and a covariate vector $x$.

The specification of $\eta_{j}(x)$ results in ACL models with particular
properties. The PO version of the ACL model has

$\eta_{j}(x)=\alpha_{j}+\beta^{\top}x\,,$ (8)

and $p+k-1$ scalar model parameters
$\theta=(\alpha_{1},\ldots,\alpha_{k-1},\beta_{1},\ldots,\beta_{p})^{\top}$.
Straightforward algebra starting from (7) gives that

$\frac{\pi_{j}(x_{2})}{\pi_{j+1}(x_{2})}=\exp\\{\beta^{\top}(x_{2}-x_{1})\\}\frac{\pi_{j}(x_{1})}{\pi_{j+1}(x_{1})}\quad\text{for
any }x_{1},x_{2}\in\Re^{p}\text{ and }j\in\\{1,\ldots,k-1\\}\,.$ (9)

As a result, adjacent-categories odds are indeed proportional with a constant
of proportionality that does not depend on the category. The NPO version of
the ACL model has

$\eta_{j}(x)=\alpha_{j}+\beta_{j}^{\top}x\,,$ (10)

with $(k-1)(p+1)$ scalar model parameters
$\theta=(\alpha_{1},\ldots,\alpha_{k-1},\beta_{1}^{\top},\ldots,\beta_{k-1}^{\top})^{\top}$,
where $\beta_{j}=(\beta_{j1},\ldots,\beta_{jp})^{\top}$. Figure 2 shows the
adjacent-categories log-odds for two distinct categories under the PO and the
NPO versions of the ACL model, for $x\in\Re$. Note that under the PO version
of the model the log-odds for distinct categories are parallel lines, which in
turn implies (9) for any pair of categories. On the other hand, under the NPO
version of the model the log-odds are not parallel lines, so (9) is not
generally satisfied.

Figure 2: The adjacent-categories log-odds for categories $j$ and $k$, $j\neq
k$, under the proportional odds (left) and the non-proportional odds (right)
versions of the model, for $x\in\Re$. The probability for category $j$ at
covariate value $x$ is denoted
$\pi_{j}(x)$.$x_{1}$$x_{2}$$\log\frac{\pi_{j}(x)}{\pi_{j+1}(x)}$$\beta_{j}(x_{2}-x_{1})$$\log\frac{\pi_{k}(x)}{\pi_{k+1}(x)}$$\beta_{k}(x_{2}-x_{1})$$x_{1}$$x_{2}$$\log\frac{\pi_{j}(x)}{\pi_{j+1}(x)}$$\beta(x_{2}-x_{1})$$\log\frac{\pi_{k}(x)}{\pi_{k+1}(x)}$$\beta(x_{2}-x_{1})$

The most general version of the ACL model is the partial proportional odds
model with

$\eta_{j}(x)=\alpha_{j}+\xi_{j}^{\top}x^{(np)}+\rho^{\top}x^{(p)}\,,$

where $x^{(np)}$ and $x^{(p)}$ are sub-vectors of $x$ with distinct components
characterizing the PO and NPO effects, respectively. All subsequent
derivations, results, and discussions can be written in terms of the more
general partial proportional odds version, and then PO and NPO can be
presented as special cases. Nevertheless, we focus on the PO and NPO versions
separately, to keep the notation concise, and because some of the following
results are specific to PO and not to NPO.

Expressions (8) and (10) immediately imply that the ACL model provides valid
category probabilities across the parameter space and regardless of whether
the local odds $\pi_{j}(x)/\pi_{j+1}(x)$ are modelled as proportional or non-
proportional. This is in contrast to other popular ordinal-response regression
models, like cumulative-logit models (McCullagh, 1980), whose NPO versions
(Peterson and Harrell, 1990) may provide invalid category probabilities in
subsets of the parameter space and covariate space, and, hence, result in
hard-to-circumvent issues with estimation, inference, and prediction.

###  7.2 Equivalence with baseline-category logit models

Writing
$\log\\{\pi_{j}(x)/\pi_{k}(x)\\}=\sum_{l=j}^{k-1}\log\\{\pi_{l}(x)/\pi_{l+1}(x)\\}$,
it is simple to show that both the PO and NPO versions of the ACL model for
ordinal responses can be written as BCL models for nominal responses (see
Agresti, 2010, Section 4.1) where the $k$ category is used as reference.

In particular, the NPO version of the ACL model in (10) is equivalent to a BCL
model with

$\log\frac{\pi_{j}(x)}{\pi_{k}(x)}=\gamma_{j}+\delta_{j}^{\top}x\quad(j=1,\ldots,k-1)\,,$
(11)

where $\gamma_{j}=\sum_{l=j}^{k-1}\alpha_{l}$ and
$\delta_{j}=\sum_{l=j}^{k-1}\beta_{l}$. The PO version of the ACL model in (8)
is equivalent to a BCL model with

$\log\frac{\pi_{j}(x)}{\pi_{k}(x)}=\gamma_{j}+(k-j)\zeta^{\top}x\quad(j=1,\ldots,k-1)\,,$
(12)

where $\gamma_{j}=\sum_{l=j}^{k-1}\alpha_{l}$ and $\beta=\zeta$.

###  7.3 Maximum likelihood estimation

A consequence of the equivalence between the BCL and ACL models is that we can
estimate the latter using the ML estimates for the former. The equivariance of
the maximum ML estimator under one-to-one transformations of the model
parameters guarantees that after computing the ML estimates for the parameters
of BCL model (11), the model parameters of the NPO version of the ACL can be
estimated as $\hat{\alpha}_{j}=\hat{\gamma}_{j}-\hat{\gamma}_{j+1}$ and
$\hat{\beta}_{j}=\hat{\delta}_{j}-\hat{\delta}_{j+1}$ $(j=1,\ldots,k-1)$ with
$\hat{\gamma}_{k}=0$ and $\hat{\beta}_{k}=0_{p}$, where $0_{p}$ is a
$p$-vector of zeros. Correspondingly, once the ML estimates for the parameters
of BCL model (12) have been obtained, the model parameters of the PO version
of the ACL model can be estimated as
$\hat{\alpha}_{j}=\hat{\gamma}_{j}-\hat{\gamma}_{j+1}$ and
$\hat{\beta}=\hat{\zeta}$ $(j=1,\ldots,k-1)$.

So, ML estimation of ACL models can be performed using ready ML
implementations for fitting the BCL models (11) and (12), like the multinom()
function of the nnet R package (Venables and Ripley, 2002) that exploits the
equivalence of BCL models with neural networks, and the brmultinom() function
of the brglm2 R package (Kosmidis, 2021) that exploits the equivalence of BCL
models with Poisson log-linear models.

###  7.4 Exponential families

The BCL model is a full exponential family distribution with natural
parameters $\gamma_{j}$ and $\delta_{j}$ for the NPO version (11) of the ACL
model, and $\gamma_{j}$ and $\zeta$ for the PO version (12) of the ACL model
$(j=1,\ldots,k-1)$. Hence, another consequence of the equivalence of ACL
models to BCL models is that both the PO and NPO versions of the ACL model are
full exponential families. Specifically, the sufficient statistics in the NPO
parameterization are $\sum_{l=1}^{j}\sum_{i=1}^{n}y_{il}$ for $\alpha_{j}$,
$\sum_{l=1}^{j}\sum_{i=1}^{n}y_{il}x_{i}$ for $\beta_{j}$, and
$\sum_{l=1}^{j}\sum_{i=1}^{n}y_{il}$ for $\alpha_{j}$ and
$\sum_{j=1}^{k-1}\sum_{i=1}^{n}(k-j)y_{ij}x_{i}$ for $\beta$ in the PO
parameterization $(j=1,\ldots,k-1)$ (see, also, Agresti, 2010, Section 4.1).

###  7.5 Infinite maximum likelihood estimates

As is the case for their equivalent BCL models, depending on the data
configuration, the ML estimates of ACL models can have infinite components,
resulting in issues for both iterative estimation procedures and for first-
order inference about the parameters. In fact, infinite ML estimates for the
PO and NPO versions of the ACL model result if, and only if, separation occurs
for the equivalent BCL models. Example 3 below uses a real data set to
illustrate the consequences that separation can have in the estimation of, and
inference from, ACL models.

######  Example 3:

Infinite ML estimates in ACL models The data set in Table 2 comes from Randall
(1989) and concerns an experiment for investigating factors that affect the
bitterness of white wine. There are two factors in the experiment, namely
temperature at the time of crushing the grapes (with two levels, “cold” and
“warm”) and contact of the juice with the skin (with two levels “Yes” and
“No”). For each combination of factors two bottles were rated on their
bitterness by a panel of 9 judges. The responses of the judges on the
bitterness of the wine were taken on a continuous scale in the interval from 0
(“None”) to 100 (“Intense”) and then they were grouped correspondingly into
$5$ ordered categories, labelled as “1”, “2”, “3”, “4”, and “5”.

Table 2: The wine tasting data (Randall, 1989). Temperature | Contact | Bitterness rating
---|---|---
1 | 2 | 3 | 4 | 5
Cold | No | 4 | 9 | 5 | 0 | 0
Cold | Yes | 1 | 7 | 8 | 2 | 0
Warm | No | 0 | 5 | 8 | 3 | 2
Warm | Yes | 0 | 1 | 5 | 7 | 5
Figure 3: The empirical adjacent logits
$\log\\{(y_{ij}+1/2)/(y_{ij+1}+1/2)\\}$ $(j=1,\ldots,4)$ for the bitterness
rating for all combinations of levels for temperature and contact.

Figure 3 shows the empirical adjacent logits
$\log\\{(y_{ij}+1/2)/(y_{ij+1}+1/2)\\}$ $(j=1,\ldots,4)$ for the bitterness
rating for all combinations of temperature and contact. Note that $1/2$ has
been added to all frequencies as a means of getting estimates of the adjacent-
categories logits with second-order mean bias (see, for example, Haldane,
1955), avoiding infinite estimates in the process.

There seems to be evidence that the adjacent logits for the combinations of
temperature and contact are parallel (see, also, Figure 2), or in other words,
the adjacent odds ratios across temperature and/or contact levels do not
depend on the rating. The latter hypothesis can be formally tested by
estimating the NPO version of the ACL model

$\log\frac{\pi_{j}(t,c)}{\pi_{j+1}(t,c)}=\alpha_{j}+\beta_{1j}t+\beta_{2j}c\quad(j=1,\ldots,4)\,,$
(13)

where $t$ is $1$ if temperature is warm and $0$ otherwise, $c$ is $1$ if
contact is yes and $0$ otherwise, and $\pi_{j}(t,c)$ is the probability of a
bitterness rating $j$ at $t$ and $c$. The hypotheses of parallel adjacent
logits can then be written in terms of the model parameters as
$\beta_{11}=\ldots=\beta_{14}=\beta_{1}$ and
$\beta_{21}=\ldots=\beta_{24}=\beta_{2}$, and tested using the value of the
Wald statistic

$W=\hat{\theta}^{\top}C^{\top}\left\\{Ci(\hat{\theta})^{-1}C^{\top}\right\\}^{-1}C\hat{\theta}\,,$
(14)

where $\hat{\theta}$ is the ML estimate of
$\theta=(\alpha_{1},\ldots,\alpha_{4},\beta_{11},\ldots,\beta_{14},\beta_{21},\ldots,\beta_{24})^{\top}$
for model (13), and $i(\theta)$ is the expected information matrix at
$\theta$. The contrast matrix $C$ we use in (14) has the form

$C=\left[\begin{array}[]{ccc}0_{3\times 4}&C_{1}&0_{3\times 4}\\\ 0_{3\times
4}&0_{3\times 4}&C_{1}\end{array}\right]\quad\text{with}\quad
C_{1}=\left[\begin{array}[]{cccc}1&0&0&-1\\\ 0&1&0&-1\\\
0&0&1&-1\end{array}\right]\,,$

where $0_{a\times b}$ is an $a\times b$ matrix of zeros. General results about
the limiting distribution of the ML estimator under mild regularity conditions
(see, for example, McCullagh, 2018, Section 7.1 and Section 7.2 and Cox and
Hinkley, 1974, Section 9.1) can be used to show that the Wald statistic has
asymptotically a $\chi^{2}_{6}$ distribution.

Table 3: Top: ML estimates and estimated standard errors (in parenthesis) from fitting the ACL model in (13) on the data in Table 2. The estimates are obtained using the vglm() function of the VGAM R package (Yee, 2021) version 1.1-5 with default converge criteria (epsilon = $10^{-7}$ in vglm.control()). Bottom: ML estimates and estimated standard errors using stricter convergence criteria (epsilon = $10^{-9}$ in vglm.control()). The estimated standard errors are computed as the square roots of the diagonal of the inverse of the expected information matrix at the ML estimates. The column $\ell(\hat{\theta})$ gives the maximized log-likelihood for each fit. epsilon | $\ell(\hat{\theta})$ | rating ($j$) | $\hat{\alpha}_{j}$ | $\hat{\beta}_{1j}$ | $\hat{\beta}_{2j}$
---|---|---|---|---|---
$10^{-7}$ | -15.29 | 1 | -0.83 | (0.59) | -20.26 | (10047.96) | -1.10 | (1.21)
2 | 0.67 | (0.52) | -1.21 | (0.66) | -0.87 | (0.64)
3 | 3.08 | (1.05) | -1.98 | (0.92) | -1.54 | (0.83)
4 | 20.22 | (10732.18) | -19.89 | (10732.18) | 0.04 | (1.08)
$10^{-9}$ | -15.29 | 1 | -0.83 | (0.59) | -25.26 | (122409.18) | -1.10 | (1.21)
2 | 0.67 | (0.52) | -1.21 | (0.66) | -0.87 | (0.64)
3 | 3.08 | (1.05) | -1.98 | (0.92) | -1.54 | (0.83)
4 | 25.22 | (130748.2) | -24.89 | (130748.24) | 0.04 | (1.08)

Table 3 shows the ML estimates of the ACL model in (13), as computed using the
vglm() function of the VGAM R package (Yee, 2021). No warnings or errors were
returned when fitting the model. As has been the case in the logistic
regression model of Example 1, the estimates and estimated standard errors for
$\alpha_{4}$, $\beta_{11}$ and $\beta_{14}$ are atypically large in absolute
value. It is also clear that these estimates and estimated standard errors
increase in absolute value as the convergence criteria get stricter, while the
maximized log-likelihood value remains the same to the displayed accuracy.

These issues are not due to the implementation of the vglm() function; instead
they are consequences of quasi-complete separation for this particular
combination of data and model (13). The ML estimates $\hat{\alpha}_{4}$,
$\hat{\beta}_{11}$ and $\hat{\beta}_{14}$ in Table 3 are formally $\infty$,
$-\infty$ and $-\infty$, the corresponding estimated standard errors are all
$\infty$, and the likelihood surface has an asymptote at $-15.29$ as
$\alpha_{4}$, $\beta_{11}$ and $\beta_{14}$ diverge to $\infty$, $-\infty$ and
$-\infty$, respectively, along a ray in the parameter space.

Note here that the estimated standard errors appear to diverge faster than the
ML estimates do as the convergence criteria get stricter. As a result, the
typically reported $Z$-statistics for individual hypothesis tests about the
parameters will tend to be spuriously small in absolute value regardless of
the strength of the evidence against the hypotheses. Hence, the naive use of
the computer output for inference about the parameters of ACL models is likely
to lead to invalid conclusions when data separation occurs. More importantly,
having estimates on the boundary of the parameter space violates the
assumptions required for the asymptotic $\chi^{2}$ distribution of (14).
Consequently, it is hard to justify the performance and validity of the Wald
statistic in that case.

##  8 Mean and median bias reduction for ACL models

A consequence of the ACL models being full exponential family distributions
(see Section 7.4) is that mean BR can be implemented by maximizing the
penalized likelihood in (2). Nevertheless, as for ML, mean BR estimates for
ACL models can be conveniently computed through a ready implementation for
mean BR in BCL models coupled with the equivariance of the mean BR estimator
under linear transformations (see Section 6.2).

Kosmidis and Firth (2011) prove that the equivalence of BCL models and Poisson
log-linear models (see, also, Palmgren, 1981 and Baker, 1994 for authoritative
descriptions of that equivalence) extends to the mBR estimates, and describe a
simple algorithm for mBR estimation of BCL models, each iteration of which
consists of the following steps:

  * P1

Rescale the Poisson means to match the observed multinomial totals.

  * P2

Add half a leverage based on the rescaled means to the observed multinomial
frequencies.

  * P3

Estimate, using ML, the equivalent Poisson log-linear model to the adjusted
frequencies.

Iteration stops when the differences between successive estimates or,
alternatively, the mean BR adjusted scores in (1) are smaller than a pre-
determined, small positive constant. An alternative criterion can be based on
the change of the mean BR penalized likelihood (2) between successive
iterations. mBR estimates for ACL models can then be computed as follows

  * S1

Compute mBR estimates of the parameters $\gamma_{j}$ and $\delta_{j}$ of the
BCL model in (11) for the NPO version (or $\gamma_{j}$ and $\zeta$ of the BCL
model in (12) for the PO version) $(j=1,\ldots,q)$ by iterating steps P1, P2,
and P3.

  * S2

Calculate the mBR estimates for the NPO version of the ACL model as
$\alpha_{j}^{*}=\gamma_{j}^{*}-\gamma_{j+1}^{*}$ and
$\beta_{j}^{*}=\delta_{j}^{*}-\delta_{j+1}^{*}$ (or
$\alpha_{j}^{*}=\gamma_{j}^{*}-\gamma_{j+1}^{*}$ and $\beta^{*}=\zeta^{*}$ for
the PO version) $(j=1,\ldots,q)$, with $\gamma_{k}^{*}=0$ and
$\beta_{k}^{*}=0_{p}$.

Implementation of mdBR for ACL models is not as direct as that of mBR. A
maximum penalized likelihood interpretation of mdBR does not exist for general
ACL models, like it does for mBR. Also, since contrasts of parameters are not
component-wise transformations, algorithms for mdBR for BCL models (see
Kosmidis et al., 2020, Section 6 for extensions of the results in Kosmidis and
Firth, 2011) can only be used to get mdBR estimates $\beta^{\dagger}$ of
$\beta$ in the PO version of the ACL model. In other words, the estimates
$\gamma_{j}^{\dagger}-\gamma_{j+1}^{\dagger}$ and
$\beta_{j}^{\dagger}=\delta_{j}^{\dagger}-\delta_{j+1}^{\dagger}$
$(j=1,\ldots,q)$ are not mdBR estimates, unless $k=2$. Hence, for general ACL
models, computing the mdBR estimates $\theta^{\dagger}$ must rely on
implementing and solving the mdBR adjusted score equations (3). That can
certainly be done (using, for example, the quasi-Fisher scoring iteration
(4)), with the only effort being in deriving $P_{t}(\theta)$ using the
expressions for mBR in BCL models in Kosmidis (2007, Appendix B.5).

######  Example 4:

Infinite ML estimates in ACL models (continued) Table 4 gives the mBR
estimates from fitting the ACL model in (13) on the data in Table 2. The mBR
estimates are computed using the bracl() function of the brglm2 R package,
which implements mBR through the corresponding Poisson log-linear model, as
detailed earlier. No convergence issues have been reported; the absolute
values of the components of the adjusted score functions in (1) at the mBR
estimates are all less than $10^{-6}$, and all estimates and estimated
standard errors remain unchanged to the reported accuracy as the convergence
criteria get stricter.

The Wald statistic (14) when $\hat{\theta}$ is replaced by $\theta^{*}$ has
value $1.067$, which is small compared to the value of the $95\%$ quantile of
a $\chi^{2}_{6}$ distribution ($12.592$), providing no evidence against the
simpler PO model with $\beta_{11}=\ldots=\beta_{14}=\beta_{1}$ and
$\beta_{21}=\ldots=\beta_{24}=\beta_{2}$.

Table 4: Mean BR estimates and estimated standard errors (in parenthesis) from fitting the ACL model in (13) on the data in Table 2. The estimates are obtained using the bracl() function of the brglm2 R package (Kosmidis, 2021) version 0.7.2 with default convergence criteria. The estimated standard errors are computed as the square roots of the diagonal of $i(\theta^{*})^{-1}$. rating ($j$) | $\alpha_{j}^{*}$ | $\beta_{1j}^{*}$ | $\beta_{2j}^{*}$
---|---|---|---
1 | -0.76 | (0.59) | -1.65 | (1.60) | -0.82 | (1.08)
2 | 0.62 | (0.52) | -1.12 | (0.66) | -0.80 | (0.64)
3 | 2.73 | (0.99) | -1.75 | (0.87) | -1.38 | (0.81)
4 | 1.53 | (1.83) | -1.26 | (1.68) | 0.07 | (1.03)

Comparing the mBR estimates in Table 4 to the ML ones in Table 3, we notice
that the mBR estimates are shrunken towards zero relative to ML ones. As a
result, the fitted multinomial probabilities at the mBR estimates are closer
to $(1/5,1/5,1/5,1/5,1/5)^{\top}$ than ones at the ML estimates. In other
words, mBR shrinks the model towards equi-probability across observations.
This is a generalization of the shrinkage effect of mBR we observed in Example
2 and that Kosmidis and Firth (2021) study theoretically in the special case
of logistic regression ($k=2$).

It is interesting to note that the shrinkage direction of mBR in cumulative
logit models for global cumulative odds (Kosmidis, 2014) is rather different;
the fitted multinomial probabilities at the mBR estimates for the PO version
of the cumulative logit model would be closer to $(1/2,0,0,0,1/2)^{\top}$ than
the ones at the ML estimates. In other words, mBR shrinks the cumulative logit
model towards a logistic regression model for the end categories.

The ML, mdBR, and mBR estimates for $\beta_{1}$ for the PO version of the ACL
model are $-1.69$, $-1.61$, and $-1.56$, respectively, with corresponding
estimated standard errors $0.41$, $0.39$, and $0.38$. The respective estimates
for $\beta_{2}$ are $-0.96$, $-0.92$, and $-0.90$, respectively, with
corresponding estimated standard errors $0.32$, $0.31$, and $0.31$. The
shrinkage towards equi-probability that mBR delivers is also apparent in the
estimates for the PO version of the ACL model. As is the case in logistic
regression, mdBR also tends to shrink estimates towards zero, but that
shrinkage effect is less strong than from mBR.

##  9 Mean bias reduction of ordinal superiority summaries

The mean BR estimates for ACL models can be used to get improved estimates of
other model summaries by using the bias of transformations of the mean RB
estimator in expression (5).

A prominent example of such a summary are the ordinal, model-based superiority
measures for comparing distributions of two groups, adjusted for covariates
that are introduced in Agresti and Kateri (2017). In ordinal-response models
with a latent variable interpretation, such as cumulative-link models
(McCullagh, 1980), ordinal superiority measures can be defined directly on the
latent scale, which results in exact (for probit, log-log, and complementary
log-log link) or approximate expressions (for logit link) that are functions
of only the coefficient of the indicator variable characterizing the two
groups being compared. This fact has been exploited in Gioia et al. (2021),
who used the equivariance properties of the mdBR estimator (see Section 6.3)
to directly transform the mdBR estimates of the group indicator parameter to
deliver mdBR estimates of ordinal superiority measures.

In more general models for ordinal responses that may also lack a latent
variable interpretation (like ACL models), ordinal superiority measures are
instead defined in terms of category probabilities that necessarily depend on
all model parameters. Suppose that the covariate vector is
$(w^{\top},z)^{\top}$, where $z$ is a group indicator variable taking value
$0$ for group $1$ and value $0$ for group $2$, and denote by $\pi_{j}(w,1)$
and $\pi_{j}(w,0)$ $(j=1,\ldots,k)$ the model-based probabilities of category
$j$ at covariate values $w$, for group 1 and group 2, respectively. The
dependence of the probabilities on the model parameters has been suppressed
here for notational convenience.

Agresti and Kateri (2017) propose comparing the distribution of the ordinal
response at group 1 to that at group 2, at covariate values $w$, through the
ordinal superiority measure

$\Delta(w;\theta)=\sum_{r>s}\pi_{r}(w,1)\pi_{s}(w,0)-\sum_{s>r}\pi_{r}(w,1)\pi_{s}(w,0)\,.$
(15)

If the two distributions are identical then $\Delta(w;\theta)=0$. Positive
values of $\Delta(w;\theta)$ indicate that for covariates $w$, it is more
likely to observe higher response categories in group 1 than in group 2, and
vice versa for negative values. A related ordinal superiority measure is

$\gamma(w;\theta)=2\Delta(w;\theta)-1\,,$ (16)

which takes values between $0$ and $1$, and is interpreted as the probability
that the response category in group 1 is higher than the response category in
group 2, while adjusting for covariates $w$ (see Klotz, 1966, for details). In
practice, the covariate setting $w$ can be taken to be a representative value
from a sample of covariate values $w_{1},\ldots,w_{n}$, e.g.
$\bar{w}=\sum_{i=1}^{n}w_{i}/n$. Alternatively, if the sample of covariate
values is representative of the population of interest then summary ordinal
superiority measures can be defined as

$\bar{\Delta}(\theta)=\frac{1}{n}\sum\Delta(w_{i};\theta)\quad\text{and}\quad\bar{\gamma}(\theta)=\frac{1}{n}\sum\gamma(w_{i};\theta)\,.$
(17)

Agresti and Kateri (2017) propose estimating the ordinal superiority measures
by replacing $\theta$ in expressions (15), (16), and (17) by the ML estimator
$\hat{\theta}$, and use the delta method to construct inferences about those
measures. Note here that because of the specific equivariance properties of
the mBR and mdBR estimator (see Section 6), replacing $\theta$ by the mBR
estimator $\theta^{*}$ or a mdBR estimator $\theta^{\dagger}$ does not, in
general, result in mBR or mdBR estimators of the measures. In fact, despite it
being the case that the resulting estimators will be consistent under the same
conditions that their ML counterparts are, they may end up having much worse
finite-sample mean and/or median bias properties than the ML version does.

mdBR estimators of (15), (16), and (17) are not easy to construct. In
contrast, an easy-to-compute mBR estimator of $\Delta(w;\theta)$ and of the
other ordinary superiority measures can be derived using expression (5). In
particular, an mBR estimator of $\Delta(w;\theta)$ is

$\Delta^{*}(w;\theta^{*})=\Delta(w;\theta^{*})-B^{*}(w;\theta^{*})\,.$

where

$B^{*}(w;\theta)=\frac{1}{2}{\mathop{\rm
trace}}\left\\{i(\theta)^{-1}\nabla\nabla^{\top}\Delta(w;\theta)\right\\}\,,$

is the first term in the right-hand side of expression (5). Computing
$\Delta^{*}(w;\theta^{*})$ requires only the mBR estimator $\theta^{*}$ that
can be obtained using the procedures in Section 8, the corresponding estimated
category probabilities at $(w^{\top},1)$ and $(w^{\top},0)$, the matrix
$i(\theta^{*})^{-1}$, and the hessian
$\nabla\nabla^{\top}\Delta(w;\theta^{*})$. All these quantities, except
$\nabla\nabla^{\top}\Delta(w;\theta^{*})$, are readily available or can be
readily computed once the model has been estimated using mBR, as is done, for
example, in Section 3 for ACL models and in Kosmidis (2014) for cumulative
link models. For specific ordinal-response models, the hessian
$\nabla\nabla^{\top}\Delta(w;\theta)$ can be analytically obtained with some
algebraic effort. For example, if $\pi_{j}(w,z)$ is based on cumulative link
models one can work with the expressions for the derivatives of
$\Delta(w;\theta)$ in Agresti and Kateri (2017, Web appendix A).
Alternatively, a very accurate approximation of
$\nabla\nabla^{\top}\Delta(w;\theta^{*})$ can be obtained for general models
using a ready implementation of $\Delta(w;\theta)$ and numerical
differentiation routines, like the ones provided in the numDeriv R package
(Gilbert and Varadhan, 2019). This is the route that the ordinal_superiority()
method of the brglm2 R package takes.

Due to the equivariance properties of mBR estimation in Section 6.2 under
linear transformations, mBR estimators of $\gamma(w;\theta)$,
$\Delta^{\dagger}(\theta)$, and $\gamma^{\dagger}(\theta)$ are readily
obtained by replacing $\Delta(w;\theta)$ by $\Delta^{*}(w;\theta^{*})$ in
expressions (16) and (17). Wald-type inferences about the mBR estimators of
the ordinal superiority measures can be constructed as proposed in Agresti and
Kateri (2017, Section 5), using the mBR estimates of the ordinal superiority
measures along with estimated standard errors obtained using the delta method,
based on $i(\theta^{*})$, and numerical gradients.

######  Example 5:

mBR for ordinal superiority measures from ACL model In order to assess the
finite sample properties of the mBR estimator of ordinal superiority scores in
ACL models, we consider the example in Christensen (2019, Section 4.3), where
the bitterness ratings “2”, “3”, and “4” in Table 2 are merged into a single
rating “2-4”. Like the PO version of the cumulative logit model (see
Christensen, 2019, Section 4.8), the ML estimates for $\alpha_{2-4}$ and
$\beta_{1}$ for the PO version of the ACL model (13) are $+\infty$ and
$-\infty$ respectively. The mBR estimates of
$\theta=(\alpha_{1},\alpha_{2-4},\beta_{1},\beta_{2})^{\top}$, on the other
hand, take the finite values $\theta^{*}=(-1.247$, $5.331$, $-3.291$,
$-1.181)^{\top}$. If $\gamma(w,\theta)$ is the ordinal superiority measure for
temperature setting $w$ ($w=0$ for cold and $w=1$ for warm), and $z$ indicates
contact ($z=1$) or not ($z=0$) of the juice with the skin, then
$\gamma(0,\theta^{*})=0.594$ and $\gamma(1,\theta^{*})=0.575$, indicating that
there is almost $60\%$ chance of higher bitterness ratings when there is
contact of the juice with the skin.

We simulate $10,000$ samples from the PO version of the ACL model at
$\bar{\theta}$, and we compute $\gamma(w,\hat{\theta})$ and
$\gamma^{*}(w,\theta^{*})$ for each sample. The simulation-based estimates of
the finite-sample relative biases of $\gamma(w,\hat{\theta})$ are $0.84\%$ and
$1.56\%$ for $w=0$ and $w=1$, respectively. As expected, the mBR version
$\gamma^{*}(w,\theta^{*})$ is found to have smaller finite-sample relative
biases at $0.13\%$ and $-0.02\%$ for $w=0$ and $w=1$, respectively. The
corresponding percentages of underestimation are $48.48\%$ and $44.69\%$ for
$\gamma(w,\hat{\theta})$, and $52.12\%$ and $51.14\%$ for
$\gamma^{*}(w,\theta^{*})$. Hence, in this case, mBR also results in
improvements in median bias. Finally, both estimators appear to perform
satisfactorily in terms of Wald-type inferences based on them. The coverage
probability of the nominally $95\%$ Wald-type confidence intervals based on
$\gamma(w,\hat{\theta})$ are $94.8\%$ ($w=0$) and $94.6\%$ ($w=1$), and
$94.7\%$ ($w=0$) and $95.1\%$ ($w=1$) for those based on
$\gamma^{*}(w,\theta^{*})$.

##  10 Supplementary material

The supplementary material consists of three scripts that replicate all the
numerical results and graphics reported in the paper, and is available at
https://ikosmidis.com/files/bracl_supplementary_v0.2.zip. The results are
exactly reproducible in R version 4.1.2, and with the following packages: VGAM
version 1.1-5 (Yee, 2021), tibble 3.1.6 (Müller and Wickham, 2021), dplyr
1.0.7 (Wickham et al., 2021), ggplot2 3.3.5 (Wickham, 2016), colorspace 2.0-2
(Zeileis et al., 2020), and ordinal 2019.12-10 (Christensen, 2019), brglm2
0.8.2 (Kosmidis, 2021), enrichwith 0.3.1 (Kosmidis, 2020), and
detectseparation 0.2 (Kosmidis and Schumacher, 2021).

##  11 Acknowledgements

The author greatly appreciates the constructive discussions with Alan Agresti
and Anestis Touloumis during the Challenges for Categorical Data Analysis 2018
Workshop in Aachen University on the equivalence between ACL and BCL models,
which informed this work. Ioannis Kosmidis has been partially supported by The
Alan Turing Institute under the EPSRC grant EP/N510129/1.

## References

  * Agresti (2002) Agresti, A. (2002).  Categorical Data Analysis (2nd ed ed.).  New York: Wiley-Interscience.
  * Agresti (2010) Agresti, A. (2010).  Analysis of ordinal categorical data (2nd ed.).  Hoboken, NJ: Wiley.
  * Agresti and Kateri (2017) Agresti, A. and M. Kateri (2017).  Ordinal probability effect measures for group comparisons in multinomial cumulative link models.  Biometrics 73(1), 214–219.
  * Albert and Anderson (1984) Albert, A. and J. Anderson (1984).  On the existence of maximum likelihood estimates in logistic regression models.  Biometrika 71(1), 1–10.
  * Baker (1994) Baker, S. G. (1994).  The multinomial-Poisson transformation.  Journal of the Royal Statistical Society: Series D (The Statistician) 43(4), 495.
  * Bull et al. (2002) Bull, S. B., C. Mak, and C. M. T. Greenwood (2002).  A modified score function estimator for multinomial logistic regression in small samples.  Computational Statistics and Data Analysis 39, 57–74.
  * Christensen (2019) Christensen, R. H. B. (2019).  ordinal—regression models for ordinal data.  R package version 2019.12-10. https://CRAN.R-project.org/package=ordinal.
  * Cox and Hinkley (1974) Cox, D. R. and D. V. Hinkley (1974).  Theoretical Statistics.  London: Chapman & Hall Ltd.
  * Di Caterina and Kosmidis (2019) Di Caterina, C. and I. Kosmidis (2019).  Location-adjusted Wald statistics for scalar parameters.  Computational Statistics & Data Analysis 138, 126–142.
  * Firth (1993) Firth, D. (1993).  Bias reduction of maximum likelihood estimates.  Biometrika 80(1), 27–38.
  * Gilbert and Varadhan (2019) Gilbert, P. and R. Varadhan (2019).  numDeriv: Accurate Numerical Derivatives.  R package version 2016.8-1.1.
  * Gioia et al. (2021) Gioia, V., E. C. Kenne Pagui, and A. Salvan (2021).  Median bias reduction in cumulative link models.  Communications in Statistics, 1–17.
  * Haldane (1955) Haldane, J. (1955).  The estimation of the logarithm of a ratio of frequencies.  Annals of Human Genetics 20, 309–311.
  * Heinze and Schemper (2002) Heinze, G. and M. Schemper (2002, August).  A solution to the problem of separation in logistic regression.  Statistics in Medicine 21(16), 2409–2419.
  * Kenne Pagui et al. (2017) Kenne Pagui, E. C., A. Salvan, and N. Sartori (2017).  Median bias reduction of maximum likelihood estimates.  Biometrika 104(4), 923–938.
  * Klotz (1966) Klotz, J. H. (1966).  The Wilcoxon, Ties, and the Computer.  Journal of the American Statistical Association 61(315), 772–787.
  * Konis (2007) Konis, K. (2007).  Linear Programming Algorithms for Detecting Separated Data in Binary Logistic Regression Models.  Ph. D. thesis, University of Oxford.
  * Kosmidis (2007) Kosmidis, I. (2007).  Bias reduction in exponential family nonlinear models.  Ph. D. thesis, University of Warwick.
  * Kosmidis (2014) Kosmidis, I. (2014).  Improved estimation in cumulative link models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76(1), 169–196.
  * Kosmidis (2020) Kosmidis, I. (2020).  enrichwith: Methods to enrich list-like R objects with extra components.  R package version 0.3.1.
  * Kosmidis (2021) Kosmidis, I. (2021).  brglm2: Bias Reduction in Generalized Linear Models.  R package version 0.7.2.
  * Kosmidis and Firth (2009) Kosmidis, I. and D. Firth (2009).  Bias reduction in exponential family nonlinear models.  Biometrika 96(4), 793–804.
  * Kosmidis and Firth (2011) Kosmidis, I. and D. Firth (2011).  Multinomial logit bias reduction via the poisson log-linear model.  Biometrika 98(3), 755–759.
  * Kosmidis and Firth (2021) Kosmidis, I. and D. Firth (2021).  Jeffreys-prior penalty, finiteness and shrinkage in binomial-response generalized linear models.  Biometrika 108(1), 71–82.
  * Kosmidis et al. (2020) Kosmidis, I., E. C. Kenne Pagui, and N. Sartori (2020).  Mean and median bias reduction in generalized linear models.  Statistics and Computing (to appear) 30, 43–59.
  * Kosmidis and Schumacher (2021) Kosmidis, I. and D. Schumacher (2021).  detectseparation: Detect and Check for Separation and Infinite Maximum Likelihood Estimates.  R package version 0.2.
  * Lesaffre and Albert (1989) Lesaffre, E. and A. Albert (1989).  Partial separation in logistic discrimination.  Journal of the Royal Statistical Society. Series B (Methodological) 51(1), 109–116.
  * Lyles et al. (2012) Lyles, R. H., Y. Guo, and S. Greenland (2012).  Reducing bias and mean squared error associated with regression-based odds ratio estimators.  Journal of Statistical Planning and Inference 142(12), 3235–3241.
  * McCullagh (1980) McCullagh, P. (1980).  Regression models for ordinal data.  Journal of the Royal Statistical Society. Series B (Methodological) 42, 109–142.
  * McCullagh (2018) McCullagh, P. (2018).  Tensor methods in statistics (2nd ed.).  Mineola, NY: Dover Publications.
  * Mehrabi and Matthews (1995) Mehrabi, Y. and J. N. S. Matthews (1995).  Likelihood-based methods for bias reduction in limiting dilution assays.  Biometrics 51, 1543–1549.
  * Müller and Wickham (2021) Müller, K. and H. Wickham (2021).  tibble: Simple Data Frames.  R package version 3.1.6.
  * Pace and Salvan (1997) Pace, L. and A. Salvan (1997).  Principles of Statistical Inference from a Neo-Fisherian Perspective.  World Scientific.
  * Palmgren (1981) Palmgren, J. (1981).  The Fisher information matrix for log linear models arguing conditionally on observed explanatory variables.  Biometrika 68(2), 563.
  * Peterson and Harrell (1990) Peterson, B. and J. Harrell, Frank E. (1990).  Partial proportional odds models for ordinal response variables.  Applied Statistics 39, 205–217.
  * R Core Team (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Randall (1989) Randall, J. H. (1989).  The analysis of sensory data by generalised linear model.  Biometrical Journal 7, 781–793.
  * Silvapulle (1981) Silvapulle, M. J. (1981).  On the existence of maximum likelihood estimators for the binomial response models.  Journal of the Royal Statistical Society. Series B (Methodological) 43(3), 310–313.
  * Venables and Ripley (2002) Venables, W. N. and B. D. Ripley (2002).  Modern Applied Statistics with S (Fourth ed.).  New York: Springer.  ISBN 0-387-95457-0.
  * Wickham (2016) Wickham, H. (2016).  ggplot2: Elegant Graphics for Data Analysis.  Springer-Verlag New York.
  * Wickham et al. (2021) Wickham, H., R. François, L. Henry, and K. Müller (2021).  dplyr: A Grammar of Data Manipulation.  R package version 1.0.7.
  * Yee (2021) Yee, T. W. (2021).  VGAM: Vector Generalized Linear and Additive Models.  R package version 1.1-5.
  * Zeileis et al. (2020) Zeileis, A., J. C. Fisher, K. Hornik, R. Ihaka, C. D. McWhite, P. Murrell, R. Stauffer, and C. O. Wilke (2020).  colorspace: A toolbox for manipulating and assessing colors and palettes.  Journal of Statistical Software 96(1), 1–49.


# Lightwave-controlled relativistic plasma mirrors

Marie Ouillé1,2 Jaismeen Kaur1 Zhao Cheng1 Stefan Haessler1
<EMAIL_ADDRESS>Rodrigo Lopez-Martens1 rodrigo.lopez-
<EMAIL_ADDRESS>1Laboratoire d’Optique Appliquée, Institut
Polytechnique de Paris, ENSTA-Paris, Ecole Polytechnique, CNRS, 91120
Palaiseau, France,  2Ardop Engineering, 6 All. Annie Fratellini, 33140
Villenave-d’Ornon, France

###### Abstract

We report on attosecond-scale control of high-harmonic and electron emission
from plasma mirrors driven by relativistic-intensity near-single-cycle
lightwaves at kHz repetition rate. By controlling the waveform of the intense
light transient, we reproducibly form a sub-cycle temporal intensity gate at
the plasma mirror surface, leading to the observation of extreme ultraviolet
spectral continua, characteristic of isolated attosecond pulse generation. We
also observe the correlated emission of a waveform-dependent relativistic
electron beam, paving the way towards fully lightwave-controlled dynamics of
relativistic plasma mirrors.

The generation of attosecond light pulses has enabled the real-time
observation of the fastest electronic processes in matter. Although new
generation schemes have emerged [1, 2], the currently preferred route to
produce attosecond pulses remains High-Harmonic Generation (HHG) from the gas
phase [3]. The mechanism is now well understood, the technology is mature (as
indicated by the emergence of turn-key, tabletop HHG beamlines), and it has
been successfully used in several pioneering applications unveiling previously
unobservable attosecond-scale dynamics [4, 5, 6]. However, it is intrinsically
limited to low intensity ($\sim 10^{15}\,\rm W/cm^{2}$) driving pulses and the
maximum laser-to-XUV conversion efficiency is currently of the order of
$10^{-4}$ [7, 8]. To circumvent those limitations, HHG from relativistic
plasma mirrors has long been recognized as a promising route [9] owing to the
predicted percent-level conversion efficiencies [10] (although the measured
efficiencies are $\sim 10^{-4}$ so far [11, 12]), their compatibility with
relativistic-intensity ($>10^{18}\,\rm W/cm^{2}$) driving lasers, as well as
the production of chirp-free attosecond pulses [13] with high spatial
coherence [14, 15].

The nonlinear interaction leading to the emission of XUV light being repeated
once (in the case of surface harmonics) or twice (in the case of gas
harmonics) per optical cycle, trains of attosecond pulses spaced by one or
half a laser period, respectively, are usually generated. To gate the process
and generate isolated attosecond pulses (IAP), different techniques based on
largely available multi-cycle drivers have been explored, such as polarization
gating [16, 17], wavefront rotation [18, 19] and two-color gating [20, 21,
22]. However, these methods are costly in terms of achievable focused laser
intensity, and it is much more efficient to gate the process using single-
cycle duration laser pulses with controlled carrier-envelope phase (CEP).
Routinely used for IAP generation from the gas phase [23] for well over a
decade now, this gating technique has thus far only been applied to
relativistic plasma mirrors through the binning and post-processing of data
generated with randomly fluctuating laser waveforms [24, 12, 25].

HHG from plasma mirrors is also known to be accompanied by the emission of
electrons [26, 27], injected at specific phases in the reflected laser pulse
[28], in the form of attosecond bunches [29]. Typically, in the case of multi-
cycle driving pulses, depending on their initial energy and time of injection
during the pulse, these electrons can either be isotropically scattered due to
the laser ponderomotive potential, leading to the observation of a doughnut-
shaped electron distribution centered about the specular direction far away
from the interaction point, or undergo direct vacuum laser acceleration (VLA)
and further gain kinetic energy before being ejected along the polarization
plane, leading to the observation of an electron density peak slightly offset
from the specular direction [30].

In this letter, we present experimental results obtained on the correlated
high-harmonic and electron emission from relativistic plasma mirrors driven by
CEP-controlled sub-1.5-cycle laser pulses: for the first time, via CEP gating,
we are able to reproducibly generate XUV continua and control relativistic
electron beaming with attosecond precision, thereby paving the way towards the
reliable generation of synchronized isolated attosecond light and electron
bunches from relativistic lightwave-driven plasma mirrors.

Experiments were conducted using the Salle Noire laser driver [31], developed
in-house at Laboratoire d’Optique Appliquée, which routinely delivers multi-mJ
energy sub-4 fs pulses (FWHM) at a central wavelength $\lambda_{0}\approx 780$
nm, corresponding to less than 1.5 optical cycles, with a temporal contrast
ratio $>10^{10}$ up to 15 ps and better than $10^{7}$ up to 3 ps before the
pulse peak. A shot-to-shot CEP stability of $\approx$ 380 mrad rms can be
maintained for extended ($\sim$ hour) periods of time [32]. The CEP, which is
always inherently measured relative to an unknown offset value, can thus be
varied at will while keeping all other laser parameters identical.

Figure 1: Schematic drawing of the experimental setup. Insert (a) shows the
experimentally measured focused beam profiles on target for both the Prepulse
(gray-color scale) and Main-pulse (multi-color scale). Insert (b) shows a
typical angle-resolved XUV spectrogram (top: vertically integrated spectrum).
Insert (c) shows a typical electron angular charge distribution recorded with
the wide-angle Lanex screen. The dashed-black line indicates the specular
direction.

The experimental setup [33] is shown in fig. 1. The laser pulses are made to
interact with plasma mirrors generated at the surface of rotating solid fused
silica targets [34]. We use a time-delayed Prepulse, obtained by picking-off a
small fraction ($\approx$ 4 $\%$ of the energy) of the Main-pulse, in order to
initiate and control plasma expansion at the target surface. Both p-polarized
pulses are tightly focused on target at an incidence angle of $55^{\circ}$ by
a f/1.3, $30^{\circ}$ off-axis parabola. As shown in the fig. 1(a), the
Prepulse and Main-pulse are focused down to $\approx 10\>\mu$m and $\approx
1.5\times 1.7\>\mu$m FWHM spot size, respectively. Considering the Main-pulse
energy on target to be approximately 2.5 mJ and assuming there are no spatio-
temporal couplings, this corresponds to a relativistic peak intensity of
$I_{\mathrm{peak}}=1.8\times 10^{19}\,\rm W/cm^{2}$ i.e. a normalized vector
potential
$a_{0}=\lambda_{0}[\mu\rm{m}]{\sqrt{I_{\rm{peak}}\rm{[W\,cm^{-2}]}/1.37\times
10^{18}}}\approx 2.8$.

The spatio-spectral XUV content of the reflected beam (photon energies ranging
from 10 eV to 40 eV) is measured in the specular direction
($\phi_{az}=55^{\circ}$ with respect to the target-normal) with a flat-field
concave XUV grating and an MCP-phosphor screen assembly imaged on a CCD camera
(fig. 1(b)). A scintillating screen (Lanex) covered with a 13-µm aluminum foil
and imaged onto another CCD camera allows simultaneous detection of electrons
with kinetic energy $>150\>$keV in the range,
$30^{\circ}<\phi_{az}<48^{\circ}$ with respect to the target normal. A 0.5-mm
pinhole and a pair of permanent magnets can be inserted in front of the Lanex
screen in order to measure the energy spectra of electrons emitted at
$\phi_{az}\approx 38^{\circ}$ and $\theta_{pol}\approx 0^{\circ}$, which
allows simultaneous XUV detection. Alternatively, a larger Lanex screen can be
inserted to image the full electron angular charge distribution covering
$18^{\circ}<\phi_{az}<72^{\circ}$ with respect to the target normal (fig.
1(c)), which now blocks XUV detection and prevents any direct correlation
measurement with XUV generation. Nevertheless, different experimental reports
confirm that electron emission and HHG are correlated in the relativistic
regime [26, 27, 33]. No electron spectra could be recorded in this
configuration due to the size of the Lanex assembly.

Figure 2: Vertically-integrated HHG spectra as a function of relative CEP for
Prepulse lead times of: -1 ps (a), 1 ps (b), and 2 ps (c). Above each frame
are two line-out spectra measured for opposite laser waveforms ($\pi$ CEP
offset). Each data acquisition is averaged over 100 consecutive laser shots
(100 ms). To the right of each frame, we show the temporal auto-correlation of
the HHG emission in the shown spectral range. The white dashed line marks the
position of the satellite peak near 1 laser period. The scans were made in
succession while the CEP-locking loop remained active. The relative CEPs thus
correspond to the same driving waveforms in each panel. Figure 3: Simultaneous
acquisition of XUV spectra (top panel), electron beam profiles (middle panel),
and electron spectra (bottom panel) for different relative CEP values. Each
data acquisition is averaged over 100 consecutive laser shots (100 ms). The
red dot in the middle panel indicates the point of measurement of the electron
spectrum. The dashed boxes indicate the two opposite waveforms in the electron
push-pull scenario: at the top, the angle-integrated HHG spectrum is shown for
the two opposite driving waveforms and corresponding push-pull scenarios at
absolute CEP values, $\psi_{\textrm{CEP}}=\pi/2$ and $3\pi/2$. The formation
of an XUV continuum can clearly be resolved in the case of a single push-pull
sequence for a given relative CEP (step size $<0.5\>$rad).

In order to estimate the surface plasma gradient scale length
$L_{g}=L_{0}+c_{s}\tau$, where $L_{0}\approx 0.01\lambda_{0}$ is the estimated
contribution of the Main-pulse finite temporal contrast, and $\tau$ the delay
between the Prepulse and the Main-pulse, we experimentally measured the plasma
expansion speed $c_{s}$ using spatial domain interferometry [35, 33] and found
$c_{s}=22\ $nm/ps. In fig. 2, we show the vertically-integrated XUV spectra
recorded as a function of relative CEP for three different values of $\tau$:
-1 ps, 1 ps, and 2 ps, corresponding to respectively $L_{g}\approx L_{0}$,
$0.04\lambda_{0}$ and $0.07\lambda_{0}$. We observe clear CEP and
$L_{g}$-dependent spectral structures of the harmonics. Assuming a nearly flat
spectral phase [13], these can be linked to the periodicity of the
corresponding attosecond temporal profiles, apparent in the temporal auto-
correlations obtained as the modulus of the Fourier transform of each HHG
spectrum in the shown spectral range and also shown in fig. 2. For
$L_{g}\approx L_{0}$ (fig. 2a), the slower spectral modulation on top of the
harmonic orders probably results from the presence of both coherent wake
emission [36] and relativistic HHG. Their different emission times [13] then
create 2 attosecond pulses per $\approx 2.5$-fs laser period, which appears as
an additional peak near 1.3 fs in the auto-correlations. As the plasma
gradient softens to $L_{g}\approx 0.04\lambda_{0}$ (fig. 2b), relativistic HHG
becomes dominant, and only the satellite peak near the $\approx 2.5$-fs laser
period remains in the auto-correlation. Its CEP-dependent shift over $\approx
0.3\>$fs likely results from a CEP-dependent plasma denting evolution over the
pulse duration [12]. Finally, in panel (c), there’s interestingly no clear
temporal shift of the satellite peak (i.e., the harmonic modulation period is
nearly constant).

Finally, we find that a 2-ps prepulse delay ($L_{g}\approx 0.07\lambda_{0}$,
fig. 2c) maximizes the spectral extent of the harmonics (from $\tau=3\>$ps on,
both harmonic signal yield and spectral extension drop), which is close to the
gradient optimum of $0.1\lambda_{0}$ typically observed for relativistic HHG
[37]. Only for this optimum plasma gradient scale length do we observe clear
XUV spectral continua for a given relative CEP value ($0\>$rad in fig.2c),
which is characteristic of IAP generation. Clearly modulated spectra are
observed when the relative CEP is offset by $\pi$. Interestingly, the
corresponding $\approx 2.5$-fs satellite peak in the temporal auto-
correlations shifts by no more than $\approx 0.1\>$fs before falling below the
noise level.

We repeated these measurements at the optimal density gradient ($L_{g}\approx
0.07\lambda_{0}$), but this time simultaneously recording both the XUV
spectrum and the electron angular charge distribution as a function of
relative CEP (fig. 3). The locked CEP was varied randomly to eliminate any
systematic experimental bias. Here again, an XUV spectral continuum is clearly
produced for a certain relative CEP value ($0~{}$rad), whereas the most
modulated spectrum appears upon adding a $\approx\pi$ phase offset
($-3~{}$rad). This $\pi$ periodicity is reminiscent of our previous
observations via CEP-tagging and is very well reproduced by numerical
simulations predicting compression of the lightwave into a sub-0.4 fs IAP,
potentially containing a third of the incident pulse energy [25], without any
spectral filtering required. Physically, this CEP-dependence can be understood
through a three-step process [28, 25]: (1) push, (2) pull, and (3) emission of
XUV light, which occurs only once in the optimal-gating-case of
$\psi_{\textrm{CEP}}=3\pi/2$ and twice in the case of
$\psi_{\textrm{CEP}}=\pi/2$ (see fig.3). We also clearly observe, for the
first time, a closely correlated waveform-dependence of the electron emission:
the highest electron charge is obtained in correlation with the XUV continuum
in which case a clear charge density peak appears close to the specular
direction, while the lowest charge is produced by the opposite laser waveform
($\pi$ CEP offset), in which case no density peak is observed.

Figure 4: Angular electron charge distribution scanned over randomized CEP
values, measured at the $\tau=2\>$ps delay between the Main-pulse and
Prepulse. Each electron image is averaged over 100 consecutive laser shots
(100 ms). The dashed black line indicates the specular direction at
$\phi_{az}=55^{\circ}$ w.r.t. target normal.

Measuring the electron energy spectrum (point of measurement indicated by the
red dot in the middle panel of fig. 3) confirms the clear influence of the CEP
on the electron acceleration mechanism. The dashed boxes in fig. 3 highlight
the acquired electron spectra under similar conditions for the two opposite
laser waveforms: In the right box, for a given relative CEP ($0\>$rad), we are
in the scenario in which the push-pull mechanism only occurs once, as
indicated by the production of an XUV continuum. Interestingly, this
corresponds to the production of an electron energy spectrum with a single
broad peak at $\approx$ 1 MeV. As the CEP is shifted by $\pi$ (left box,
CEP$=-3~{}$rad), a more modulated electron spectrum is observed, with the
appearance of a shoulder at lower energies. These features could be the
indication of the phase-locked injection and subsequent VLA of either one (in
the optimal-gating case of $\psi_{\textrm{CEP}}=3\pi/2$) or two (in the case
of $\psi_{\textrm{CEP}}=\pi/2$) attosecond relativistic electron bunches.
These preliminary results should be confirmed in a more detailed experimental
study. Unfortunately, our setup only allowed us to measure the spectrum at the
periphery of the electron beam, as indicated by the red dot in the middle
panel in fig. 3. We expect higher energies in the beam center, in particular
in the optimal-gating scenario (right box), and intend to confirm this in
future experiments with a new design for the diagnostics. Our observation is
nonetheless the first experimental evidence of CEP-controlled emission of
relativistic electrons from plasma mirrors.

Fig. 4 shows the full electron beam measured under similar conditions with the
wide-angle Lanex screen, as the relative CEP is cycled through $2\pi$. The
electron angular charge distribution exhibits the same CEP dependence as in
the previous measurements (middle panel in fig. 3), with maximum charge peak
appearing at a relative CEP of $0~{}$rad and a minimally charged beam
(approximately half the total integrated charge) when the CEP is shifted by
$\pi$. It is interesting to note that although the distribution and direction
of the electron beam do not change significantly with CEP, there is, however,
a small sign of ponderomotive blowout around the specular direction in the
case of relative CEP of 3 rad, which would indicate that when
$\psi_{\textrm{CEP}}=\pi/2$ one of the electron bunches is not injected with
the adequate phase-energy combination for efficient VLA. We plan to measure
the spatially-resolved electron spectrum across the full beam in order to
visualize CEP-dependent fine structure that could provide more insight into
the attosecond-scale electron acceleration dynamics.

In conclusion, we report the first experimental evidence of the CEP dependence
and control of correlated high-harmonic and relativistic electron emission
from plasma mirrors driven at kHz repetition rate by controlled near-single-
cycle lightwaves. In particular, we can produce XUV spectral continua on
demand by forming a singular sub-cycle temporal intensity gate at the plasma
mirror surface, which in light of previously reported experimental
observations of the very high spatio-temporal quality of the emitted
attosecond pulses [13, 14, 15] represents a major step towards the realization
of isolated XUV attosecond pulses with record-intensity at kHz repetition
rate, which we aim at verifying in the future.

###### Acknowledgements.

This work was supported by the Agence Nationale pour la Recherche
(ANR-11-EQPX-005-ATTOLAB, ANR-14-CE32-0011-03 APERO, ANR-22-CE30-0005-01
BANDITO); European Research Council (ERC FEMTOELEC 306708, ERC ExCoMet
694596); LASERLAB-EUROPE (H2020-EU.1.4.1.2. grant agreement ID 654148), Région
Ile-de-France (SESAME 2012-ATTOLITE).

## References

  * Hassan  _et al._ [2016] M. T. Hassan, T. T. Luu, A. Moulet, O. Raskazovskaya, P. Zhokhov, M. Garg, N. Karpowicz, A. M. Zheltikov, V. Pervak, F. Krausz, and E. Goulielmakis, Optical attosecond pulses and tracking the nonlinear response of bound electrons, Nature 530, 66 (2016).
  * Franz  _et al._ [2024] P. Franz, S. Li, T. Driver, R. R. Robles, D. Cesar, E. Isele, Z. Guo, J. Wang, J. P. Duris, K. Larsen, J. M. Glownia, X. Cheng, M. C. Hoffmann, X. Li, M.-F. Lin, A. Kamalov, R. Obaid, A. Summers, N. Sudar, E. Thierstein, Z. Zhang, M. F. Kling, Z. Huang, J. P. Cryan, and A. Marinelli, Terawatt-scale attosecond x-ray pulses from a cascaded superradiant free-electron laser, Nature Photonics 10.1038/s41566-024-01427-w (2024).
  * Hentschel  _et al._ [2001] M. Hentschel, R. Kienberger, C. Spielmann, G. A. Reider, N. Milosevic, T. Brabec, P. Corkum, U. Heinzmann, M. Drescher, and F. Krausz, Attosecond metrology, Nature 414, 509 (2001).
  * Drescher  _et al._ [2002] M. Drescher, M. Hentschel, R. Kienberger, M. Uiberacker, V. Yakovlev, A. Scrinzi, T. Westerwalbesloh, U. Kleineberg, U. Heinzmann, and F. Krausz, Time-resolved atomic inner-shell spectroscopy, Nature 419, 803 (2002).
  * Sansone  _et al._ [2010] G. Sansone, F. Kelkensberg, J. F. Pérez-Torres, F. Morales, M. F. Kling, W. Siu, O. Ghafur, P. Johnsson, M. Swoboda, E. Benedetti, F. Ferrari, F. Lépine, J. L. Sanz-Vicario, S. Zherebtsov, I. Znakovskaya, A. L’Huillier, M. Y. Ivanov, M. Nisoli, F. Martín, and M. J. J. Vrakking, Electron localization following attosecond molecular photoionization, Nature 465, 763 (2010).
  * Calegari  _et al._ [2014] F. Calegari, D. Ayuso, A. Trabattoni, L. Belshaw, S. D. Camillis, S. Anumula, F. Frassetto, L. Poletto, A. Palacios, P. Decleva, J. B. Greenwood, F. Martín, and M. Nisoli, Ultrafast electron dynamics in phenylalanine initiated by attosecond pulses, Science 346, 336 (2014), https://www.science.org/doi/pdf/10.1126/science.1254061 .
  * Takahashi  _et al._ [2013] E. J. Takahashi, P. Lan, O. D. Mücke, Y. Nabekawa, and K. Midorikawa, Attosecond nonlinear optics using gigawatt-scale isolated attosecond pulses, Nature Communications 4, 2691 (2013).
  * Nayak  _et al._ [2018] A. Nayak, I. Orfanos, I. Makos, M. Dumergue, S. Kühn, E. Skantzakis, B. Bodi, K. Varju, C. Kalpouzos, H. I. B. Banks, A. Emmanouilidou, D. Charalambidis, and P. Tzallas, Multiple ionization of argon via multi-xuv-photon absorption induced by 20-gw high-order harmonic laser pulses, Phys. Rev. A 98, 023426 (2018).
  * Heissler  _et al._ [2012] P. Heissler, R. Hörlein, J. M. Mikhailova, L. Waldecker, P. Tzallas, A. Buck, K. Schmid, C. M. S. Sears, F. Krausz, L. Veisz, M. Zepf, and G. D. Tsakiris, Few-cycle driven relativistically oscillating plasma mirrors: A source of intense isolated attosecond pulses, Phys. Rev. Lett. 108, 235003 (2012).
  * Mikhailova  _et al._ [2012] J. M. Mikhailova, M. V. Fedorov, N. Karpowicz, P. Gibbon, V. T. Platonenko, A. M. Zheltikov, and F. Krausz, Isolated attosecond pulses from laser-driven synchrotron radiation, Phys. Rev. Lett. 109, 245005 (2012).
  * Rödel  _et al._ [2012] C. Rödel, D. an der Brügge, J. Bierbach, M. Yeung, T. Hahn, B. Dromey, S. Herzer, S. Fuchs, A. G. Pour, E. Eckner, M. Behmke, M. Cerchez, O. Jäckel, D. Hemmers, T. Toncian, M. C. Kaluza, A. Belyanin, G. Pretzler, O. Willi, A. Pukhov, M. Zepf, and G. G. Paulus, Harmonic generation from relativistic plasma surfaces in ultrasteep plasma density gradients, Phys. Rev. Lett. 109, 125002 (2012).
  * Jahn  _et al._ [2019] O. Jahn, V. E. Leshchenko, P. Tzallas, A. Kessel, M. Krüger, A. Münzer, S. A. Trushin, G. D. Tsakiris, S. Kahaly, D. Kormin, L. Veisz, V. Pervak, F. Krausz, Z. Major, and S. Karsch, Towards intense isolated attosecond pulses from relativistic surface high harmonics, Optica 6, 280 (2019).
  * Chopineau  _et al._ [2021] L. Chopineau, A. Denoeud, A. Leblanc, E. Porat, P. Martin, H. Vincenti, and F. Quéré, Spatio-temporal characterization of attosecond pulses from plasma mirrors, Nature Physics 17, 968 (2021).
  * Dromey  _et al._ [2009] B. Dromey, D. Adams, R. Hörlein, Y. Nomura, S. G. Rykovanov, D. C. Carroll, P. S. Foster, S. Kar, K. Markey, P. McKenna, D. Neely, M. Geissler, G. D. Tsakiris, and M. Zepf, Diffraction-limited performance and focusing of high harmonics from relativistic plasmas, Nature Physics 5, 146 (2009).
  * Leblanc  _et al._ [2017] A. Leblanc, S. Monchocé, H. Vincenti, S. Kahaly, J.-L. Vay, and F. Quéré, Spatial properties of high-order harmonic beams from plasma mirrors: A ptychographic study, Phys. Rev. Lett. 119, 155001 (2017).
  * Baeva  _et al._ [2006] T. Baeva, S. Gordienko, and A. Pukhov, Relativistic plasma control for single attosecond x-ray burst generation, Phys. Rev. E 74, 065401 (2006).
  * Yeung  _et al._ [2015] M. Yeung, J. Bierbach, E. Eckner, S. Rykovanov, S. Kuschel, A. Sävert, M. Förster, C. Rödel, G. G. Paulus, S. Cousens, M. Coughlan, B. Dromey, and M. Zepf, Noncollinear polarization gating of attosecond pulse trains in the relativistic regime, Phys. Rev. Lett. 115, 193903 (2015).
  * Wheeler  _et al._ [2012] J. A. Wheeler, A. Borot, S. Monchocé, H. Vincenti, A. Ricci, A. Malvache, R. Lopez-Martens, and F. Quéré, Attosecond lighthouses from plasma mirrors, Nature Photonics 6, 829 (2012).
  * Vincenti and Quéré [2012] H. Vincenti and F. Quéré, Attosecond lighthouses: How to use spatiotemporally coupled light fields to generate isolated attosecond pulses, Phys. Rev. Lett. 108, 113904 (2012).
  * Zhang  _et al._ [2020] Y. X. Zhang, S. Rykovanov, M. Shi, C. L. Zhong, X. T. He, B. Qiao, and M. Zepf, Giant isolated attosecond pulses from two-color laser-plasma interactions, Phys. Rev. Lett. 124, 114802 (2020).
  * Smith  _et al._ [2021] R. E. Smith, E. White, P. M. Gilmore, B. Dromey, and M. Yeung, Double optical gating of high order harmonics from plasma surfaces, New Journal of Physics 23, 093010 (2021).
  * Blanco and Flores-Arias [2017] M. Blanco and M. Flores-Arias, Frequency gating to isolate single attosecond pulses with overdense plasmas using particle-in-cell simulations, Opt. Express 25, 13372 (2017).
  * Goulielmakis  _et al._ [2008] E. Goulielmakis, M. Schultze, M. Hofstetter, V. S. Yakovlev, J. Gagnon, M. Uiberacker, A. L. Aquila, E. M. Gullikson, D. T. Attwood, R. Kienberger, F. Krausz, and U. Kleineberg, Single-cycle nonlinear optics, Science 320, 1614 (2008), https://www.science.org/doi/pdf/10.1126/science.1157846 .
  * Kormin  _et al._ [2018] D. Kormin, A. Borot, G. Ma, W. Dallari, B. Bergues, M. Aladi, I. B. Földes, and L. Veisz, Spectral interferometry with waveform-dependent relativistic high-order harmonics from plasma surfaces, Nature Communications 9 (2018).
  * Böhle  _et al._ [2020] F. Böhle, M. Thévenet, M. Bocoum, A. Vernier, S. Haessler, and R. Lopez-Martens, Generation of XUV spectral continua from relativistic plasma mirrors driven in the near-single-cycle limit, Journal of Physics: Photonics 2, 034010 (2020).
  * Chopineau  _et al._ [2019] L. Chopineau, A. Leblanc, G. Blaclard, A. Denoeud, M. Thévenet, J.-L. Vay, G. Bonnaud, P. Martin, H. Vincenti, and F. Quéré, Identification of coupling mechanisms between ultraintense laser light and dense plasmas, Phys. Rev. X 9, 011050 (2019).
  * [27] S. Haessler, M. Ouillé, J. Kaur, M. Bocoum, F. Böhle, D. Levy, L. Daniault, A. Vernier, J. Faure, and R. Lopez-Martens, High-harmonic generation and correlated electron emission from relativistic plasma mirrors at 1 kHz repetition rate, Ultrafast Science 2022, 9893418.
  * Thévenet  _et al._ [2016] M. Thévenet, H. Vincenti, and J. Faure, On the physics of electron ejection from laser-irradiated overdense plasmas, Physics of Plasmas 23, 063119 (2016), https://doi.org/10.1063/1.4954822 .
  * Zhou  _et al._ [2021] C. Zhou, Y. Bai, L. Song, Y. Zeng, Y. Xu, D. Zhang, X. Lu, Y. Leng, J. Liu, Y. Tian, R. Li, and Z. Xu, Direct mapping of attosecond electron dynamics, Nature Photonics 15, 216 (2021).
  * Thévenet  _et al._ [2016] M. Thévenet, A. Leblanc, S. Kahaly, H. Vincenti, A. Vernier, F. Quéré, and J. Faure, Vacuum laser acceleration of relativistic electrons using plasma mirror injectors, Nature Physics 12, 355 (2016).
  * Ouillé  _et al._ [2020] M. Ouillé, A. Vernier, F. Böhle, M. Bocoum, A. Jullien, M. Lozano, J.-P. Rousseau, Z. Cheng, D. Gustas, A. Blumenstein, P. Simon, S. Haessler, J. Faure, T. Nagy, and R. Lopez-Martens, Relativistic-intensity near-single-cycle light waveforms at kHz repetition rate, Light Science & Applications 9, 2047 (2020).
  * Huijts  _et al._ [2022] J. Huijts, L. Rovige, I. A. Andriyash, A. Vernier, M. Ouillé, J. Kaur, Z. Cheng, R. Lopez-Martens, and J. Faure, Waveform control of relativistic electron dynamics in laser-plasma acceleration, Phys. Rev. X 12, 011036 (2022).
  * Kaur  _et al._ [2023] J. Kaur, M. Ouillé, D. Levy, L. Daniault, A. Robbes, N. Zaïm, A. Flacco, E. Kroupp, V. Malka, S. Haessler, and R. Lopez-Martens, High repetition rate relativistic laser–solid–plasma interaction platform featuring simultaneous particle and radiation detection, Rev. Sci. Instr. 94, 113002 (2023).
  * Borot  _et al._ [2014] A. Borot, D. Douillet, G. Iaquaniello, T. Lefrou, P. Audebert, J.-P. Geindre, and R. Lopez-Martens, High repetition rate plasma mirror device for attosecond science, Review of Scientific Instruments 85, 013104 (2014), https://doi.org/10.1063/1.4860035 .
  * Bocoum  _et al._ [2015] M. Bocoum, F. Böhle, A. Vernier, A. Jullien, J. Faure, and R. Lopez-Martens, Spatial-domain interferometer for measuring plasma mirror expansion, Optics Letters 40, 3009 (2015).
  * Quéré  _et al._ [2006] F. Quéré, C. Thaury, P. Monot, S. Dobosz, Martin, Ph., J. P. Geindre, and P. Audebert, Coherent Wake Emission of High-Order Harmonics from Overdense Plasmas, Physical Review Letters 96, 125004 (2006).
  * Kahaly  _et al._ [2013] S. Kahaly, S. Monchocé, H. Vincenti, T. Dzelzainis, B. Dromey, M. Zepf, P. Martin, and F. Quéré, Direct observation of density-gradient effects in harmonic generation from plasma mirrors, Phys. Rev. Lett. 110, 175001 (2013).