# Planning a method for covariate adjustment in individually-randomised
trials: a practical guide

T.P.Tim P. Morris  A.S.A. Sarah Walker  E.J.Elizabeth J. Williamson  I.R.Ian
R. White  MRC Clinical Trials Unit at UCL, London, UK  Department of Medical
Statistics, LSHTM, London, UK

###### Abstract

Background It has long been advised to account for baseline covariates in the
analysis of confirmatory randomised trials, with the main statistical
justifications being that this increases power and, when a randomisation
scheme balanced covariates, permits a valid estimate of experimental error.
There are various methods available to account for covariates.

Methods We consider how, at the point of writing a statistical analysis plan,
to choose between three broad approaches: direct adjustment, standardisation
and inverse-probability-of-treatment weighting (IPTW), which are in our view
the most promising methods. Using the GetTested trial, a randomised trial
designed to assess the effectiveness of an electonic STI (sexually transmitted
infection) testing and results service, we illustrate how a method might be
chosen in advance and show some of the anticipated issues in action.

Results The choice of approach is not straightforward, particularly with
models for binary outcome measures, where we focus most of our attention. We
compare the properties of the three broad approaches in terms of the quantity
they target (estimand), how a method performs under model misspecification,
convergence issues, handling designed balance, precision of estimators,
estimation of standard errors, and finally clarify some issues around handling
of missing data.

Conclusions We conclude that no single approach is always best and explain why
the choice will depend on the trial context but encourage trialists to
consider the three methods more routinely.

Covariate adjustment,

Estimands,

Standardisation,

Inverse probability of treatment weighting,

Randomised controlled trials,

Clinical Trials,

Missing data,

###### keywords:

Methodology

## Background

Randomised controlled trials are designed to estimate average treatment
effects. This article considers the handling of covariates in the analysis of
individually randomised trials. By covariates, we mean measurements on
participants recorded at baseline that are thought to be prognostic. Typical
examples are age at randomisation and disease severity at randomisation.

Confounding does not affect trials that are properly randomised, since
confounding is a systematic bias and any imbalances in covariates in
randomised trials are due to chance, which introduces non-systematic error
that is reflected in inference[1]. However, it remains desirable from a
statistical perspective to account for covariates in the analysis of a trial.
There are two principal reasons for this.

First, covariate adjustment is desirable because, if a covariate predicts
outcome, accounting for its effect on outcome will improve power to detect a
treatment effect[2, 3, 4] unless none of the covariates in a model are
prognostic[2]. This is sometimes explained as ‘accounting for chance
imbalance’, though we view this differently: a large imbalance (as might be
expected in an observational study) can make inference less precise;
adjustment gains power by acknowledging chance covariate balance.

Second, adjustment is desirable because, if a covariate is balanced by the
randomisation scheme, for example by using stratified block randomisation or
minimisation, adjustment is necessary to obtain a valid estimate of
experimental error. An unadjusted analysis mistakenly assumes that chance
imbalance in covariates could have occurred – an extremely useful property for
covariates that were unmeasured. However, when randomisation has been balanced
according to measured covariates, the fact that imbalance could not have
occurred must be acknowledged in our statistical inference[5]. This was
pointed out by Fisher almost a century ago but is frequently ignored in
practice[6, 7].

A premise of this article is therefore that accounting for covariates in the
analysis of trials is a good idea. We are concerned with how trial
statisticians and clinical investigators should agree on a method before
seeing the data. This perspective is taken to help statisticians working on
trials make sensible, informed decisions when writing a statistical analysis
plan.

We will focus most heavily on binary outcome data since this is where some of
the issues are most acute. When the outcome measure is continuous and the
analysis aims to estimate a difference in means, some – though not all – of
the considerations of this paper become redundant. In particular,
collapsibility of the mean difference means that none of the discussion around
collapsibility is relevant.

For adjustment to be worthwhile, the covariates to be included in the analysis
must be prognostic: adjustment for a covariate that is not prognostic is
essentially equivalent to an unadjusted analysis, though it can lose power in
small samples. It is not the purpose of this article to discuss which
covariates to choose. The Committee for Proprietary Medicinal Products Points
to consider on adjustment for baseline covariates document gives some
guidance[8]. It does caution against approaches that select covariates most
strongly associated with the outcome in the trial; however, subsequent work
(for example [9, 10]) has shown that such a procedure can be pre-specified in
a principled manner.

Part of the motivation for this article is that statistical research papers
have frequently recommended covariate adjustment due to improvements in power,
with little thought given to the implications of different adjustment methods
(TPM being one culpable author). The approach generally favoured in clinical
trials is direct adjustment using an outcome regression model. This is so
prevalent that reviews of practice have not needed to discuss which methods
were actually used[11, 7]. This article works through the implications and
aims to contrast its properties with two other methods of adjustment better
known in the epidemiology literature.

## Methods

### Motivating examples: the GetTested trial

The GetTested trial was designed to assess the effectiveness of an Internet-
accessed sexually transmitted infection (STI) testing and results service
(chlamydia, gonorrhoea, HIV, and syphilis) on STI testing uptake and STI cases
diagnosed.

Briefly, 2,072 participants were recruited in the London boroughs of Southwark
and Lambeth. Participants were randomised to an invitation to use an internet-
based sexual health service (intervention) or a standard test from a walk-in
sexual health clinic (control).

Treatment allocation involved a minimisation procedure balancing for gender
(male/female/trans), number of sexual partners in the 12 months before
randomisation (categorised as 1, 2+; note that one-or-more was part of the
eligibility criteria) and sexual orientation (men who have sex with men vs.
other groups). Each of the covariates was weighted equally when determining
marginal imbalance and intervention or control were assigned at random with
80% probability of assignment to the favoured arm.

The two outcomes of principal interest were both binary. The first was whether
participants took an STI test within six weeks of randomisation, with the
control arm proportion anticipated to be around 10%. The second outcome, and
the outcome for which sample size was calculated, was STI diagnosis (following
a positive STI test) within six weeks of randomisation, with control arm
proportion anticipated to be around 0.6%. The primary analysis was planned to
account for the following covariates: gender, age, number of sexual partners
in the 12 months before randomisation (10 categories where the final category
is $>10$), sexual orientation (four categories) and ethnicity (five
categories)[12]. The chosen covariates were all assumed to be prognostic.

### Three broad approaches to accounting for covariates

Three broad approaches to covariate adjustment are described below. We outline
the generic procedures for estimation in this section and discuss their
properties in the next section. Note that the estimand targeted is discussed
there.

#### Direct adjustment

‘Direct adjustment’ refers to fitting an outcome regression model including
terms for randomised treatment $Z$ (an indicator equal to 1 if assigned to
treatment and 0 if assigned to control) and covariates $X$, but no interaction
between $Z$ and $X$. The treatment effect and its standard error are estimated
directly as the treatment coefficient from the model. This might be done, for
example, using a generalised linear model, which may be standard (with
canonical link function) or non-standard, or a Cox model[13, 14].

#### Standardisation

Standardisation fits an outcome regression model including $Z$ and $X$,
possibly including interactions with $Z$, and then standardises the results by
summing or integrating over the distribution of covariates observed in the
trial. One intuitive way to achieve this is by making predictions if all
participants were assigned to the intervention arm and then if all
participants were assigned to the control arm and forming a suitable contrast
of the predictions[15, 16, 17]. Standardisation builds on outcome regression
by fixing the summary of the treatment effect but allowing flexible modelling
to estimate it; for example, a risk ratio can be estimated using logistic
regression. Standardisation is sometimes termed marginalisation or
G-computation, and is implemented in Stata’s margins command and R’s stdReg
package[18].

#### Inverse probability of treatment weighting

Inverse probability of treatment weighting (IPTW) involves fitting two models.
The first uses participants’ covariate values to predict the probability of
being randomised to the arm they were in fact randomised to. The second then
fits a simple model contrasting the treatment arms, weighted according to the
inverse probability of treatment estimated by the first model. This
effectively creates a weighted trial sample (pseudo-population) in which both
trial arms have the same distribution of observed covariate values. This may
sound odd since, if we have used simple randomisation, we know that the ‘true’
weight for all participants is identical, so trying to predict it with
covariates appears futile[19]. However, the goal is not to obtain an estimate
of this probability but instead to either reweight to a better balanced trial
or acknowledge the balance observed.

The method is computationally the same as using propensity scores but the
justification is different (just as the justification for covariate adjustment
is different in randomised trials vs. cohort studies).

#### Comments on methods

There are many possible specific ways to implement each of the three broad
approaches. For example, Tsiatis and colleagues used a form of standardisation
but the within-arm outcome models are defined to yield ‘as good predictions as
possible without concerns over bias’. However, for the purposes of this
article, we will generally refer to or use simple implementations of the three
broad approaches. For standardisation with binary outcomes we will (by
default) use a logistic regression working model fitting main effects of
covariates to produce predictions. For IPTW we will use a logistic regression
including only the main effects of covariates to model the probability of
randomised group given covariates. For direct adjustment, we will use
generalised linear models with a link function that permits parameter
estimation on the scale of the summary measure of interest, for example a
binomial generalised linear model with identity link function for a risk
difference.

Besides the three broad approaches there exist hybrid methods; in particular,
estimators based on the influence function target marginal summaries and
include both a model for treatment (as does IPTW) and a model for outcome (as
does standardisation). Further, since IPTW relates to propensity scores, this
suggests other estimators, such as matching rather than weighting, would be
possible. However, the idea of discarding data from randomised individuals
makes it unattractive. We do not focus further on these approaches.

## Results

We now work through some properties that are important to consider when
writing a statistical analysis plan. These properties are mainly – but not
solely – statistical. We also report results of applying the different methods
to the GetTested trial.

### Summary measure of the estimand: marginal or conditional

The topic of estimands has become increasingly prominent in clinical trials
since the publication of the ICH E9(R1) addendum[20]. The addendum lists five
attributes of an estimand: the treatment condition of interest, target
population, outcome variable, handling of intercurrent events and population-
level summary. Here, we focus on the population-level summary (the later
section on handling missing data will focus on the target population). Table 1
lists some population-level summaries commonly used in clinical trials and, in
particular, notes whether conditional and marginal estimands coincide.

It is not our purpose to argue for any particular summary measure but it would
be remiss to pass over how the choice should be made. Some statisticians
assume that the correct approach must be to choose a measure as a parameter of
a model that might have generated the data. Permutt argues that the choice of
scale should be ‘linear in utility’: a hypothetical value of treatment effect
of 2 should be twice as attractive as 1 whatever the potential outcome on
control.[21] Others argue for a measure that can be easily interpreted. We
regard the first approach as misled even if the true model were known and the
second as too strict, since it rules out any relative measure. The choice of
measure should be a trial-specific tradeoff between ease of interpretation,
close relation to average patient benefit and potential transportability to
other settings.

The following discussion relates to non-collapsible summaries[22]. For readers
not familiar with this term, we provide a numerical example in table 2 and
note . The reader should suppose that the frequencies given are ‘true’ in the
sense that, had we recruited a very large sample size, the cells of table 2
would contain exact multiples of the frequencies shown.

Consider a trial in a condition that includes participants from two measurable
strata, A and B, which have a substantial effect on prognosis. The trial team
recruits 40 participants – 20 from each stratum – and, within strata,
randomises 10 to intervention and 10 to control. In stratum A, the odds of
dying on the control arm is $5/5=1$. In stratum B, the outlook is far more
favourable, with the odds of dying on the control arm just $1/9$. Despite
these differences, the treatment effect (a conditional odds ratio), is 9 in
each stratum. We might have conducted a trial recruiting patients from just
one of the strata or from both. If we put all 40 people together, as shown at
the right-most block of table 2 (hence the term marginal), the odds of dying
in the control arm becomes $(1+5)/(5+9)=3/7$ and the marginal odds ratio is
$5.4$.

At first sight this is astonishing! Both strata have odds ratios equal to 9
but, when put together, the odds ratio changes. It is not a weighted average
of the within-stratum odds ratios. Treatment is exactly balanced within
strata, so this is not due to imbalance. Neither is it effect modification,
since the log odds ratio is identical in both strata.

The general phenomenon is known as ‘non-collapsibility’, which describes the
relationship between the marginal and conditional summary measure: the true
marginal odds ratio is attenuated towards 1 compared with the conditional odds
ratio. It occurs because the average of the logit is not the logit of the
average. While the odds ratio is non-collapsible, the risk ratio and risk
difference are collapsible. For insights into why non-collapsibility occurs,
see Daniel, Zhang and Farewell[22]; see also [23] and [24] for more technical
discussions, particularly into the relation between collapsibility and
confounding.

Rather than collapsing our strata into ‘both’, as in the right-hand panel of
table 2, it is possible to adjust for strata and recover the conditional odds
ratio of 9 (using logistic regression adjusted for strata, or stratified
Mantel–Haenszel).

Note that non-collapsibility is not ‘bias’ as sometimes supposed, but a case
of different estimands. The odds ratio of 9 is a within-stratum or conditional
estimand formed by comparing the effect of treatment within a stratum,
relating to the question ‘what would be the odds ratio for treatment comparing
people in the same stratum?’ The odds ratio of 5.4 is a between-strata or
‘marginal’ estimand formed by comparing the effect of treatment for groups
made up of half stratum A and half B, relating to the question ‘what would be
the odds ratio in a population made up of half stratum A and half stratum B’.

Had the notional trial recruited a different proportion of patients from each
stratum, the conditional odds ratio would have remained 9 but the marginal
odds ratio would not have remained 5.4. For example, suppose that in stratum A
we had recruited three times as many people (so that the numbers in that row
are all multiplied by three). The marginal odds ratio would then be 6. Its
true value depended on the proportion of participants in each stratum. In
general, a marginal summary depends on the covariate distribution when we have
non-collapsibility. This is uncomfortable; after all, regression models
condition on covariate values rather than modelling their distribution, but we
see that the distribution nonetheless matters.

This discomfort may lead us to conclude that the conditional odds ratio is
obviously preferable. This is misled for two reasons. First, suppose there
were a second covariate to stratify on table 2 that also strongly predicted
outcome but it was unknown (or at least unmeasured). Then true value of the
conditional odds ratio still depends on the distribution of a covariate, but
since the covariate is unmeasured its distribution is unknown: we cannot know
what we are marginalising over. Second, this example is contrived such that
the conditional odds ratios within strata were identical. When there is effect
modification on the scale of the population summary measure, the true value of
the conditional measure will also depend on the distribution of observed
covariates.

Either a marginal or conditional estimand may be desirable and this depends on
context. For example, a patient may wish to know ‘what would happen if someone
similar to me were to choose this intervention vs. not?’ Meanwhile, for policy
makers, the average difference a drug would make if offered to a group of
people might be of more interest, though they might equally wish to know about
the effect for specific groups. Interestingly, marginal estimands appear to be
favoured for causal inference from observational data: Hernán and Robins
define a population causal effect as ‘a contrast of any functional of the
marginal distributions of counterfactual outcomes under different actions or
treatment values’ (emphasis added)[15]. By this definition, the within-stratum
odds ratio of 9 does not appear to target a population causal estimand, since
it is a contrast of conditional distributions, though the word marginal is a
choice rather than necessary to the definition.

In non-inferiority studies, where the null is a non-zero difference between
arms on some scale, it is frequently argued that the margin of non-inferiority
is easier to understand, define and interpret on a marginal than a conditional
scale. This is in line with our own experiences in collaborations.

The choice of marginal or conditional estimand is clearly not simple: the true
value of the estimand will depend either on the distribution of observed
covariates (always marginal and sometimes conditional), on which covariates
are conditioned-on in the model (conditional), and further on the distribution
of omitted prognostic covariates (both). Note that these aspects have
implications for: 1) the quantities being combined in meta-analysis; and 2)
apparent heterogeneity in meta-analysis, though we will not comment further on
these points here.

For non-inferiority studies with a non-collapsible summary measure, it is
worth noting which scale the non-inferiority margin is defined on. Suppose the
margin is specified as a marginal hazard ratio, then the corresponding non-
inferiority margin on the scale of the conditional hazard ratio is further
from 1. If this fact is forgotten and a conditional hazard ratio is estimated
without changing the margin, we could expect to lose power compared with an
unadjusted analysis (since this targets a marginal estimand).

In terms of the three broad methods considered, direct adjustment always
targets a conditional summary measure; standardisation typically – but not
necessarily – targets a marginal summary measure; and IPTW always targets a
marginal summary measure. Standardisation and IPTW are rarely used in trials
but receive more attention in the epidemiological literature. One
justification given for this seems be the notion that a ‘target trial’ would
always target a marginal summary[15]. Ironically, trials which do adjust for
covariates tend to use direct adjustment and so target a conditional summary.

### Convergence

Having defined an estimand, we require an estimator to compute an estimate.
For many estimators, parameter estimation proceeds through some iterative
technique. In maximum likelihood estimation, for example, an algorithm is used
to find parameters that maximize the likelihood of the data. This involves
finding parameter values that maximize a function. Once an algorithm has found
a maximum, it is said to have ‘converged’. It sometimes happens that the
algorithm fails to converge to a maximum or that the maximum to which is
converges is local (that is, a small bump rather than the true maximum) or not
unique. This is clearly an issue.

In the analysis of randomised trials, non-convergence tends to occur for one
of two reasons: first, problems that occur with certain models (e.g.
generalised linear models with binomial outcome distribution and identity- or
log-link function); second, including too many parameters compared with the
effective sample size (e.g. fitting fixed centre intercepts with few
participants per centre)[25]. Essentially, the observed data are not
consistent with a model that fits within the given constraints.

We are concerned with choosing a procedure for analysis prior to seeing data
so would be unwise to jeopardise the analysis by choosing a procedure that may
not converge[26, 27]. While it may be possible to specify a backup procedure,
it would need to target the same estimand. This may raise the question ‘why
not specify the backup procedure as the first choice?’ One good reason may be
due to lower power.

Direct adjustment and standardisation may involve fitting different outcome
models. For example, to estimate a risk difference, direct adjustment may use
a generalised linear model with binomial outcome with identity link function.
In terms of convergence, this would be a risky plan. It is possible to
estimate an adjusted risk difference using other methods. This could use
standardisation following estimation through a logistic regression model,
which comes with the guarantee of converging to a unique maximum. A popular
technique for estimating the risk ratio without incurring convergence problems
is to use a poisson model with robust standard errors[28], where convergence
is likely due to the canonical link function.

IPTW involves specifying a model for treatment $P(Z~{}|~{}X)$ to estimate
weights. This may be any model for binary data (regardless of the trial
outcome type). Due to randomisation, it is always true that
$P(Z~{}|~{}X)=P(Z)$. This means that the model for a binary treatment would
not be misspecified regardless of how covariates $X$ are modelled, provided
parameters were not constrained to be wrong. Allocation ratios in trials are
most frequently 1:1 but rarely more extreme than 1:2. For 1:1 allocation, the
‘outcome’ proportion in the treatment model will be approximately 50%, and
never near 0 or 1, and its distribution given covariates is random, making
‘separation’ (or ‘perfect prediction’) unlikely[29]. All this means that the
treatment model has a good chance of converging. The subsequent outcome model
has no covariates and so is certain to converge.

IPTW therefore seems to be the safest broad approach if convergence is
anticipated to be an issue, while standardisation may mitigate possible issues
associated with direct adjustment (given the same estimand such as a risk
ratio).

### Efficiency/precision/power

A key reason to account for covariates in the analysis is to increase power.
Note that for non-collapsible summary measures it is wrong to attempt to
compare precision of marginal and conditional estimators but in general it is
possible to compare power when the null is zero difference, since
collapsibility then holds.[22] Marginal adjusted estimators have been shown to
be more efficient than marginal unadjusted estimators. [30, 31]

Because adjustment separates the effect of a treatment from the effects of
covariates, we can typically infer the effect of a treatment with a little
more precision, though it is possible to lose precision in small samples with
non-prognostic covariates. It is therefore usually desirable to use an
efficient method of accounting for covariate effects, or the potential gains
in power may not be fully realised.

While it is sometimes argued that weighting estimators are inefficient,
Williamson, White and Forbes showed that, in the trials context with a
continuous outcome measure, IPTW is asymptotically as efficient as direct
adjustment[19]. Any ‘inefficiency’ of IPTW tends to arise due to extreme
weights, just as the variance reduction it achieves is a result of similar
estimated weights for all individuals. As with convergence, thinking about the
weighting model makes clear that this will not be a problem when using this
method in the analysis of trials. Note that a closely-related method, overlap
weighting, has recently been shown to be more efficient in finite samples and
is worth further consideration for covariate adjustment[32, 33].

### Handling covariates balanced by design

There are many methods of balancing covariates at the design stage. The most
popular seem to be stratified blocks and minimisation[7]. When a covariate-
balancing method is used, it is necessary to account for the randomisation
scheme in the analysis, or the estimated standard error for the treatment
effect will be biased upwards, producing confidence intervals that are too
wide (meaning they have greater than $1-\alpha$ coverage) and miscalibrated
p-values[5].

For some intuition, suppose a trial uses stratified blocks, with
stratification by a single, binary, prognostic covariate. That covariate will
then always be distributed equally across the randomised groups (provided each
block is completed). It is then impossible for any difference seen to be due
to this covariate. Effectively, the slight variability in the treatment effect
due to possible chance imbalances is eliminated, since imbalances can never
occur under this design. An analysis that ignores this systematic balance will
assume that imbalance in a covariate could have occurred by chance and
calculate a standard error accordingly. This would be too large, since an
imbalance in this particular covariate could not in fact have occurred[1].
Adjustment for the covariate separates the effect of covariate/s on outcome
from the effect of treatment on outcome, and this is acknowledged in
calculation of the standard error.

There is literature on this going back at least to Fisher, who seemed to
regard the point as obvious for analysis of variance in agricultural
experiments[6]. It is generally accepted that direct adjustment and
standardisation can be specified to provide valid standard errors. This also
holds for IPTW, though has not previously been commented on. Suppose we again
have a single binary covariate which is perfectly balanced across treatment
groups and estimate $\text{P}(Z=z~{}|~{}X=x)$. Then this probability will be
identical for every individual in the trial. Fitting the weighted regression
to contrast the treatment effect will then return an identical estimate to a
model that ignores the covariate. What is perhaps surprising is that the IPTW
estimator still has a smaller standard error than the unweighted model. The
variance formula, which does not say anything about the design, effectively
‘sees’ and acknowledges the balance after estimating the weights and rewards
itself accordingly[19]. It is also clear that this happens by analogy to
direct adjustment, since asymptotically the two methods have the same standard
errors.

### Variance estimation

Some variance formulas rely on approximations and some are asymptotic. For
direct adjustment based on maximum likelihood estimation, formulas are
available for all commonly used models. For standardisation, the standard
error from a fitted model is transformed using the delta method (asymptotic).
Due to non-linearity, this could lead to p-values and confidence intervals
that do not quite agree, if the p-value is taken from the model on the
estimation scale; one possibility is to use test-based confidence intervals
for measures with the same null. For example, suppose the outcome model is a
logistic regression and the summary measure of the treatment effect is a risk
difference. The logistic regression returns a test-statistic of $z$ for
treatment. A 95% confidence interval for the risk difference can then be
constructed by taking $\pm 1.96/z$ times the distance between the estimated
risk difference and 0, and adding the result to the estimated risk difference.

IPTW estimators use robust standard errors that acknowledge the estimation of
weights in the first step. These robust standard errors are asymptotically
valid but recent work has demonstrated that they can produce slight
undercoverage in small samples[34]. In this type of setting, non-parametric
bootstrap may be required, until a closed variance formula with a small-sample
correction has been developed.

While non-parametric bootstrap is a useful tool, we regard it as not being
ideal due to inherent Monte Carlo error – though it may sometimes be the only
option. Monte Carlo error can be made small with a suitably large number of
bootstrap repetitions for important analyses. An important but often neglected
point about the bootstrap is that the resampling procedure must mimic the
sampling used in the study itself. A simple bootstrap procedure invokes simple
randomisation, and will return upwardly biased standard errors or confidence
intervals that are too wide if the trial did not in fact use simple
randomisation. So, if the trial used blocked randomisation within strata, the
bootstrap procedure should be restricted to do the same, else it targets the
wrong sampling distribution. Ensuring that the bootstrap procedure mimics the
design actually used may be awkward if for example the trial used
minimisation.

### Misspecification of the covariate model

We consider misspecification of the mean function relating covariates to
outcome rather than misspecification more generally. We illustrate the ideas
using a continuous outcome, which lends itself to this visual explanation, but
expect similar results for other outcome types.

Consider a study with a single covariate $X_{i}$, randomised treatment $Z_{i}$
and the model that generates outcomes $Y_{i}$ is

$Y_{i}=\alpha+\theta Z_{i}+\gamma X_{i}^{2}.$ (1)

Note that there is no residual error here; the mean function determines the
outcome exactly. The analyst plans to fit to the observed data a model (which
is misspecified) with mean function

$y_{i}=\hat{\alpha}+\hat{\theta}z_{i}+\hat{\lambda}x_{i}.$ (2)

Suppose that this notional study is run and the observed $x$ among those
recruited is perfectly uniform on $(-0.5,0.5)$, as depicted by the first
horizontal grey bar at the top of figure 1. When the misspecified model is
fitted, $\hat{\lambda}=0$. Next consider a trial where observed $x_{i}$ values
are uniform on $(0,1)$ or $(0.5,1.5)$, also depicted by grey horizontal bars
in figure 1. It is now clear that when the analyst fits their model 2,
$\hat{\lambda}>0$. In the first case the sample correlation of $x$ with
$x^{2}$ is zero, but in second and third cases it is greater than zero. The
analyst’s adjustment for $x$ thus partially adjusts for $x^{2}$ despite the
model being misspecified. This will generally be true when a covariate
actually adjusted for is correlated with covariates not adjusted for.

The lower panel of figure 1 gives the estimated standard errors after linear
adjustment, showing that linear adjustment is always as efficient as no
adjustment.

When using direct adjustment using the data in figure 1, the model is
misspecified. Meanwhile, when using IPTW, the model used to form weights is by
definition correctly specified. Despite this, the two return nearly identical
results. The criteria for good specification of IPTW are slightly different
than usual: what matters is not the correctness of the specification of the
model for $Z~{}|~{}X$ but how well the model for $Z~{}|~{}X$ models the
predictors of $Y$; doing this better will result in a more suitable
‘rebalancing’. By attempting to balance $X$ instead of $X^{2}$, the covariate
will be well balanced at certain points but less so at (for example)
particularly high or low values of $X$.

### Handling data missing at random with the three adjustment approaches

Some data will inevitably be missing for some participants in the majority of
randomised trials. This is most often in the outcome (unless ‘missing’ is
somehow part of the outcome definition) but sometimes occurs in one or more
covariates. Meanwhile, data on randomised arm should never be missing. We will
first discuss the issues when outcomes are incomplete and then when covariates
are incomplete, along with some solutions for each of the three broad
approaches. As with any inference from incomplete data, it is important to
understand the mechanisms under which bias will and will not be introduced and
so we discuss missingness dependent on randomised arm, covariates and outcome
separately.

With missing outcome data, a good starting point is to consider implications
of missingness for the simplest analysis: including only those individuals
with observed outcomes (complete-case analysis). When the probability of data
being missing depends only on randomised arm, a complete-case estimator is
unbiased and efficient. When missingness depends on the outcome, complete-case
analysis and multiple imputation under missing-at-random are biased in
general; we will return to this point under sensitivity analysis.

When outcome missingness depends on the observed value/s of covariates, a
complete-case estimator may or may not be biased. If the covariate/s were not
adjusted for, data would be missing not at random, which is a more difficult
statistical problem; if the covariate/s causing missingness are adjusted for
in the analysis, this becomes a missing-at-random problem. It is sometimes
said that multiple imputation is not needed when outcomes are missing-at-
random. However, estimation based on the complete cases is unbiased for a
population represented by the complete cases. This estimand does not in
general equal the estimand that targets the population actually randomised,
unless the treatment effect is the same for these two populations; potentially
a strong assumption. Supposing that the all-randomised population is of
interest, a complete-cases estimator is then potentially biased, though the
magnitude of bias will be small in practice. For a worked numerical example
explaining this point, see the supplementary material.

If the aim of accounting for covariates in the analysis is simply to increase
precision or to estimate a conditional summary measure, and not to target the
all-randomised population, the remainder of this section and the appendices
can be skipped.

With outcomes missing according to this mechanism (depending on covariates),
performing multiple imputation by-randomised-arm and then analysing the
imputed data sets by any of the three approaches to adjustment targets the
all-randomised population[35]. For direct adjustment, this is straightforward
but there are subtleties in terms of statistical inference for the other two
approaches:

  * •

For standardisation, the question is whether to apply Rubin’s rules before or
after the standardisation step. Since approximate normality is more likely on
the estimation scale (for example log-odds) than the summary scale (for
example risk difference), this is likely to be the appropriate scale for
combining.

  * •

For IPTW, it is however possible that Rubin’s variance formula will be
inconsistent due to uncongeniality[36]. Further, attempting to use multiple
imputation may involve fitting the direct adjustment model – the first step of
standardisation – so using multiple imputation may imply that IPTW is not
needed. We would lose, for example, the advantages in terms of convergence.

There are alternatives to multiple imputation. Under covariate-dependent
missingness, standardisation can be applied to the all-randomised sample
rather than only those with complete outcome data. Meanwhile, IPTW can be
combined with inverse probability of missingness weighting, with missingness
predicted from covariates separately by randomised arm. These are not our
primary recommendations because we want a principal analysis that can be
readily extended to principled sensitivity analyses.

Missing covariate values are not inevitable in randomised trials but do
sometimes occur. As with missing outcome data, analysis based only on the
complete cases may inadvertently target a complete-cases population rather
than all randomised, and may be biased if the all-randomised population is
intended. Unlike with missing outcome data, discarding individuals with
observed data on treatment and outcome does not follow the intention-to-treat
principle and does throw away information. However, it is sometimes simple to
target the all-randomised population: for any method that targets a marginal
or collapsible summary measure, simple mean imputation (across arms, not
within) and the missing indicator method are generally appropriate
methods[37]. When using direct adjustment with a non-collapsible summary
measure, it is more difficult to deal with incomplete covariate data, and this
typically requires a correct model[37]: If covariates are missing not at
random, this can be very difficult. The message is then that more care is
required to collect all covariate data if direct adjustment is to be used with
a non-collapsible summary measure.

### Sensitivity analysis with outcomes missing-not-at-random

There are rarely cases where we know the true missingness mechanism. The
assumption of missing-at-random depending on randomised arm and covariates is
a convenient starting point but it is important to examine the extent to which
inferences are sensitive to alternative missingness mechanisms. This prompts
sensitivity analysis.[26] To obtain valid inference, missingness mechanisms
that represent departures from missing-at-random then need to be explicitly
invoked.

In this situation, we view multiple imputation as a general and convenient
framework for statistical inference under various departures for each of the
three broad approaches, though not the only one.[38] Suppose for example that
the planned approach to covariate adjustment was standardisation, and under
missing-at-random we planned to standardise to the all-randomised sample
(which is valid). There is no extension of this concept under missing-not-at-
random mechanisms. Multiple imputation by-arm is a convenient way to do
sensitivity analysis. It makes little sense to have a mismatch between the
primary analysis and sensitivity analyses other than the missingness mechanism
invoked, since we want to ensure that sensitivity of results are attributable
to the change of missingness mechanism rather than the change of method. It is
possible that multiple imputation under missing-at-random may have delivered a
slightly different result to all-randomised standardisation. Sensitivity
analyses are important enough that we regard coherence between the primary and
sensitivity analysis as worthy of consideration.

## Analyses of GetTested

This article focuses on planning a statistical analysis but it is nonetheless
instructive to consider some of the issues discussed when different approaches
are used. This helps to illustrate what may happen and prompts us to reflect
on how we might plan.

As described previously, our analysis of the GetTested trial is for two
outcome measures: any test, which occurred in $35\%$, and any diagnosis, which
occurred in $1.6\%$.[39] For both outcomes, two summary measures were of
interest to the investigators: the risk ratio and the risk difference. Recall
that there are five categorical covariates to adjust for, which in the direct
adjustment model use $18$ parameters in addition to the intercept and
randomised arm. Table 3 gives the estimated treatment effect and standard
error from various analyses estimating each measure on the two outcomes. For
the Any test and Any diagnosis outcomes, there were $612$ and $27$ events
respectively from $1,739$ observed outcomes ($324$ had missing outcomes).
Table 3 presents results for the complete-cases population. Results targeting
the all-randomised estimand are presented and discussed in appendix 2.
Covariate data were fully observed.

The results presented are intended to compare adjustment methods for the same
outcome and summary measure. For the risk ratio, two direct adjustment methods
are used: a log-binomial model and a poisson model. In the analyses presented,
the methods used included ‘main’ effects of covariates only. Of course,
interactions between covariates could have been included for any method, and
interactions between covariates and randomised arm could have been included
for standardisation. Given the low number of diagnosis events – both
anticipated and observed – including these interactions would have been
inadvisable for that outcome. For these illustrative analyses, missing outcome
data were assumed to be missing-at-random given covariates.

Three of the 14 analyses in table 3 failed to produce any sensible estimate.
Two of these instances were due to the use of an identity-link-function
binomial model to estimate a risk difference, leading to non-convergence,
which happened for both outcomes. This emphasises the point that it would have
been unwise to plan this as the adjustment model. The last was when directly
estimating the log risk ratio for any test using a poisson model. It did
produce an estimate, which was 541, which indicated separation of outcome and
a clearly untrustworthy estimate.

Of the methods that did converge, the estimated treatment effects and standard
errors tended to be similar across methods. The most notable difference in
estimates is the complete-cases log risk ratio for any diagnosis, where the
direct and standardisation analyses estimated a larger value than the IPTW
analysis. This turns out to be due to the covariate men who have sex with men,
which in the direct model based on poisson regression estimated extremely
large risk ratios (around 15). Removing this covariate from the two models
resulted in very similar estimates, close to the IPTW result shown.

## Discussion

We have compared the properties of three broad methods for estimating and
adjusted treatment effect: direct adjustment, standardisation and inverse
probability of treatment weighting. Our impression is that direct adjustment
is the most commonly used approach in clinical trials and that standardisation
and IPTW are better appreciated in observational epidemiology and warrant more
consideration by trialists. In particular, it is clear that many investigators
are interested in summarising the treatment effect as a risk ratio or risk
difference. Direct adjustment is notoriously unstable for both measures and so
an unadjusted estimate is frequently reported, which will be inefficient.
Standardisation and IPTW warrant more use in clinical trials.

Having discussed some properties of the three approaches, we provide table 4
for reference, listing some of the points to consider for each approach.

We have described the methods in general terms without giving recommendations
about specific implementations. In our analysis of GetTested, the
standardisation estimates were computed from the direct adjustment model
(logistic regression including main effects of covariate) but could have been
computed from a model including interactions among covariates and of
covariates with randomised group. Similarly, we used IPTW on the propensity
score, which was estimated using logistic regression.

Tension between methods is typically greatest with binary outcome data, so we
chose to focus on this. In particular, for continous data analysed using
linear models, none of the issues around collapsibility were relevant.
Covariate adjustment using analysis of covariance is then typically suitable
(direct regression of the outcome on treatment and covariates) unless
treatment–covariate interactions are to be modelled, when standardisation
would be the appropriate approach.

Binary outcomes are common in non-inferiority trials and the scale of the
population summary should match the scale on which the non-inferiority margin
is defined. In our experience, the margin is never defined as a conditional
odds ratio, suggesting that the use of standardisation or IPTW will be
necessary.

For time-to-event outcome data, the issues will depend on the chosen summary
measure, as with binary outcome data. In clinical trials, this is frequently
the hazard ratio, which is non-collapsible. A further important issue specific
to hazard ratios is that even if hazards are proportional on the conditional
scale they may not be proportional on the marginal scale, so an adjusted
marginal hazard ratio (estimated by standardisation or IPTW) may be an
inappropriate summary of the treatment effect, though this is also true of an
unadjusted hazard ratio.

It is possible to estimate the covariate-adjusted difference in survival
proportion at a time $t$ using standardisation or IPTW, which some trials are
beginning to do in statistical analysis plans, or the difference in
restricted-mean survival times[40]. A variance estimator has been derived for
standardisation[18], though not (to our knowledge) for IPTW.

## Conclusions

We hope that this work stimulates statisticians working in trial teams to
think carefully about adjustment methods, particularly by placing the estimand
– which requires clinical investigators’ input – first, followed consideration
of the more statistical aspects. None of the three approaches is always best
and the choice will depend on the trial context. Standardisation and IPTW are
largely unused in trials, but have many advantages which mean they warrant
routine consideration.

## Appendix 1: targeting the complete-case and all-randomised populations in
the presence of missing data

The received wisdom is that, when missingness depends only on covariates, a
complete-case estimator is unbiased and efficient. There is however a subtlety
to this, which impacts on our three adjustment methods: The estimand no longer
targets a population with the empirical distribution of covariates among those
randomised but among the complete cases. This is true unless the treatment
effect is identical within strata on the chosen summary scale, as we will
show.

As with non-collapsibility, the fact that the value of the estimand for all-
randomised and complete-cases may differ is not bias but a case of targeting
different estimands: the target population attribute of the estimands differ.
It is perhaps difficult to contrive an argument for the complete-case
estimand.

To demonstrate this point and outline some solutions, we will work through a
simple numerical example with covariate-dependent missingness, show that a
standard complete-case estimator changes the estimand, and note one possible
advantage of standardisation over the other methods in this respect. The
numerical example can be regarded as deterministic in the sense that
estimators that fail to recover the exact value of the estimand seen in the
true data do target a different estimand.

Consider a trial with binary $X$, $Y$ and $Z$, where the data can be
represented in a $2\times 2\times 2$ table. Table 5(a) gives the full data
from a notional randomised trial in which $2,000$ participants are recruited.
Of these, $1,000$ have the covariate $X=1$ and $1,000$ have $X=0$. Treatment
$Z$ is stratified by covariate $X$ so that $500$ participants have each
combination of values of $X$ and $Z$. Finally, table 5(c) gives the within-
stratum odds ratios, risk ratios and risk differences. Note that we have
ensured there is some effect modification by $X$ on each of these summary
scales, since we will not in general be able to choose a scale on which there
is no effect modification (particularly before seeing the data).

Suppose now that among those with $X=0$ exactly half are a complete case,
while among those with $X=1$ all are a complete case. The complete cases are
depicted in table 5(b). The probability of being a complete case depends only
on the covariate $X$, meaning that the frequencies in the cells of the
$2\times 2$ table for $X=1$ are exactly half what they were in table 5(a).
Table 5(c) thus correctly represents the true values of the within-stratum
summaries from table 5(b) and 5(a).

As previously, we do not dictate which summary of the treatment effect should
be used for the trial. We consider four: a risk ratio, a risk difference, a
conditional odds ratio, a marginal odds ratio. In these examples, conditional
odds ratios are estimated by logistic regression, while marginal odds ratios,
risk ratios and risk differences are estimated in two ways: first by
standardisation after logistic regression, and second by IPTW. The
standardisation analysis does not include an interaction of $X$ with $Z$ in
the logistic regression model, though it could do. Note that when two or more
approaches can be used for one summary measure (the marginal odds ratio, risk
ratio and risk difference are all estimated using standardisation and IPTW),
the estimates produced will be identical (see [15], chapters 12 and 13).

Table 6 gives results of analyses for all randomised and for the complete
cases. The key point is that the true value of the complete-cases estimand
differs from the true value of the all-randomised estimand for all summary
measures. This is because there is effect modification by $X$ on each scale
and the relative proportion with $X=1$ changes in the complete cases under
covariate-dependent missingness, so the relative contribution of each stratum-
specific effect to the overall effect changes.

It would seem that there is little to choose between the methods in terms of
handling of missing data. However, there are different modifications to the
adjustment methods that we can used to target the all-randomised estimand in
several cases.

We have until now talked about ‘complete cases’ without specifying whether $X$
or $Y$ is incomplete. With $X$ incomplete, the simple-mean-imputation or
missing indicator methods describd by White and Thompson can be used to target
the all-randomised estimand for marginal summary measures[37]. Moving beyond
complete cases is imperative here since the incomplete cases have important
information in observed $Y$. However, neither method is appropriate if the
summary measure of interest is conditional and non-collapsible and data are
missing-not-at-random, in which case there is no method to target the all-
randomised estimand besides correct modelling of the not-at-random missingness
mechanism.

With covariate-dependent missingness and outcome $Y$ incomplete, a simple
method can be used to target standardisation to the all-randomised estimand:

  1. 1.

Retain in the data all randomised individuals;

  2. 2.

Fit the estimation model for standardisation in the complete cases (possibly
including additional or different interactions between $Z$ and covariates);

  3. 3.

Standardise over the all-randomised distribution of $X$ (implemented in
Stata’s margins command with the noesample option).

The result of applying this method to the numerical example gives results
identical to the left hand column of table 6, indicating that the procedure
can target the all-randomised estimand. However, as seen in the analysis of
the GetTested study (see Appendix 2), this neat ‘trick’ is not fool-proof.

For direct adjustment and IPTW, estimators of the all-randomised estimand are
slightly less straightforward, but can be achieved by using multiple
imputation with a separate imputation model for each randomised arm or using
inverse probability of missingness weighting with the weighting model based on
randomised arm and covariates. It would seem most natural to pair direct
adjustment with by-arm multiple imputation, because both posit a model for the
outcome data. IPTW therefore pairs more naturally with inverse probability of
missingness weighting, where neither of the weighting models involve the value
of the outcome.

## Appendix 2: Issues with estimation for the all-randomised population in
GetTested

We now consider further results from re-analysis of GetTested. Table 3
considered an estimand for the complete-case population, but here we consider
the all-randomised population. To target this population, we retain the
assumption that outcomes are missing at random and then use a different method
for each adjustment approach.

For direct adjustment, we target the all-randomised population using multiple
imputation ‘by-arm’. Each outcome was multiply imputed separately (not
jointly), since they are always observed or missing simultaneously and so
there is not auxiliary information. The imputation procedure used a logistic
regression model, separately for each arm, to impute outcome, including main
effects of covariates. The imputation model for both outcomes was ‘augmented’
to handle separation of outcome given covariates[41]. Ten imputations were
used. Note that it is the separate imputation model by-arm that targets the
all-randomised population. Simply including randomised arm as a covariate in
the imputation model would not achieve this unless there were no
treatment–covariate interactions on the scale of the summary measure.

For standardisation, the all-randomised estimand was targeted by standardising
over the all-randomised population, as described in appendix 1.

For IPTW, inverse-probability-of-missingness weighting was used. The
probability of missingness was estimated using a logistic regression model for
missingness including all covariates that are present in the primary analysis,
and interactions with randomised arm. Since the two outcomes were always
missing simultaneously, estimated missingness probabilities were identical for
both outcomes. This approach is straightforward to implement due to the
monotone missingness pattern (baseline covariates were fully observed and
outcomes were missing on the same individuals).

Additional file 2 contains the Stata code used, which can be inspected to
resolve any ambiguity in the description of the three methods. Note that,
multiple imputation or inverse-probability-of-missingness weighting could have
been used for any of the three adjustment methods.

Table 7 shows the results from these analyses. The main point to note is that
some analyses which returned an estimate for the complete-case population (see
table 3) do not for the all-randomised. This turns out to be due to
collinearity or perfect prediction, described below. However, this affects the
methods in different ways, as we see below.

For direct adjustment, the multiple imputation step suffered from perfect
prediction due to collinearity and the imputation model for any diagnosis had
to be augmented[41]. The subsequent estimation of risk differences with an
identity-link binomial model failed to converge for both outcomes. For the any
diagnosis risk ratios, the log-link binomial model also failed to converge;
the poisson model did return an estimate but we do not regard this as a valid
all-randomised estimate because many parameters had to be dropped from the
imputation models in order for them to fit. For any test, there was severe
collinearity in the imputation model. Even after augmentation, only 916
individuals could be used in the imputation model. This again represents a
failure to return an estimate relevant to the all-randomised population.

For standardisation, no estimate was immediately returned for either outcome
or summary measure. This was due to collinearity in the logistic regression
models predicting outcome, which meant that predictions were not produced for
four transgender individuals. By omitting these individuals in the
standardisation step, an almost-all-randomised estimate could be obtained for
any test but not for any diagnosis (marked with asterisks in table 7).

IPTW returned estimates for both summary measures and both outcomes. The main
threat to it obtaining an estimate for the all-randomised population is
perfect prediction of missingness in the missingness model.

Given the above issues, there is litte to say in terms of comparing the
estimated log risk ratios obtained across approaches. When approaches did
return an estimate and standard error, they were very similar. It is hard to
say what differences can be attributed to: the populations targeted, the
specific assumptions about missing data, or the approach to modelling
assumptions specific to the approach to covariate adjustment. Further
investigation of methods targeting the all-randomised population is ongoing.
The ideal solution is of course to avoid missing outcome data as far as
possible during the conduct of a trial.

## Acknowledgements

We are grateful to Jonathan Bartlett, Richard Riley, Leanne McCabe, Brennan
Kahan, Andrew Althouse, Maarten van Smeden and Babak Choodari-Oskooei for
discussions relating to this work. Our acknowledgement of these individuals
does not imply their endorsement of this article.

## Funding

TPM, ASW and IRW are funded by the MRC grants MC_UU_12023/21 and
MC_UU_12023/29. ASW is also funded by the MRC grant MC_UU_12023/22. ASW is an
National Institute for Health Research Senior Investigator; as such, the views
expressed are those of the authors and not necessarily those of the NHS, the
NIHR, or the Department of Health. EJW was supported by the MRC Network of
Hubs for Trials Methodology HTMR Award MR/L004933/2/N96 and MRC project grant
MR/S01442X/1.

## Abbreviations

IPTW – inverse probability of treatment weighting
STI – sexually transmitted infection

## Availability of data and materials

Text for this section…

## Ethics approval and consent to participate

N/A

## Competing interests

Tim Morris consults for Kite Pharma, Inc. Ian White has provided consultancy
services or courses to Exelixis, AstraZeneca, GSK and Novartis, for which his
employer has received funding. Elizabeth Williamson declares personal income
from providing training to AstraZeneca.

## Consent for publication

N/A

## Authors’ contributions

TPM, ASW, EJW and IRW conceived of the article, planned the work and
interpreted the results. TPM drafted the article. All authors have approved
the submitted version.

## References

  * [1] S. Senn.  Seven myths of randomisation in clinical trials.  Statistics in Medicine, 32(9):1439–1450, 2013.
  * [2] B. C. Kahan, V. Jairath, C. J. Doré, and T. P. Morris.  The risks and rewards of covariate adjustment in randomized trials: an assessment of 12 outcomes from 8 studies.  Trials, 15(1):139+, 2014.
  * [3] A. V. Hernández, M. J. C. Eijkemans, and E. W. Steyerberg.  Randomized controlled trials with time-to-event outcomes: How much does prespecified covariate adjustment increase power?  Annals of Epidemiology, 16(1):41–48, 2006.
  * [4] E. L. Turner, P. Perel, T. Clayton, P. Edwards, A. V. Hernández, I. Roberts, H. Shakur, and E. W. Steyerberg.  Covariate adjustment increased power in randomized controlled trials: an example in traumatic brain injury.  Journal of Clinical Epidemiology, 65(5):474–481, 2012.
  * [5] B. C. Kahan and T. P. Morris.  Improper analysis of trials randomised using stratified blocks or minimisation.  Statistics in Medicine, 31(4):328–340, 2012.
  * [6] R. A. Fisher.  The arrangement of field experiments.  Journal of the ministry of agriculture of Great Britain, 33:83–94, 1926.
  * [7] B. C. Kahan and T. P. Morris.  Reporting and analysis of trials using stratified randomisation in leading medical journals: review and reanalysis.  BMJ (Clinical research edition), 345, 2012.
  * [8] Committee for Proprietary Medicinal Products.  Committee for proprietary medicinal products (CPMP) points to consider on adjustment for baseline covariates.  Statistics in Medicine, 23(5):701–709, 2004.
  * [9] A. A. Tsiatis, M. Davidian, M. Zhang, and X. Lu.  Covariate adjustment for two-sample treatment comparisons in randomized clinical trials: A principled yet flexible approach.  Statistics in Medicine, 27(23):4658–4677, 2008.
  * [10] L. B. Balzer, M. J. van der Laan, M. L. Petersen, and the SEARCH Collaboration.  Adaptive pre-specification in randomized trials with and without pair-matching.  Statistics in Medicine, 35(25):4528–4545, jul 2016.
  * [11] J. D. Ciolino, H. L. Palac, A. Yang, M. Vaca, and H. M. Belli.  Ideal vs. real: a systematic review on handling covariates in randomized controlled trials.  BMC Medical Research Methodology, 19(1), jul 2019.
  * [12] E. Wilson, C. Free, T. P. Morris, M. G. Kenward, J. Syred, and P. Baraitser.  Can internet-based sexual health services increase diagnoses of sexually transmitted infections (STI)? protocol for a randomized evaluation of an internet-based STI testing and results service.  JMIR Research Protocols, 5(1):e9, 2016.
  * [13] P. McCullagh and J. A. Nelder.  Generalized Linear Models.  Chapman and Hall, London, 2nd edition, 1989.
  * [14] D. R. Cox.  Regression models and life tables.  Journal of the Royal Statistical Society, series B, 34:187–220, 1972.
  * [15] M. A. Hernán and J. M. Robins.  Causal Inference: What If.  Boca Raton: Chapman & Hall/CRC, 2020.
  * [16] P. Cummings.  The relative merits of risk ratios and odds ratios.  Archives of Pediatrics & Adolescent Medicine, 163(5):438, 2009\.
  * [17] J. M. Snowden, S. Rose, and K. M. Mortimer.  Implementation of g-computation on a simulated data set: Demonstration of a causal inference technique.  American Journal of Epidemiology, 173(7):731–738, mar 2011.
  * [18] A. Sjölander.  Regression standardization with the R package stdReg.  European Journal of Epidemiology, 31(6):563–574, 2016.
  * [19] E. J. Williamson, A. Forbes, and I. R. White.  Variance reduction in randomised trials by inverse probability weighting using the propensity score.  Statistics in Medicine, 33(5):721–737, February 2014.
  * [20] ICH E9 working group.  Addendum on estimands and sensitivity analysis in clinical trials to the guideline on statistical principles for clinical trials, 2019.
  * [21] T. Permutt.  Do covariates change the estimand?  Statistics in Biopharmaceutical Research, 12(1):45–53, 2020.
  * [22] R. Daniel, J. Zhang, and D. Farewell.  Making apples from oranges: Comparing noncollapsible effect estimators and their standard errors after adjustment for different covariate sets.  Biometrical Journal, 63(3):528–557, 2021.
  * [23] V. Didelez and M. J. Stensrud.  On the logic of collapsibility for causal effect measures.  Biometrical Journal, feb 2021.
  * [24] A. Huitfeldt, M. J. Stensrud, and E. Suzuki.  On the collapsibility of measures of effect in the counterfactual causal framework.  Emerging Themes in Epidemiology, 16(1), jan 2019.
  * [25] J. Kim, A. B. Troxel, S. D. Halpern, K. G. Volpp, B. C. Kahan, T. P. Morris, and M. O. Harhay.  Analysis of multicenter clinical trials with very low event rates.  Trials, 21(1):917, 2020.
  * [26] T. P. Morris, B. C. Kahan, and I. R. White.  Choosing sensitivity analyses for randomised trials: principles.  BMC medical research methodology, 14(1):11+, 2014.
  * [27] C. Gamble, A. Krishan, D. Stocken, S. Lewis, E. Juszczak, C. Doré, P. R. Williamson, D. G. Altman, A. Montgomery, P. Lim, J. Berlin, S. Senn, S. Day, Y. Barbachano, and E. Loder.  Guidelines for the content of statistical analysis plans in clinical trials.  JAMA, 318(23):2337, 2017.
  * [28] L N Yelland, A B Salter, and P Ryan.  Relative risk estimation in randomized controlled trials: A comparison of methods for independent observations.  The International Journal of Biostatistics, 7(1):1–31, 2011.
  * [29] G. Heinze and M. Schemper.  A solution to the problem of separation in logistic regression.  Statistics in Medicine, 21:2409–2419, 2002.
  * [30] E. Colantuoni and M. Rosenblum.  Leveraging prognostic baseline variables to gain precision in randomized trials.  Statistics in Medicine, 34(18):2602–2617, August 2015.
  * [31] K. L. Moore and M. J. van der Laan.  Covariate adjustment in randomized trials with binary outcomes: Targeted maximum likelihood estimation.  Statistics in Medicine, 28(1):39–64, jan 2009.
  * [32] S. Zeng, F. Li, R. Wang, and F. Li.  Propensity score weighting for covariate adjustment in randomized clinical trials.  Statistics in Medicine, 2020.
  * [33] R. J. Desai and J. M. Franklin.  Alternative approaches for confounding adjustment in observational studies using weighting based on the propensity score: a primer for practitioners.  BMJ, page l5657, oct 2019.
  * [34] H. Raad, V. Cornelius, S. Chan, E. Williamson, and S. Cro.  An evaluation of inverse probability weighting using the propensity score for baseline covariate adjustment in smaller population randomised controlled trials with a continuous outcome.  BMC Medical Research Methodology, 20(1), 2020.
  * [35] T. R. Sullivan, I. R. White, A. B. Salter, P. Ryan, and K. J. Lee.  Should multiple imputation be the method of choice for handling missing data in randomized trials?  Statistical Methods in Medical Research, 27(9):2610–2626, dec 2016\.
  * [36] X. L. Meng.  Multiple-imputation inferences with uncongenial sources of input.  Statistical Science, 9:538–558, 1994.
  * [37] I. R. White and S. G. Thompson.  Adjusting for partially missing baseline measurements in randomised trials.  Statistics in Medicine, 24(7):993–1007, 2005.
  * [38] S. Cro, T. P. Morris, M. G. Kenward, and J. R. Carpenter.  Sensitivity analysis for clinical trials with missing continuous outcome data using controlled multiple imputation: A practical guide.  Statistics in Medicine, 2020.
  * [39] E. Wilson, C. Free, T. P. Morris, J. Syred, I. Ahamed, A. S. Menon-Johansson, M. J. Palmer, S. Barnard, E. Rezel, and P. Baraitser.  Internet-accessed sexually transmitted infection (e-STI) testing and results service: A randomised, single-blind, controlled trial.  PLOS Medicine, 14(12):e1002479, 2017.
  * [40] P. Royston and D. Altman.  External validation of a Cox prognostic model: principles and methods.  BMC Medical Research Methodology, 13(1):33+, 2013.
  * [41] I. R. White, R. Daniel, and P. Royston.  Avoiding bias due to perfect prediction in multiple imputation of incomplete categorical variables.  Computational Statistics & Data Analysis, 54(10):2267–2275, oct 2010.

## Figures

Figure 1: Upper panel: Data from four notional trials where individuals
recruited have different distributions of $X$. The two quadratic curves show
the data in the two arms. Lower panel: SE after no adjustment and after linear
adjustment for each of the three trials.

## Tables

Table 1: Some population-level summaries commonly used in clinical trials with
binary outcome measures

Outcome type | Summary measures | Collapsible?⋆
---|---|---
Continuous | Mean difference | Yes
Binary | Risk difference | Yes
| Risk ratio | Yes
| Odds ratio | No
Time-to-event | Hazard ratio | No
| Restricted mean survival | Yes
| time difference |

⋆Do conditional and marginal summary measures always coincide?

Table 2: An illustration of non-collapsibility of the odds ratio  | Stratum
---|---
| A | B | Both
Allocation | Dead | Alive | Dead | Alive | Dead | Alive
Intervention | 9 | 1 | 5 | 5 | 14 | 6
Control | 5 | 5 | 1 | 9 | 6 | 14
Odds ratio | 9 | 9 | 5.4
Table 3: Results of analyses of the GetTested trial. All models included main effects only. Link function is canonical unless otherwise specified. The dash symbol - means model did not converge, except for * where the log risk ratio estimated as 541, indicating separation for one or more covariates. Outcome | Summary | Adjustment | Model (variable modelled) | Treatment effect
---|---|---|---|---
measure | measure | method |  | estimate (SE)
Any test | Risk | Direct | Identity-link binomial (outcome) | -
(occurred | difference | Standardisation | Logistic (outcome) | 0.260 (0.021)
in 35%) |  | IPTW | Logistic (treatment) | 0.262 (0.021)
| Log risk | Direct | Poisson, robust SE (outcome) | -*
| ratio | Direct | Log-link binomial (outcome) | 0.797 (0.075)
|  | Standardisation | Logistic (outcome) | 0.796 (0.074)
|  | IPTW | Logistic (treatment) | 0.806 (0.075)
Any diagnosis | Risk | Direct | Identity-link binomial (outcome) | -
(occurred | difference | Standardisation | Logistic (outcome) | 0.015 (0.006)
in 1.6%) |  | IPTW | Logistic (treatment) | 0.013 (0.006)
| Log risk | Direct | Poisson, robust SE (outcome) | 0.972 (0.433)
| ratio | Direct | Log-link binomial (outcome) | 0.915 (0.412)
|  | Standardisation | Logistic (outcome) | 0.959 (0.412)
|  | IPTW | Logistic (treatment) | 0.855 (0.412)
Table 4: Points to consider on properties of each approach  Issue |  Direct |  Standardisation |  Inverse probability
---|---|---|---
|  adjustment |  |  weighting
Estimand for non-collapsible summary measures |  Conditional |  Marginal |  Marginal
For non-collapsible summary measures, true $\beta$ depends on… |  Covariates conditioned on in outcome model |  In-trial distribution of covariates |  In-trial distribution of covariates
Misspecification of covariate effects |  Loses efficiency vs. correctly specified model but expected to gain vs. no adjustment. True $\beta$ changes under non-collapsibility |  Loses efficiency vs. correctly specified model but expected to gain vs. no adjustment |  Loses efficiency vs. correctly specified model but expected to gain vs. no adjustment
Convergence |  Vulnerable |  Reasonable (but see GetTested experience) |  Solid
Stratification / minimisation handled by variance estimator |  Yes |  Yes |  Yes
Efficiency |  Asymptotically optimal |  Asymptotically optimal |  Asymptotically optimal
Standard error calculation |  Direct |  Delta method |  Robust, accounting for estimation of weights via joint estimating equations. Standard error can be biased downwards in small samples[34]
Handling of missing covariate data in order to target all-randomised population |  Missing indicator or single mean imputation (though neither is suitable with non-collapsible population summary measures) |  Missing indicator or single mean imputation |  Missing indicator or single mean imputation
Handling of missing outcome data in order to target all-randomised population |  Multiple imputation by-arm (or inverse probability of missingness weighting) |  Standardisation to all-randomised rather than complete-case sample; alternatively multiple imputation by-arm or inverse probability of missingness weighting |  Inverse probability of missingness weighting (or multiple imputation by-arm )

Table 5: a) Full data from a notional randomised trial. b) Complete cases,
where all individuals with $X=0$ and half of individuals with X=1 are complete
cases. c) True value of summary measure within levels of $X$.

a) Full data
---
| $X=0$ | $X=1$
| $Z=0$ | $Z=1$ | $Z=0$ | $Z=1$
$Y=1$ | 42 | 26 | 180 | 140
$Y=0$ | 458 | 474 | 320 | 360
| 500 | 500 | 500 | 500

b) Complete cases
---
| $X=0$ | $X=1$
| $Z=0$ | $Z=1$ | $Z=0$ | $Z=1$
$Y=1$ | 21 | 13 | 180 | 140
$Y=0$ | 229 | 237 | 320 | 360
| 250 | 250 | 500 | 500

c) Value of summary measures within $X$ (both
---
in all randomised and among complete cases)
Summary measure | $X=0$ | $X=1$
Odds ratio | $0.598$ | $0.691$
Risk ratio | $0.619$ | $0.778$
Risk difference | $-0.032$ | $-0.080$

Table 6: Summary of the true value of the estimand in all randomised and among
the complete cases under covariate-dependent missingness.

Summary measure | All-randomised | Complete cases
---|---|---
Conditional odds ratio | $0.670$ | $0.679$
Marginal odds ratio | $0.698$ | $0.700$
Risk ratio | $0.748$ | $0.761$
Risk difference | $-0.056$ | $-0.064$

*This method is possible with complete $X$ and incomplete $Y$

Table 7: Appendix 2: Results of analyses of the GetTested trial targeting the
all-randomised population. All models included main effects only. Link
function is canonical unless otherwise specified. The dash symbol - means
model did not converge, with reasons described in the text of appendix 2.

Outcome | Summary | Adjustment | Model (variable modelled) | Treatment effect
---|---|---|---|---
measure | measure | method |  | estimate (SE)
Any test | Risk | Direct | Identity-link binomial (outcome) | -
(occurred | difference | Standardisation | Logistic (outcome) | 0.258 (0.021)*
in 35%) |  | IPTW | Logistic (treatment) | 0.259 (0.021)
| Log risk | Direct | Poisson, robust SE (outcome) | 0.804 (0.069)
| ratio | Direct | Log-link binomial (outcome) | -
|  | Standardisation | Logistic (outcome) | 0.799 (0.075)*
|  | IPTW | Logistic (treatment) | 0.802 (0.075)
Any diagnosis | Risk | Direct | Identity-link binomial (outcome) | -
(occurred | difference | Standardisation | Logistic (outcome) | -
in 1.6%) |  | IPTW | Logistic (treatment) | 0.013 (0.006)
| Log risk | Direct | Poisson, robust SE (outcome) | -
| ratio | Direct | Identity-link binomial (outcome) | -
|  | Standardisation | Logistic (outcome) | -
|  | IPTW | Logistic (treatment) | 0.866 (0.414)

*Almost-all-randomised. Estimate was returned only after omitting four participants affected by collinearity.

## Additional Files

### Additional file 1 – Stata code to generate figure 1

### Additional file 2 – Stata code for the analysis of the GetTested trial

(Assumes the data file journal.pmed.1002479.s001.xls has been downloaded from
https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.1002479#sec020)


# Optical spectroscopy of the $5p_{3/2}\to 6p_{1/2}$ electric dipole forbidden
transition in atomic rubidium.

F. Ponciano-Ojeda  C. Mojica-Casique  S. Hernández-Gómez  O. López-Hernández
J. Flores-Mijangos  F. Ramírez-Martínez<EMAIL_ADDRESS>Instituto de
Ciencias Nucleares, UNAM. Circuito Exterior, Ciudad Universitaria, 04510
México City, México.  D. Sahagún  R. Jáuregui  Instituto de Física, UNAM,
Ciudad Universitaria, 04510 México City, México.  J. Jiménez-Mier
<EMAIL_ADDRESS>Instituto de Ciencias Nucleares, UNAM. Circuito
Exterior, Ciudad Universitaria, 04510 México City, México.

###### Abstract

We present the first evidence of excitation of the $5p_{3/2}\rightarrow
6p_{1/2}$ electric dipole-forbidden transition in atomic rubidium. The
experiments were carried out in a rubidium vapor cell using Doppler-free
optical-optical double-resonance spectroscopy with counter-propagating beams.
A $5s_{1/2}\rightarrow 5p_{3/2}$ electric dipole preparation step using a
diode laser locked to the $F=3\to 4$ cyclic transition of the $D2$ line in
85Rb is used to prepare the atoms in the first excited state. This is then
followed by the $5p_{3/2}\,F_{2}=4\rightarrow 6p_{1/2}\,F_{3}$ dipole-
forbidden excitation ($\lambda\approx 917.5$ nm) to establish a two-photon
ladder ($\Xi$) excitation scheme. Production of atoms in the $6p_{1/2}$
excited state is verified by detection of the $421$ nm fluorescence that
results from direct decay into the $5s_{1/2}$ ground state. The polarization
dependence of the relative intensities of the lines of the decay fluorescence
is also investigated. Experimental data for different polarization
configurations of the light beams used in this two-photon spectroscopy are
compared with the results of calculations that consider a strong atom-field
coupling in the preparation step, followed by a weak electric quadrupole
excitation and the blue fluorescence decay emission. Good agreement between
experiment and this three-step model is found in the case of linear-linear
polarizations.

###### pacs:

32.70.Cs,32.70.Fw

##  I Introduction

The study of transitions beyond the electric dipole approximation, commonly
referred to as “forbidden transitions”, in the interaction between atoms and
optical radiation fields has played a significant role in the in-depth study
of atomic structures. Once limited to observation in high energy contexts,
such as the studies of astrophysical phenomena and plasmas Biémont and Zeippen
(1996), these forbidden transitions are now frequently used in other areas
like metrology van Rooij  _et al._ (2011), three-wave mixing experiments
Flusberg  _et al._ (1977a, b) and parity conservation experiments Bouchiat and
Bouchiat (1997). Thanks in part to the advent of high powered continuous-wave
(CW) and pulsed laser sources, weak absorption lines, such as those given by
forbidden transitions, have been observed and studied in alkali-metal vapors
Guéna  _et al._ (1987); Tojo  _et al._ (2005); Bayram  _et al._ (2000) and
cold atoms Bhattacharya  _et al._ (2003); Pires  _et al._ (2009); Tong  _et
al._ (2009), as well as recently with optical-optical double resonance
spectroscopy Ponciano-Ojeda  _et al._ (2015); Chan  _et al._ (2016); Mojica-
Casique  _et al._ (2016).

In particular, $p\to p$ transitions in alkali atoms are an important part of
the study of forbidden transitions. Early experiments found that the
$np\to(n+1)p$ excitation, with $n$ the ground state principal quantum number,
is strong enough to be produced with moderate CW laser power Bhattacharya  _et
al._ (2003). Also, experiments with cold atoms conclusively proved that the
$5p\to 8p$ forbidden transition in rubidium is an electric quadrupole
transition with no significant contribution from the magnetic dipole term
Pires  _et al._ (2009). Our research group recently showed Ponciano-Ojeda  _et
al._ (2015); Mojica-Casique  _et al._ (2016) that the hyperfine structure of
the $6p_{3/2}$ state could be resolved in the $5p_{3/2}\to 6p_{3/2}$ forbidden
spectra of room temperature rubidium atoms. This is possible because it is
necessary first to prepare atoms in the $5p_{3/2}$ first excited state, so
these are naturally collinear optical-optical double resonance experiments
which allow the study of atomic transitions free of Doppler broadening
Demtröder (2015). In this experimental scheme one of the laser beams prepares
atomic populations in a first excited state while the second maps the
resulting population distribution onto a second excited state. Fixing the
frequency of the first of the aforementioned lasers also allows velocity-
selective spectra Pinard  _et al._ (1979) to be obtained based on the
direction of propagation and the wavelength difference of the second beam
relative to the fixed one. Furthermore, we also showed that this forbidden
transition is very sensitive to the polarization states of both preparation
and excitation lasers, and that the choice of polarization configurations
allows a direct test of electric quadrupole selection rules over the atomic
magnetic quantum numbers Mojica-Casique  _et al._ (2016).

In this paper we extend the study of forbidden transitions in room temperature
rubidium atoms to the $5p_{3/2}\to 6p_{1/2}$ electric quadrupole (E2)
excitation. This is a $p\to p$ transition to a different fine structure state
($J=1/2$) of the $6p$ manifold, which has a different hyperfine structure and
also different transition matrix elements. We make a detailed analysis of the
dependence of the E2 spectra on the relative polarization directions of the
$5s_{1/2}\to 5p_{3/2}$ linearly polarized preparation beam and the linearly
polarized beam responsible of the E2 transition. The results are then
interpreted by applying the model developed for the $5p_{3/2}\to 6p_{3/2}$
excitation Mojica-Casique  _et al._ (2016). This model assumes that a strong
preparation step establishes the relative populations of the
$5p_{3/2}\,FM_{F}$ magnetic sublevels and that the E2 transition intensity to
the different $6p_{1/2}$ hyperfine states is determined by the geometric part
of the electric quadrupole transition matrix element.

##  II Experimental setup

An energy level diagram showing the total angular momentum quantum numbers and
hyperfine splittings for the $6p_{1/2}$ state in 85Rb is shown in figure 1. In
our experiment the atoms are first prepared in the $5p_{3/2}$ excited state by
an electric dipole transition, where the populations in the $M_{F^{\prime}}$
magnetic sublevels reach a stationary state. This is the same preparation step
as used in Ponciano-Ojeda  _et al._ (2015); Mojica-Casique  _et al._ (2016). A
second laser, whose frequency is swept across several hundred MHz, pumps atoms
in the $5p_{3/2}$ excited state to the hyperfine $F^{\prime\prime}$ levels of
the $6p_{1/2}$ state. The hyperfine structure of the $6p_{1/2}$ state, as well
as secondary cross-over features due to velocity-selective resonances, are
resolved in spectra obtained by means of this excitation scheme. This
spectroscopy, when realized with different polarization configurations of the
laser beams, demonstrate the control on the preparation of the $6p_{1/2}$
hyperfine states via the E2 selection rules Mojica-Casique  _et al._ (2016).

Figure 1: Energy levels of the $6p_{1/2}$ state (85Rb). Hyperfine splittings
are shown in MHz. The electric dipole excitation (D2 line) is locked to the
$F\rightarrow F+1$ cyclic transition. The frequency of the electric quadrupole
(E2) excitation laser is swept across the $6p_{1/2}$ hyperfine manifold in
order to resolve the structure.

The setup used for this experiment is similar to that reported in Ponciano-
Ojeda  _et al._ (2015), with the main difference being the use of a commercial
titanium-saphire laser MSq  for the $917.5$ nm electric quadrupole excitation
light. A homebuilt external cavity diode laser (ECDL)Arnold  _et al._ (1998);
Hawthorn  _et al._ (2001) in resonance with the $5s_{1/2}\rightarrow 5p_{3/2}$
transition at $780$ nm (D2 line) is used to prepare atoms in the $5p_{3/2}$
hyperfine states. Spontaneous decay of the excited atoms from the $6p_{1/2}$
hyperfine levels to the $5s_{1/2}$ ground state produces fluorescence at $421$
nm that is used for detection of the $5p_{3/2}\rightarrow 6p_{1/2}$
transition. Photons produced from this decay inside the cell are collected by
a system of lenses that focuses them onto the cathode of a photomultiplier
tube (PMT), and the output is later processed with a lock-in amplifier to
enhance the signal-noise ratio. One expects from the theory that a doublet
with the frequency splitting of the well known hyperfine structure of the
$6p_{1/2}$ state Arimondo  _et al._ (1977); Feiertag and zu Putlitz (1973);
Sansonetti (2006) will appear in the obtained spectra. Both lasers are
linearly polarized. A half-wave plate in the path of the $917.5$ nm laser is
used to vary the angle between the polarization directions.

##  III Calculation of relative line intensities.

The overall geometry used in calculation of the relative intensities of the
fluorescence resonances is shown in figure 2. The $780$ and $917.5$ beams
(preparation and non-dipole excitation, respectively) propagate along the $x$
axis in opposite directions with linear polarizations along $\mathbf{E_{1}}$
and $\mathbf{E_{2}}$. The linear polarization of the preparation laser is kept
fixed along the $z$ axis. The polarization of the $917.5$ beam is rotated so
that an angle $\theta$ is formed between $\mathbf{E_{1}}$ and
$\mathbf{E_{2}}$. The detection of the fluorescence is performed along the
positive direction of the $y$ axis by the photomultiplier tube.

The relative intensities of the observed fluorescence produced in this two-
photon excitation can be calculated using the simple three-step model
developed for the $5p_{3/2}\to 6p_{3/2}$ forbidden transition in references
Ponciano-Ojeda  _et al._ (2015) and Mojica-Casique  _et al._ (2016). The
probability to observe a $421$ nm photon from the decay of the
$|6p_{1/2}F_{3}\rangle$ hyperfine state is given by the expression Ponciano-
Ojeda  _et al._ (2015); Mojica-Casique  _et al._ (2016)

$\displaystyle
P(F_{3})=\sum_{M_{2},M_{3},F^{\prime}_{1},M^{\prime}_{1},\lambda}$
$\displaystyle\sigma(F_{2},M_{2})|\langle
5p_{3/2}F_{2}M_{2}|T|6p_{1/2}F_{3}M_{3}\rangle|^{2}$
$\displaystyle\times|\langle
6p_{1/2}F_{3}M_{3}|D_{\lambda}|5s_{1/2}F^{\prime}_{1}M^{\prime}_{1}\rangle|^{2}$
(1)

where $\sigma(F_{2},M_{2})$ is the population of the $5p_{3/2}\,F_{2}M_{2}$
states produced by the strong, electric dipole $5s_{1/2}\to 5p_{3/2}$
preparation step, $\langle
5p_{3/2}F_{2}M_{2}|T|6p_{1/2}F_{3}M_{3}\rangle|^{2}$ is the weak, electric
quadrupole transition probability, and $|\langle
6p_{1/2}F_{3}M_{3}|D_{\lambda}|5s_{1/2}F^{\prime}_{1}M^{\prime}_{1}\rangle|^{2}$
is the probability of decay from the $F_{3}M_{3}$ magnetic state of the
$6p_{1/2}$ manifold. It is important to notice that the preparation step in
this experiment is the one used to study the $5p_{3/2}\to 6p_{3/2}$ transition
and thus the values of the $\sigma(F_{2},M_{2})$ populations used in this work
are the same as the ones used previously Ponciano-Ojeda  _et al._ (2015);
Mojica-Casique  _et al._ (2016). Hyperfine and Zeeman pumping effects in this
preparation step establish the magnetic state populations of the $5p_{3/2}$
hyperfine manifold which then determine the response of the atoms to the non-
dipole transition. In the transition studied here the only changes necessary
are to consider the $J=1/2$ level in the $6p$ fine structure and its
respective hyperfine states in the electric quadrupole excitation and the
$6p_{1/2}\to 5s_{1/2}$ electric dipole decay transition matrix elements. All
the approximations made in the case of the $5p_{3/2}\rightarrow 6p_{3/2}$
transition can also be seen to hold for the present case.

Figure 2: (Color online) Geometry used for the calculation of relative
intensities of the $421$ fluorescence using the three step model presented in
Ponciano-Ojeda  _et al._ (2015); Mojica-Casique  _et al._ (2016).
$\mathbf{E}_{1}$ is the polarization direction of the $780$ nm preparation
laser; $\mathbf{E}_{2}$ is the polarization direction of the $917.5$ nm laser
that excites the quadrupole transition. Both beams propagate along the $x$
axis and the $421$ nm fluorescence is detected along the $y$ axis. See text
for further details.

It is also important to recall that different polarization configurations give
place to different selections rules for the electric dipole forbidden
quadrupole step. In ref. Mojica-Casique  _et al._ (2016) it was shown that for
parallel linear polarizations we have the selection rule $\Delta M=\pm 2$
while for the perpendicular case we have $\Delta M=\pm 1$.For intermediate
values of $\theta$ this model predicts spectra that are superposition of these
two extremes, with weighing factors proportional to $\cos^{2}\theta$ for the
parallel spectrum and $\sin^{2}\theta$ for the perpendicular spectrum.
Therefore the angular dependence can be written Mojica-Casique  _et al._
(2016) in terms of the second order Legendre polynomial of $\cos\theta$. The
fluorescence signal recorded at the ”magic angle”
($\theta_{M}=\cos^{-1}\sqrt{1/3}\approx 54.7^{\circ}$) is given by
$I(\theta_{M})=I_{\parallel}+2I_{\perp}$ which is a natural normalization
factor for a direct comparison between experiment and theory.

##  IV Results. Comparison between experiment and theory.

Typical spectra showing the $421$ nm decay fluorescence as the frequency of
the $917.5$ nm laser is scanned are shown in figure 3. Here we only show
spectra for parallel and perpendicular polarizations, but data were obtained
for an entire angular distribution study. The frequency scale of these plots
was obtained by using the calculated hyperfine splitting for the $6p_{1/2}$
two main peaks to adjust the distance between peak centers Feiertag and zu
Putlitz (1973); Sansonetti (2006). The frequency zero-reference was then
shifted to the center of gravity of the $6p_{1/2}$ hyperfine manifold of the
isotope in question. A least squares fit of line profiles was done for the
main observable peaks in each spectrum. The peaks were all adjusted with the
same value for the width whilst the centers and heights varied from spectrum
to spectrum. Finally, the intensity of the $F=2$, $3$ spectra recorded for
parallel and perpendicular polarizations were used to calculate the
normalization factor

$\sum_{F}I(\theta_{M})=\sum_{F}(I_{\parallel}+2I_{\perp})$ (2)

Figure 3: Spectra showing the $421$ nm decay fluorescence from the $6p_{1/2}$
state in 85Rb for parallel polarizations (bottom) and perpendicular
polarizations (top). The center for each of the peaks, corresponding to the
calculated position for each hyperfine level $F$, is indicated with vertical
lines. Voigt fits for each of the main peaks are indicated by the thick (red;
color online) lines, with a common width of $7.6\pm 0.1$ MHz. The expected
position of the velocity-selected peaks is also indicated.

For the two spectra one clearly observes the two expected lines that result
from an excitation sequence $5s_{1/2},F=3\rightarrow 5p_{3/2},F=4\rightarrow
6p_{1/2},F_{3}$, ($F_{3}=2$ and $3$) for zero velocity atoms. Due to the
experiments being carried out in a room-temperature vapor cell there are also
groups of atoms with non-zero velocity projections that are excited by the
preparation beam. This results in secondary peaks, which are due to the
partial compensation of the Doppler shift of the preparation beam by the
$917.5$ nm counter-propagating beam, and give place to the dipole-forbidden
transitions appearing at a different frequency to those obtained with the
maximum $F_{2}$ preparation Ponciano-Ojeda  _et al._ (2018). The strongest of
these velocity-selected non-dipole transitions results from the $F\rightarrow
F\rightarrow F$ excitation chain ($2\rightarrow 2\rightarrow 2$ in 87Rb and
$3\rightarrow 3\rightarrow 3$ in 85Rb). In the 85Rb spectra small shoulders
are observed at $\approx 19$ MHz above the $F=2$, $3$ peaks, in good agreement
with the position of the velocity-selected transition expected to appear at
$18.1$ MHz above the zero velocity excitation. One can also observe an
additional small feature above the first velocity-selected transition for the
$F=3$ peak that corresponds to the $F\rightarrow F-1\rightarrow F$ excitation
chain, expected at $27.6$ MHz above the zero velocity excitation. The
corresponding feature for the $F=2$ peak is not visible due to the noise.

The fits carried out on the spectra also give some information about the
relative intensities of the hyperfine lines. Using the expressions for the
angular dependence of the transition probabilities given in Mojica-Casique
_et al._ (2016) one can calculate approximate values for the relative
intensities. These values are compared to the experimental ones in table 1.

| Parallel | Perpendicular
---|---|---
| Experimental | Calculated | Experimental | Calculated
$F=2$ | 35.8 | 39.1 | 8.9 | 5.9
$F=3$ | 19.8 | 23.1 | 13.4 | 12.9
Table 1: Calculated and experimental percentages for the relative intensities
of the hyperfine peaks in spectra for 85Rb as a function of the polarization
configuration of the excitation beams. We have taken the normalization factor
given by eq. 2 (see text).

The result of the measurements of the angular dependence of the fluorescence
line intensities as functions of the angle $\theta$ between polarization
directions is shown in figure 4. The figure also shows the results of the
angular distribution calculated using the three step model Mojica-Casique  _et
al._ (2016). For this comparison the experimental data were normalized to the
total intensity of peaks $F=2$ and $3$ (eq. 2). Also, a common background was
subtracted from all spectra. As for the $5p_{3/2}\to 6p_{3/2}$ electric
quadrupole transition Mojica-Casique  _et al._ (2016), there is very good
agreement between experiment and theory. The maxima in Fig. 4 for the electric
quadrupole transition into both hyperfine states occur for parallel
polarizations, with minima for perpendicular polarizations. Both experiment
and theory give a much stronger angular dependence of the transition to the
$F=2$ state. These results confirm the validity of the theoretical model, and
provide an indication of the ability to use the polarization of the light to
control the population of hyperfine states excited by electric quadrupole
transitions.

Figure 4: Angular dependence of relative intensities of the $F=3$ and $F=2$
hyperfine lines in 85Rb. The experimental data are compared with the results
of the three step model calculation using the expressions obtained in Mojica-
Casique  _et al._ (2016).

##  V Conclusions

We have presented a new $p\to p$ electric dipole forbidden transition in
thermal rubidium atoms, namely, the $5p_{3/2}\to 6p_{1/2}$ electric quadrupole
transition. It was also demonstrated that the three step model is a very good
approximation of the excitation dynamics and plays a key role in the correct
interpretation of the experimental results. Both experimental data and the
model show that the relative intensities of the forbidden spectra strongly
depend on the angle $\theta$ between the linear polarization directions of
preparation and excitation beams. The model shows that this is mainly the
result of different population distributions among the magnetic sublevels of
both $5p_{3/2}$ and $6p_{1/2}$ excited states. The results obtained serve as
complementary information to the $5p_{3/2}\rightarrow 6p_{3/2}$ transition
reported in refs. Ponciano-Ojeda  _et al._ (2015); Mojica-Casique  _et al._
(2016), thus giving a more complete picture of the $5p_{3/2}\rightarrow
6p_{J}$ electric dipole-forbidden transitions in atomic rubidium.

###### Acknowledgements.

We thank J. Rangel for his help in the construction of the diode laser. This
work was supported by DGAPA-UNAM, México, under projects PAPIIT Nos. IN116309,
IN110812, and IA101012, by CONACyT, México, under Basic Research project No.
44986 and National Laboratory project LN260704.

## References

  * Biémont and Zeippen (1996) E. Biémont and C. J. Zeippen, Phys. Scr. T65, 192 (1996).
  * van Rooij  _et al._ (2011) R. van Rooij, J. S. Borbely, J. Simonet, M. D. Hoogerland, K. S. E. Eikema, R. A. Rozendaal, and W. Vassen, Science 333, 196 (2011).
  * Flusberg  _et al._ (1977a) A. Flusberg, T. Mossberg, and S. R. Hartmann, Phys. Rev. Lett. 38, 59 (1977a).
  * Flusberg  _et al._ (1977b) A. Flusberg, T. Mossberg, and S. R. Hartmann, Phys. Rev. Lett. 38, 694 (1977b).
  * Bouchiat and Bouchiat (1997) M.-A. Bouchiat and C. Bouchiat, Rep. Prog. Phys. 60, 1351 (1997).
  * Guéna  _et al._ (1987) J. Guéna, M. Lintz, P. Jacquier, L. Pottier, and M. Bouchiat, Opt. Comm. 62, 97 (1987).
  * Tojo  _et al._ (2005) S. Tojo, T. Fujimoto, and M. Hasuo, Phys. Rev. A 71, 012507 (2005).
  * Bayram  _et al._ (2000) S. B. Bayram, M. D. Havey, D. V. Kupriyanov, and I. M. Sokolov, Phys. Rev. A 62, 012503 (2000).
  * Bhattacharya  _et al._ (2003) M. Bhattacharya, C. Haimberger, and N. P. Bigelow, Phys. Rev. Lett. 91, 213004 (2003).
  * Pires  _et al._ (2009) R. Pires, M. Ascoli, E. E. Eyler, P. L. Gould, and A. Derevianko, Phys. Rev. A 80, 062502 (2009).
  * Tong  _et al._ (2009) D. Tong, S. M. Farooqi, E. G. M. van Kempen, Z. Pavlovic, J. Stanojevic, R. Côté, E. E. Eyler, and P. L. Gould, Phys. Rev. A 79, 052509 (2009).
  * Ponciano-Ojeda  _et al._ (2015) F. Ponciano-Ojeda, S. Hernández-Gómez, O. López-Hernández, C. Mojica-Casique, R. Colín-Rodríguez, F. Ramírez-Martínez, J. Flores-Mijangos, D. Sahagún, R. Jáuregui, and J. Jiménez-Mier, Phys. Rev. A 92 (2015), 10.1103/PhysRevA.92.042511.
  * Chan  _et al._ (2016) E. A. Chan, S. A. Aljunid, N. I. Zheludev, D. Wilkowski, and M. Ducloy, Opt. Lett. 41, 2005 (2016).
  * Mojica-Casique  _et al._ (2016) C. Mojica-Casique, F. Ponciano-Ojeda, S. Hernández-Gómez, O. López-Hernández, J. Flores-Mijangos, F. Ramírez-Martínez, D. Sahagún, R. Jáuregui, and J. Jiménez-Mier, J. Phys. B 50, 025003 (2016).
  * Demtröder (2015) W. Demtröder,  _Laser Spectroscopy Vol. 2. Experimental Techniques._ , 5th ed. (Springer-Verlag, 2015).
  * Pinard  _et al._ (1979) M. Pinard, C. G. Aminoff, and F. Laloë, Phys. Rev. A 19, 2366 (1979).
  * (17) SolsTiS-1600-SRX-F laser system, MSquared Lasers Inc.
  * Arnold  _et al._ (1998) A. S. Arnold, J. S. Wilson, and M. G. Boshier, Rev. Sci. Instrum. 69, 1236 (1998).
  * Hawthorn  _et al._ (2001) C. J. Hawthorn, K. P. Weber, and R. E. Scholten, Rev. Sci. Instrum. 72, 4477 (2001).
  * Arimondo  _et al._ (1977) E. Arimondo, M. Inguscio, and P. Violino, Rev. Mod. Phys. 49, 31 (1977).
  * Feiertag and zu Putlitz (1973) D. Feiertag and G. zu Putlitz, Z. Physik 261, 1 (1973).
  * Sansonetti (2006) J. E. Sansonetti, J. Phys. Chem. Ref. Data 35, 301 (2006).
  * Ponciano-Ojeda  _et al._ (2018) F. Ponciano-Ojeda, S. Hernández-Gómez, C. Mojica-Casique, E. Ruiz-Martínez, O. López-Hernández, R. Colín-Rodríguez, F. Ramírez-Martínez, J. Flores-Mijangos, D. Sahagún, R. Jáuregui, and J. Jiménez-Mier, Am. J. Phys 86, 6 (2018).