arXiv:1705.00951v1 [stat.ME] 2 May 2017

A mean score method for sensitivity analysis
to departures from the missing at random
assumption in randomised trials

Ian R. White1,2,∗ , James Carpenter
1

2

3

2,3

, and Nicholas J. Horton

MRC Biostatistics Unit, Cambridge, UK

MRC Clinical Trials Unit at UCL, London, UK

London School of Hygiene and Tropical Medicine, UK
4

Amherst College, Amherst, MA, USA
∗

Email: ian.white@ucl.ac.uk
May 3, 2017

4

Abstract

Most analyses of randomised trials with incomplete outcomes make untestable
assumptions and should therefore be subjected to sensitivity analyses. However, methods for sensitivity analyses are not widely used. We propose a
mean score approach for exploring global sensitivity to departures from missing at random or other assumptions about incomplete outcome data in a
randomised trial. We assume a single outcome analysed under a generalised
linear model. One or more sensitivity parameters, specified by the user, measure the degree of departure from missing at random in a pattern mixture
model. Advantages of our method are that its sensitivity parameters are relatively easy to interpret and so can be elicited from subject matter experts;
it is fast and non-stochastic; and its point estimate, standard error and confidence interval agree perfectly with standard methods when particular values
of the sensitivity parameters make those standard methods appropriate. We
illustrate the method using data from a mental health trial.
Keywords: Intention-to-treat analysis, Longitudinal data analysis, Mean
score, Missing data, Randomised trials, Sensitivity analysis.
Running title: Mean score method for sensitivity analysis.

1

Introduction

Missing outcome data are a threat to the validity of randomised controlled
trials, and they usually require untestable assumptions to be made in the
analysis. One common assumption is that data are missing at random (MAR)
(Little and Rubin, 2002). Other possible assumptions may be less implausible in particular clinical settings. For example, in smoking cessation trials,
the outcome is binary, indicating whether an individual quit over a given
period, and it is common to assume that missing values are failures — “missing=failure” (West et al., 2005); while in weight loss trials, missing data is
sometimes assumed to be unchanged since baseline — “baseline observation
carried forward” (Ware, 2003).
The US National Research Council (2010) suggested measures that should
be taken to minimise the amount of missing outcome data in randomised trials, and described analysis strategies based on various assumptions about the
missing data. This report recommended that “Sensitivity analyses should be
part of the primary reporting of findings from clinical trials. Examining
sensitivity to the assumptions about the missing data mechanism should be
a mandatory component of reporting.” However, among “several important areas where progress is particularly needed”, the first was “methods
for sensitivity analysis and principled decision making based on the results

1

from sensitivity analyses”. Sensitivity analysis is also an essential part of an
intention-to-treat analysis strategy, which includes all randomised individuals in the analysis strategy (White et al., 2011a, 2012): even if the main
analysis is performed under MAR and hence draws no information from individuals with no outcome data, such individuals are included in sensitivity
analysis and hence in the analysis strategy.
Sensitivity analysis is often done by performing two different analyses,
such as an analysis assuming MAR and an analysis by last observation carried forward, and concluding that inference is robust if the results are similar
(Wood et al., 2004). Better is a principled sensitivity analysis, where the
data analyst typically formulates a model including unidentified ‘sensitivity
parameter(s)’ that govern the degree of departure from the main assumption
(e.g. from MAR), and explores how the estimate of interest varies as the
sensitivity parameter(s) are varied (Rotnitzky et al., 1998; Kenward et al.,
2001). We consider global sensitivity analyses where the sensitivity parameter(s) are varied over a range of numerical values that subject-matter experts
consider plausible.
Likelihood-based analyses assuming MAR can usually ignore the missing
data mechanism and simply analyse the observed data (Little and Rubin,
2002). Under a missing not at random (MNAR) assumption, however, it is
usually necessary to model the data of interest jointly with the assumed miss2

ing data mechanism. The joint model can be specified as a pattern-mixture
model, which explicitly describes the differences between profiles of patients
who complete and drop out (Little, 1993, 1994), or as a selection model,
which relates the chance of drop-out to the (possibly missing) response values either directly (Diggle and Kenward, 1994; Kenward, 1998) or indirectly
through a random effect (Follmann and Wu, 1995; Roy, 2003). Rotnitzky
et al. (1998) proposed a selection model for incomplete repeated measures
data and showed how to estimate it by inverse probability weighting, given
values of the sensitivity parameters. Scharfstein et al. (2003) adopted a
non-parametric Bayesian approach to analysing incomplete randomised trial
data, and argued that sensitivity parameters are more plausibly a priori independent of other parameters of interest in a selection model than in a
pattern-mixture model. Scharfstein et al. (2014) proposed a fully parametric approach based on a selection model. On the other hand, Daniels and
Hogan (2000) advocated a pattern-mixture framework as “a convenient and
intuitive framework for conducting sensitivity analyses”.
We use the pattern-mixture model in this paper because its sensitivity
parameters are usually more easily interpreted (White et al., 2007, 2008).
For a binary outcome, a convenient sensitivity parameter is the informative
missing odds ratio (IMOR), defined, conditional on covariates, as the odds of
positive outcome in missing values divided by the odds of positive outcome in
3

observed values (Higgins et al., 2008; Kaciroti et al., 2009). For a continuous
outcome, a convenient sensitivity parameter is the covariate-adjusted mean
difference between missing and observed outcomes (Mavridis et al., 2015).
Most estimation procedures described for general pattern-mixture models are likelihood-based (Little, 1993, 1994; Little and Yau, 1996; Hedeker
and Gibbons, 2006), while the National Research Council (2010) describes
point estimation using sensitivity parameters with bootstrap standard errors. In this paper we propose instead using the mean score method, a
computationally convenient method which was originally proposed under a
MAR assumption for incomplete outcome data (Pepe et al., 1994) and for
incomplete covariates (Reilly and Pepe, 1995). The method is particularly
useful to allow for auxiliary variables, so that outcomes can be assumed MAR
given model covariates and auxiliary variables but not necessarily MAR given
model covariates alone (Pepe et al., 1994). We are not aware of the mean
score method having been used for sensitivity analysis.
The aim of this paper is to propose methods for principled sensitivity
analysis that are fast, non-stochastic, available in statistical software, and
agree exactly with standard methods in the special cases where standard
methods are appropriate. We focus on randomised trials with outcome measured at a single time, allowing for continuous or binary outcomes, or indeed
any generalised linear model, and for covariate adjustment.
4

The paper is organised as follows. Section 2 describes our proposed
method. Section 3 proposes small-sample corrections which yield exact equivalence to standard procedures in special cases. Section 4 illustrates our
method in QUATRO, a mental health trial with outcome measured at a
single time. Section 5 describes a simulation study. Section 6 discusses the
implementation of our method, possible alternatives, limitations and extensions.

2

Mean score approach

Assume that for the ith individual (i = 1 to n) in an individually randomised
trial, there is an outcome variable yi , and let ri be an indicator of yi being
observed. Let nobs and nmis = n − nobs be the numbers of observed and
missing values of y respectively. Let xi be a vector of covariates including
the pS -dimensional fully-observed covariates xSi in the substantive model,
comprising an intercept, an indicator zi for the randomised group, and (optionally) baseline covariates. xi may also include fully-observed auxiliary
covariates xAi that are not in the substantive model but that help to predict
yi , and/or covariates xRi that are only observed in individuals with missing
yi and describe the nature of the missing data: for example, the reason for
missingness.

5

The aim of the analysis is to estimate the effect of randomised group,
adjusting for the baseline covariates. We assume the substantive model is a
generalised linear model (GLM) with canonical link,
E [yi |xSi ; β S ] = h(β TS xSi )

(1)

where h(.) is the inverse link function. We are interested in estimating βSz ,
the component of the pS -dimensional vector β S corresponding to z.
If we had complete data, we would estimate β S by solving the estimating
equation US∗ (β S ) = 0 where US∗ (β S ) =

P

i

∗
(β S ) and
USi

∗
USi
(β S ) = {yi − h(β TS xSi )}xSi .

(2)

The mean score approach (Pepe et al., 1994; Reilly and Pepe, 1995) handles missing data by replacing US∗ (β S ) with US (β S ), its expectation over
the distribution of the missing data given the observed data. We write
US (β S ) =

P

i

∗
USi (β S ) and USi (β S ) = E [USi
(β S )|xi , ri , ri yi ]. Then by the re-

∗
(β S )|xSi ] since xi includes
peated expectation rule, E [USi (β S )|xSi ] = E [USi

xSi , so US (β S ) = 0 is an unbiased estimating equation if US∗ (β S ) = 0 is.
To compute USi (β S ), we need only E [yi |xi , ri = 0], because (2) is linear
in yi . We estimate this using the pattern-mixture model




E [yi |xi , ri ; β P ] = h β TP xP i + ∆(xi )(1 − ri )

(3)

where xP i = (xSi , xAi ) of dimension pP ; the subscript P distinguishes the
6

parameters β P and covariates xP i of the pattern-mixture model from the
parameters β S and covariates xSi of the substantive model. Models (1) and
(3) are typically not both correctly specified: we return to this issue in the
simulation study.
In equation (3), ∆(xi ) is a user-specified departure from MAR for individual i. MAR in this setting means that p(ri = 1|xi , yi ) = p(ri = 1|xi ),
which implies E [yi |xi , ri ] = E [yi |xi ] and hence ∆(xi ) = 0 for all i. A simple
choice of ∆(xi ) that expresses MNAR is ∆(xi ) = δ, where the departure
from MAR is the same for all individuals. Differences in departure from
MAR between randomised groups are often plausible and can have strong
impact on the estimated treatment effect (White et al., 2007), so an alternative choice is ∆(xi ) = δzi . The departure ∆(xi ) could also depend on reasons
for missingness coded in xRi : for example, it could be 0 for individuals lost
to follow-up (if MAR seemed plausible for them), and δzi for individuals who
refused follow-up.
Putting it all together, the mean score method solves
{ỹi (β P ) − h(β TS xSi )}xSi = 0

X

(4)

i





where ỹi (β P ) is defined as yi if ri = 1 and h β TP xP i + ∆(xi )(1 − ri ) if ri = 0.

7

2.1

Estimation using full sandwich variance

The parameter β P in (3) is estimated by regressing yi on xP i in the complete
cases (ri = 1). Once β P is estimated, we calculate the ỹi (β̂ P ) using the known
values ∆(xi ) and solve (4) for β S . The whole procedure amounts to solving
the set of estimating equations U(β) = 0 where β = (β TS , β TP )T , U(β) =
(US (β)T , UP (β P )T )T , US (β) =
USi (β) =

n

ri yi + (1 −
n

UP i (β P ) = ri yi −

P

i

USi (β), UP (β P ) =

ri )h(β TP xP i

h(β TP xP i )

o

+ ∆(xi )) −

P

i

UP i (β P ),
o

h(β TS xSi )





xSi , 

(5)






xP i .

Pepe et al. (1994) derived a variance expression for β S assuming categorical xSi . To accommodate any form of xSi , we instead obtain standard errors
by the sandwich method, based on both estimating equations. The sandwich
 

estimator of var β̂ is
V = B−1 CB−T
where B = −dU/dβ evaluated at β = β̂, C =

(6)
P

i

Ui (β̂)Ui (β̂)T and Ui (β) =

(USi (β)T , UP i (β P )T )T . B and C are given in Section A of the Supplementary
Materials.

2.2

Estimation using two linear regressions

A special case arises if there are no auxiliary variables, so xP i = xSi for
all i, and h(.) is the identity function, as in linear regression. Then we can
8

rearrange (5) to give
n

o

USi (β) − UP i (β P ) = (1 − ri )∆(xi ) − (β S − β P )T xSi xSi .
Thus (β S −β P ) may be estimated by linear regression of (1−ri )∆(xi ) on xSi .
This estimate is uncorrelated with β̂ P because cov (USi (β) − UP i (β P ), UP i (β P )) =




0. This gives a direct way to estimate β S , and its variance var β̂ P +




var β̂ S − β̂ P , from standard linear regressions.
In particular, consider a two-arm trial with no covariates, and write the
coefficients of zi in (1) and (3) as βSz and βP z . Then βSz −βP z is estimated as
the difference between arms in the mean of (1 − ri )∆(xi ), which is a1 δ1 − a0 δ0
where aj is the proportion of missing data in arm j = 0, 1 and δj is the average
of ∆(xi ) over individuals with missing data in arm j = 0, 1. Therefore the
estimated parameter of interest is β̂ Sz = β̂ P z + a1 δ1 − a0 δ0 as in White et al.
(2007); the same result can be derived in other ways.

3

Equivalence to standard procedures

We now consider two special cases which can be fitted by standard procedures: (1) when MAR is assumed and there are no auxiliary variables, so
incomplete cases contribute no information and the standard procedure is an
analysis of complete cases, and (2) when ‘missing = failure’ is assumed for a
binary outcome, so the standard procedure is to replace missing values with
9

failures. Our aim is that point estimates, standard errors and confidence
intervals produced by the mean score procedure should agree exactly with
those produced by the standard procedures in these cases.
Equality of point estimates is easy to see. In case (1), we have ∆(xi ) = 0
and xP i = xSi for all i, so if β P solves UP (β P ) = 0 then β = (β P , β P ) solves
US (β) = 0. In case (2), ‘missing = failure’ can be expressed as ∆(xi ) = −∞
for all i, so the mean score procedure gives yi∗ = 0 whenever ri = 0, and
solving US (β S , β P ) = 0 gives the same point estimate as replacing missing
values with failures.
Exact equality of variances between mean score and standard procedures
depends on which finite sample corrections (if any) are applied. Many such
corrections have been proposed to reduce the small-sample bias of the sandwich variance estimator and to improve confidence interval coverage (Kauermann and Carroll, 2001; Lu et al., 2007). Here, we assume that the standard
procedures use the commonly used small-sample correction factor for the




sandwich variance vd
ar β̂ S = f V where f = n/(n − p∗ ), n is the sample size
and p∗ is the number of regression parameters (in linear regression) or 1 (in
other GLMs): this is for example the default in Stata (StataCorp, 2011).
Exact equality of confidence intervals between mean score and standard
procedures additionally depends on the distributional assumptions used to
construct confidence intervals. Here, we assume that standard procedures for
10

linear regression construct confidence intervals from the t distribution with
n − p∗ degrees of freedom, and that standard procedures for other GLMs
construct confidence intervals from the Normal distribution.
With missing data, we propose using the same small-sample correction
factor, distributional assumptions and degrees of freedom, but replacing n by
an effective sample size neff as shown below. Thus we propose forming confidence intervals for linear regression by assuming β̂ S ∼ tneff −p∗


neff
∗V
eff −p

and for other GLMs by assuming β̂ S ∼ N β S , n

3.1



n
β S , n eff
∗V
eff −p





.

Full sandwich method

For the full sandwich method of Section 2.1, we propose computing neff as
nobs + (Imis /Imis∗ )nmis where Imis is the influence of the individuals with
missing values, and Imis∗ is the same individuals’ influence if the missing
values had been observed. The comparison of the same individuals is crucial,
because missing individuals, if observed, would have different influence from
observed individuals.
To determine Imis , we consider weighted estimating equations

P

i

wi Ui (β) =

0 with solution β w . Differentiating with respect to w = (w1 , . . . , wn )T at
β)
w = 1 gives UT + dU(
dβ

dβ w
dw

= 0 where U is a n × (pS + pP ) matrix with ith

11

row Ui (β). Since also

dU(β )
dβ

= −B, we get
dβ w
= −B−1 UT .
dw

(7)

We now define the influence of observation i as
Imis,i =

dβ wS T −1 dβ wS
VS
dwi
dwi

where the S subscript denotes the elements corresponding to β S . Hence we
define the influence of the individuals with missing values as Imis =

P

i (1

−

Ri )Imis,i .
To determine Imis∗ , we let β ∗wS be the (unknown) parameter estimate
that would be obtained if the complete data had been observed, following
the pattern-mixture model (3). In this case the influence would be
∗
Imis,i

dβ ∗wS T −1 dβ ∗wS
VS
.
=
dwi
dwi

∗
over the disWe define the “full-data influence” as the expectation of Imis,i

tribution of the complete data given the observed data, under the pattern∗

mixture model (3). From (7) we get

dβ w
dwi

n

o

T
∗T
∗
= −B−1 U∗T
i and Ui = yi − h(β s xSi ) xSi ,

so that
E
and E

n

yi∗

−

h

∗
Imis,i

i

=E

o2
h(β Ts xS )

n

yi∗

−

o2
h(β Ts xSi )



−1 −1
xTSi B−T
SS VS BSS xSi

(8)



is evaluated as the squared residual plus the resid-

ual variance from model (1) for individual i. Finally, Imis∗ =
12

P

h

i

∗
i (1−ri )E Imis,i .

In case (1), Imis = 0 so neff = nobs , as in standard analysis. In case (2),
Imis∗ = Imis so neff = n, again as in standard analysis.

3.2

Two linear regressions method

For the two linear regressions method of Section 2.2, the small-sample correction to the variance is naturally applied separately to each variance in








vd
ar β̂ P + vd
ar β̂ S − β̂ P = Vsmall say. We can derive the corresponding variances without small-sample correction as

nobs −p
vd
ar
nobs







Vlarge say. To compute neff , we use the heuristic Vsmall ≈
hence we estimate neff by solving |Vsmall | =






β̂ P + n−p
vd
ar β̂ S − β̂ P =
n

neff
neff −p

p

neff
V
,
neff −p large

and

|Vlarge |. In case (1),



vd
ar β̂ S − β̂ P = 0 so neff = nobs , as in standard analysis. Case (2) does not
apply to linear regression.

4

Example: QUATRO trial

The QUATRO trial (Gray et al., 2006) was a randomised controlled trial in
people with schizophrenia, to evaluate the effectiveness of a patient-centred
intervention to improved adherence to prescribed antipsychotic medications.
The trial included 409 participants in four European centres. The primary
outcome, measured at baseline and 1 year, was participants’ quality of life,
expressed as the mental health component score (MCS) of the SF-36 (Ware,
13

1993). The MCS is designed to have mean 50 and standard deviation 10 in
a standard population, and a higher MCS score implies a better quality of
life. The data are summarised in Table 1.
We first estimate the intervention effect on MCS, adjusted for baseline
MCS and centre. Thus in the substantive model (1), h(.) is the identity
link, yi is MCS at 1 year for participant i, and xSi is a vector containing
1, randomised group zi (1 for the intervention group and 0 for the control
group), baseline MCS and dummy variables for three centres. We have no
xAi or xRi .
We also estimate the intervention effect on a binary variable, MCS dichotomised at an arbitrary value of 40, where for illustration we use baseline
MCS and dummy variables for centre as auxiliary variables xAi . Thus in the
substantive model (1), h(η) = 1/{1 + exp (−η)} is the inverse of the logit
link, yi is dichotomised MCS at 1 year for participant i, and xSi is a vector
containing 1 and zi .
In both analyses, we replace the 23 missing values of baseline MCS with
the mean baseline MCS: while such mean imputation is not valid in general,
it is appropriate and efficient in the specific case of estimating intervention effects with missing baseline covariates in randomised trials (White and
Thompson, 2005; Groenwold et al., 2012).
As expected (results not shown), the point estimate, standard error and
14

confidence interval from the mean score method agree exactly with standard
methods under MAR and under missing=failure, using the small-sample corrections of Section 3.
We consider three sets of sensitivity analyses using the mean score method
around a MAR assumption, with departures from MAR (1) in the intervention arm only (∆(xi ) = δzi ), (2) in both arms (∆(xi ) = δ), or (3) in the
control arm only (∆(xi ) = δ(1 − zi )) (White et al., 2012). For the quantitative outcome, the investigators suggested that the mean of the missing data
could plausibly be lower than the mean of the observed data by up to 10 units
(equal to nearly one standard deviation of the observed data), so we allow
δ to range from 0 to -10 (Jackson et al., 2010). The investigators were not
asked about missing values of the dichotomised outcome, so for illustrative
purposes we allow δ to range from 0 to -6, which is close to “missing=failure”.
Figure 1 shows the results of these sensitivity analysis using the two linear
regressions method for the quantitative outcome (upper panel) and using the
full sandwich method for the dichotomised outcome (lower panel). Results
are more sensitive to departures from MAR in the intervention arm because
there are more missing data in this arm. However, the finding of a nonsignificant intervention effect is unchanged over these ranges of sensitivity
analyses. This means that the main results of the trial are robust even to
quite strong departures from MAR.
15

Figure 2 shows the effective sample size for these two analyses. Effective sample size increases from 367 at MAR for both analyses. For the
dichotomised outcome it is near the total sample size of 409 when δ = −6 in
both arms. For the quantitative outcome it does not pass 370 because the
range of δ is more moderate.

5

Simulation study

We report a simulation study aiming (i) to evaluate the performance of the
mean score method when it is correctly specified, (ii) to compare the mean
score method with alternatives, and (iii) to explore the impact of the incompatibility of models (1) and (3). We assume the sensitivity parameters ∆(xi )
are correctly specified.

5.1

Data generating models

We generate data under four data generating models (DGMs), each with
four choices of parameters. We focus on the case of a binary outcome. In
DGMs 1-3, we generate data under a pattern-mixture model. In DGM 1,
there are no baseline covariates, so xi = xSi = xP i = (1, zi ). We generate a treatment indicator zi ∼ Bern(0.5); a missingness indicator ri with
logit P (ri |zi ) = α1 + αz zi ; and a binary outcome yi following model (3) with
16

logit P (yi |zi , ri ) = βP 1 + βP z zi + βP r (1 − ri ). The substantive model (1) is
then logit P (yi |zi ) = βS1 + βSz zi . Because this substantive model contains
only a single binary covariate, it is saturated and cannot be mis-specified.
Therefore both substantive model and pattern-mixture model are correctly
specified
DGM 2 extends DGM 1 by including xi ∼ N (0, 1) as a single baseline
covariate independent of zi , so xi = (1, xi , zi ). The missingness indicator
follows logit P (ri |xi , zi ) = α1 + αx xi + αz zi and the binary outcome follows
model (3) with logit P (yi |xi , zi , ri ) = βP 1 + βP x xi + βP z zi + βP r (1 − ri ). The
substantive model is as in DGM 1 with xSi = (1, zi ), and xAi = (xi ) is an
auxiliary variable in the analysis. Thus the substantive model and patternmixture model are again both correctly specified.
DGM 3 is identical to DGM 2, but now xi is included in the substantive
model, which is therefore logit P (yi |xi , zi ) = βS1 + βSx xi + βSz zi with xSi =
xP i = (1, xi , zi ). Now the substantive model is incorrectly specified while the
pattern-mixture model remains correctly specified.
DGM 4 is a selection model. Here zi and xi are generated as in DGM
2, then yi is generated following the substantive model logit P (yi |xi , zi ) =
βS1 + βSx xi + βSz zi and ri is generated following logit P (ri |xi , zi , yi ) = α1 +
αx xi + αz zi + αy yi . Here the substantive model is correctly specified while
the pattern mixture model is mis-specified.
17

For the parameter values, we consider scenarios a-d for each DGM. In
scenario a, the sample size is nobs = 500; the missingness model has αx =
αz = αy = 1 and we choose α1 to fix πobs = P (r = 1) = 0.75; and the pattern
mixture model has βP 1 = 0, βP x = βP z = 1, βP r = −1. (αx and βP x are
ignored in DGM 1, αy is ignored in DGM 1-3 and βP r is ignored in DGM
4.) Scenarios b-d vary scenario a by setting nobs = 2000, πobs = 0.5, and
βP r = −2 respectively. 1000 data sets were simulated in each case. Table 2
summarises the simulation design.

5.2

Analysis methods

The mean score (MS) method is implemented as described in sections 2 and
3, with logit link. xi is used as an auxiliary in DGM 2. In DGM 1-3, ∆(xi )
is taken to equal the known value −βP r for all individuals; in DGM 4, ∆(xi )
is not known but (for the purposes of the simulation study) is estimated by
fitting the pattern-mixture model to a data set of size 1,000,000 before data
deletion.
The MS method is compared with analysis of data before data deletion
(Full); analysis of complete cases (CC), which wrongly assumes MAR; and
two alternative methods that allow for MNAR, multiple imputation (MI)
and selection model with inverse probability weighting (SM).

18

In the MI approach (Rubin, 1987; White et al., 2011b), the imputation
model is equation (3), and data are imputed with an offset ∆(xi ) in the
imputation model. The number of imputations is fixed at 30.
In the SM approach, we use the response model logit p(ri = 1|yi , xi ) =
αT xP i +∆∗ (xi )yi where the sensitivity parameter ∆∗ (xi ) expresses departure
from MAR as the log odds ratio for response per 1-unit change in yi . In DGM
4, ∆∗ (xi ) is taken to equal the known value αy ; in DGM 1-3, ∆∗ (xi ) is not
known but (for the purposes of the simulation study) is estimated by fitting
the selection model to a data set of size 1,000,000 before data deletion. The
parameters α cannot be estimated by standard methods, since some yi are
missing, so we use a weighted estimating equation which does not involve the
missing yi ’s (Rotnitzky et al., 1998; Dufouil et al., 2004; National Research
Council, 2010):
(
X

xP i

i

)

ri
− 1 = 0.
T
h(α xP i + ∆∗ (xi )yi )

The substantive model is then fitted to the complete cases with stabilised
weights p̂(ri = 1|xSi )/p̂(ri = 1|yi , xi ), where p̂(ri = 1|xSi ) is estimated by the
same procedure as p̂(ri = 1|yi , xi ) but with no ∆∗ (xi )yi or xAi terms (Robins
et al., 2000). Variances are computed by the sandwich variance formula,
ignoring uncertainty in α̂.

19

5.3

Estimand

The estimand of interest is the coefficient βSz in the substantive model. It is
computed by fitting the substantive model to the data set of size 1,000,000
before data deletion. We explore bias, empirical and model-based standard
errors, and coverage of estimates β̂Sz .

5.4

Results

Results are shown in Table 3. CC is always biased, often inefficient, and
poorly covering. Small bias (at most 3% of the true value) is observed in
the “Full” analysis (i.e. before data deletion) in some settings: this is a
small-sample effect (Nemes et al., 2009), since as noted above, the true value
and “Full” are calculated in the same way with large and small samples
respectively. Taking “Full” as a gold standard, MS, MI and SM methods
all have minimal bias (at most 2% of true values). Precisions of MS and
MI are similar, with SM slightly inferior in DGM 2. Coverages are near
95%, with some over-coverage for SM in some settings, as a consequence of
slightly overestimated standard errors (results not shown) due to ignoring
uncertainty in α̂ (Lunceford and Davidian, 2004). The performance of MS
is not appreciably worse when the selection model is mis-specified (DGM 3)
or when the pattern-mixture model is mis-specified (DGM 4). Computation

20

times for MI are 15-18 times longer than for MS, which is 10-30% longer than
SM.

6

Discussion

We have proposed a mean score method which works well when the sensitivity
parameters are known. In practice, of course, the sensitivity parameters are
unknown, and a range of values will be used in a sensitivity analysis.
The main practical difficulty in implementing any principled sensitivity
analysis is choosing the value(s) of the sensitivity parameters. This is a subjective process requiring subject-matter knowledge and is best done by discussion between statisticians and suitable ‘experts’, typically the trial investigators. By using the pattern-mixture model, we use a sensitivity parameter
∆(xi ) that is easier to communicate with ‘experts’ than the corresponding
parameters in selection models or shared parameter models. The procedure
has been successfully applied in several trials (White, 2015). Special attention is needed to the possibility that ∆(xi ) varies between randomised
groups, because estimated treatment effects are highly sensitive to such variation (White et al., 2007). As with all aspects of trial analysis, plausible
ranges of the ∆(xi ) parameters should be defined before the data are collected or before any analysis. An alternative approach would report the

21

“tipping point”, the value of ∆(xi ) for which the main results are substantively affected, leaving the reader to make the subjective decisions about the
plausibility of more extreme values (Liublinska and Rubin, 2014). Presenting this information could be complex without subjective decisions about
the difference in ∆(xi ) between randomised groups. Effective methods are
therefore needed for eliciting sensitivity parameters.
Our method does not incorporate data on discontinuation of treatment,
unless this can be included as an auxiliary variable. Our method would
however be a suitable adjunct to estimation of effectiveness in a trial with
good follow-up after discontinuation of treatment. Further work is needed
to combine our sensitivity analysis with models for outcome before and after
discontinuation of treatment (Little and Yau, 1996). In a drug trial in which
follow-up ends on discontinuation of treatment, we see our method as estimating efficacy or a de jure estimand; if effectiveness or a de facto estimand
is required then post-discontinuation missing data in each arm may be imputed by the methods of Carpenter et al. (2013). Further work is needed to
perform a full sensitivity analysis in this setting.
Our use of parametric models makes our results susceptible to model misspecification, and indeed in many cases models (1) and (3) cannot both be
correctly specified except under MAR. However, the simulation study shows
that the impact of such model inconsistency is small relative to the impact of
22

assumptions about the missing data and the difficulty of knowing the values
of the sensitivity parameters.
We compared the mean score method with multiple imputation and inverse probability weighting. In the case of a quantitative outcome, MI can
be simplified by imputing under MAR and then adding the offset to the imputed data before fitting the substantive model and applying Rubin’s rules.
We could also impute under MAR and then use a weighted version of Rubin’s
rules to allow for MNAR (Carpenter et al., 2007). Both MI methods are subject to Monte Carlo error and so seem inferior to the mean score method. A
full likelihood-based analysis of the selection model would also be possible,
and a Bayesian analysis could directly allow for uncertainty about ∆(xi ) in
a single analysis (Mason et al., 2012). These alternative approaches are both
more computationally complex.
The proposed mean score method can be extended in various ways. We
have illustrated the method for departures from a MAR assumption, but it
can equally be used if the primary analysis with a binary outcome assumed
“missing = failure”, by varying ∆(xi ) from −∞ rather than from 0. The
method is also appropriate in observational studies, except that mean imputation for missing covariates is not appropriate in this context. The method
can be applied to a cluster-randomised trial as described in Section B of
the Supplementary Materials. Further work could allow the imputation and
23

substantive models to have different link functions, non-canonical links to
be used with suitable modification to US , and extension to trials with more
than two arms.
Sensitivity analysis should be more widely used to assess the importance
of departures from assumptions about missing data. The proposed mean
score approach provides data analysts with a fast and fully theoretically
justified way to perform the sensitivity analyses. It is implemented in a
Stata module rctmiss available from Statistical Software Components (SSC)
at https://ideas.repec.org/s/boc/bocode.html.

Supplementary materials
The supplementary materials give details of the sandwich variance in equation (6) and sketch an extension to clustered data.

Acknowledgements
IRW was funded by the Medical Research Council [Unit Programme number U105260558]. NJH was funded by NIH grant 5R01MH054693-12. JC
was funded by ESRC Research Fellowship RES-063-27-0257. We thank the
QUATRO Trial Team for access to the study data.

24

References
Carpenter, J. R., Kenward, M. G. and White, I. R. (2007) Sensitivity analysis
after multiple imputation under missing at random: a weighting approach.
Statistical Methods in Medical Research, 16(3), 259–275.
Carpenter, J. R., Roger, J. H. and Kenward, M. G. (2013) Analysis of longitudinal trials with protocol deviation: a framework for relevant, accessible
assumptions, and inference via multiple imputation. Journal of Biopharmaceutical Statistics, 23(3), 1352–71.
Daniels, M. J. and Hogan, J. W. (2000) Reparameterizing the Pattern Mixture Model for Sensitivity Analyses Under Informative Dropout. Biometrics, 56(4), 1241–1248.
Diggle, P. and Kenward, M. G. (1994) Informative drop-out in longitudinal
data analysis. Applied Statistics, 43(1), 49–93.
Dufouil, C., Brayne, C. and Clayton, D. (2004) Analysis of longitudinal
studies with death and drop-out: a case study. Statistics in Medicine,
23(14), 2215–2226.
Follmann, D. and Wu, M. (1995) An approximate generalized linear model
with random effects for informative missing data. Biometrics, 51, 151–168.

25

Gray, R., Leese, M., Bindman, J., Becker, T., Burti, L., David, A., Gournay,
K., Kikkert, M., Koeter, M., Puschner, B., Schene, A., Thornicroft, G.
and Tansella, M. (2006) Adherence therapy for people with schizophrenia:
European multicentre randomised controlled trial. The British Journal of
Psychiatry, 189(6), 508–514.
Groenwold, R. H. H., White, I. R., Donders, A. R. T., Carpenter, J. R.,
Altman, D. G. and Moons, K. G. M. (2012) Missing covariate data in
clinical research: when and when not to use the missing-indicator method
for analysis. Canadian Medical Association Journal, 184(11), 1265–1269.
Hedeker, D. R. and Gibbons, R. D. (2006) Longitudinal data analysis. John
Wiley and Sons.
Higgins, J. P., White, I. R. and Wood, A. M. (2008) Imputation methods for
missing outcome data in meta-analysis of clinical trials. Clinical Trials,
5(3), 225–239.
Jackson, D., White, I. R. and Leese, M. (2010) How much can we learn about
missing data?: an exploration of a clinical trial in psychiatry. Journal of
the Royal Statistical Society: Series A (Statistics in Society), 173(3), 593–
612.
Kaciroti, N. A., Schork, M. A., Raghunathan, T. and Julius, S. (2009) A
26

Bayesian sensitivity model for intention-to-treat analysis on binary outcomes with dropouts. Statistics in Medicine, 28(4), 572–585.
Kauermann, G. and Carroll, R. J. (2001) A note on the efficiency of sandwich covariance matrix estimation. Journal of the American Statistical
Association, 96(456), 1387–1396.
Kenward, M., Goetghebeur, E. and Molenberghs, G. (2001) Sensitivity analysis for incomplete categorical data. Statistical Modelling, 1(1), 31–48.
Kenward, M. G. (1998) Selection models for repeated measurements with
non-random dropout: an illustration of sensitivity. Statistics in Medicine,
17(23), 2723–2732.
Little, R. and Yau, L. (1996) Intent-to-treat analysis for longitudinal studies
with drop-outs. Biometrics, 52(4), 1324–1333.
Little, R. J. A. (1993) Pattern-mixture models for multivariate incomplete
data. Journal of the American Statistical Association, 88(421), 125–134.
Little, R. J. A. (1994) A class of pattern-mixture models for normal incomplete data. Biometrika, 81(3), 471–483.
Little, R. J. A. and Rubin, D. B. (2002) Statistical Analysis with Missing
Data. Hoboken, N. J.: Wiley, second edition.

27

Liublinska, V. and Rubin, D. B. (2014) Sensitivity analysis for a partially
missing binary outcome in a two-arm randomized clinical trial. Statistics
in Medicine, 33(24), 4170–4185.
Lu, B., Preisser, J. S., Qaqish, B. F., Suchindran, C., Bangdiwala, S. I.
and Wolfson, M. (2007) A comparison of two bias-corrected covariance
estimators for generalized estimating equations. Biometrics, 63(3), 935–
941.
Lunceford, J. K. and Davidian, M. (2004) Stratification and weighting via the
propensity score in estimation of causal treatment effects: a comparative
study. Statistics in Medicine, 23, 2937–2960.
Mason, A., Richardson, S., Plewis, I. and Best, N. (2012) Strategy for modelling nonrandom missing data mechanisms in observational studies using
Bayesian methods. Journal of Official Statistics, 28(2), 279–302.
Mavridis, D., White, I. R., Higgins, J. P. T., Cipriani, A. and Salanti, G.
(2015) Allowing for uncertainty due to missing continuous outcome data in
pair-wise and network meta-analysis. Statistics in Medicine, 34, 721–741.
National Research Council (2010) The Prevention and Treatment of Missing
Data in Clinical Trials. Washington, DC: Panel on Handling Missing Data

28

in Clinical Trials. Committee on National Statistics, Division of Behavioral
and Social Sciences and Education. The National Academies Press.
Nemes, S., Jonasson, J. M., Genell, A. and Steineck, G. (2009) Bias in odds
ratios by logistic regression modelling and sample size. BMC Medical Research Methodology, 9(1), 56.
Pepe, M. S., Reilly, M. and Fleming, T. R. (1994) Auxiliary outcome data
and the mean score method. Journal of Statistical Planning and Inference,
42, 137–160.
Reilly, M. and Pepe, M. S. (1995) A mean score method for missing and
auxiliary covariate data in regression models. Biometrika, 82(2), 299–314.
Robins, J. M., Hernan, M. A. and Brumback, B. (2000) Marginal structural
models and causal inference in epidemiology. Epidemiology, 11, 550–560.
Rogers, W. H. (1993) Regression standard errors in clustered samples. Stata
Technical Bulletin, 13, 19–23.
Rotnitzky, A., Robins, J. M. and Scharfstein, D. O. (1998) Semiparametric
regression for repeated outcomes with nonignorable nonresponse. Journal
of the American Statistical Association, 93(444), 1321–1339.
Roy, J. (2003) Modeling longitudinal data with nonignorable dropouts using
a latent dropout class model. Biometrics, 59(4), 829–836.
29

Rubin, D. B. (1987) Multiple Imputation for Nonresponse in Surveys. New
York: John Wiley and Sons.
Scharfstein, D., McDermott, A., Olson, W. and Wiegand, F. (2014)
Global sensitivity analysis for repeated measures studies with informative dropout: a fully parametric approach. Statistics in Biopharmaceutical
Research, 6(4), 338–348.
Scharfstein, D. O., Daniels, M. J. and Robins, J. M. (2003) Incorporating
prior beliefs about selection bias into the analysis of randomized trials with
missing outcomes. Biostatistics, 4, 495–512.
StataCorp (2011) Stata Statistical Software: Release 12. College Station,
TX: StataCorp LP.
Ware, J. E. (1993) SF-36 Health Survey: Manual and Interpretation Guide.
Boston: The Health Institute, New England Medical Center.
Ware, J. H. (2003) Interpreting incomplete data in studies of diet and weight
loss. New England Journal of Medicine, 348(21), 2136–2137.
West, R., Hajek, P., Stead, L. and Stapleton, J. (2005) Outcome criteria in
smoking cessation trials: proposal for a common standard. Addiction, 100,
299–303.

30

White, I. R. (2015) Sensitivity analysis: The elicitation and use of expert
opinion. In Handbook of Missing Data Methodology (Eds G. Molenberghs,
G. Fitzmaurice, M. G. Kenward, A. Tsiatis and G. Verbeke). Chapman
and Hall.
White, I. R. and Thompson, S. G. (2005) Adjusting for partially missing
baseline measurements in randomized trials. Statistics in Medicine, 24(7),
993–1007.
White, I. R., Carpenter, J., Evans, S. and Schroter, S. (2007) Eliciting and
using expert opinions about dropout bias in randomized controlled trials.
Clinical Trials, 4(2), 125–139.
White, I. R., Higgins, J. P. T. and Wood, A. M. (2008) Allowing for uncertainty due to missing data in meta-analysis - Part 1: Two-stage methods.
Statistics in Medicine, 27(5), 711–727.
White, I. R., Horton, N. J., Carpenter, J. and Pocock, S. J. (2011a) Strategy
for intention to treat analysis in randomised trials with missing outcome
data. British Medical Journal, 342, d40.
White, I. R., Wood, A. and Royston, P. (2011b) Tutorial in biostatistics:
Multiple imputation using chained equations: issues and guidance for practice. Stat Med, 30, 377–399.
31

White, I. R., Carpenter, J. and Horton, N. J. (2012) Including all individuals
is not enough: lessons for intention-to-treat analysis. Clinical Trials, 9,
396–407.
Wood, A. M., White, I. R. and Thompson, S. G. (2004) Are missing outcome
data adequately handled? A review of published randomized controlled
trials in major medical journals. Clinical Trials, 1(4), 368–376.

32

Table 1: QUATRO trial: data summary.
Intervention
(n=204)
Centre
Amsterdam (%)
50 (25%)
Leipzig (%)
49 (24%)
London (%)
45 (22%)
Verona (%)
60 (29%)
MCS at baseline Mean (SD)
38.4 (11.2)
Missing (%)
13 (6%)
MCS at 1 year
Mean (SD)
40.2 (12.0)
> 40 (%)
99 (57%)
Missing (%)
29 (14%)

Control
(n=205)
50 (24%)
48 (23%)
47 (23%)
60 (29%)
40.1 (12.1)
10 (5%)
41.3 (11.5)
104 (54%)
13 (6%)

Table 2: Simulation study: data generating models.
DGM
1
2
3
4

Type
Pattern-mixture
Pattern-mixture
Pattern-mixture
Selection

PM correct?
Yes
Yes
Yes
No

SM correct?
Yes
Yes
No
Yes

Auxiliary variables?
No
Yes
No
No

Parameters
Description
a
Base case
b
Larger sample size
c
More missing data
d
Larger departure from MAR

33

33

Table 3: Simulation results: bias, empirical standard error and coverage of nominal 95% confidence interval for
methods Full (before data deletion), CC (complete cases analysis), MS (mean score), MI (multiple imputation), SM
(selection model + IPW). Error denotes maximum Monte Carlo error.
DGM

34

1a
1b
1c
1d
2a
2b
2c
2d
3a
3b
3c
3d
4a
4b
4c
4d
Error

Full
0.007
0.004
0.007
0.000
0.009
0.004
0.005
0.011
0.016
0.003
0.012
0.019
0.005
0.006
0.008
0.026
0.007

Bias
CC
MS
-0.128 0.010
-0.130 0.007
-0.167 0.018
-0.177 0.003
-0.151 0.010
-0.157 0.003
-0.206 0.002
-0.175 0.011
-0.110 0.021
-0.127 0.003
-0.160 0.015
-0.163 0.020
-0.148 0.014
-0.147 0.015
-0.203 0.004
-0.118 0.045
0.011 0.011

MI
SM
Full
0.010 0.009 0.191
0.006 0.006 0.095
0.018 0.018 0.183
0.004 0.004 0.187
0.010 0.005 0.178
0.003 -0.001 0.092
-0.001 -0.004 0.185
0.010 0.007 0.186
0.020 0.015 0.208
0.003 -0.003 0.106
0.012 0.018 0.217
0.020 0.020 0.216
0.012 0.006 0.218
0.015 0.008 0.107
0.002 0.011 0.213
0.045 0.036 0.211
0.011 0.011 0.005

Empirical
CC
MS
0.223 0.218
0.113 0.111
0.267 0.258
0.223 0.203
0.219 0.203
0.109 0.101
0.296 0.256
0.227 0.200
0.249 0.245
0.122 0.120
0.327 0.316
0.254 0.238
0.255 0.251
0.124 0.122
0.345 0.336
0.250 0.246
0.008 0.008

SE
MI
0.219
0.112
0.259
0.204
0.203
0.101
0.257
0.200
0.246
0.120
0.318
0.239
0.252
0.122
0.337
0.248
0.008

SM
0.218
0.111
0.258
0.203
0.209
0.104
0.290
0.202
0.249
0.121
0.323
0.249
0.253
0.123
0.343
0.249
0.008

Full
95.0
94.3
96.5
95.1
96.3
94.6
94.9
95.0
96.5
94.8
95.1
94.2
94.2
94.0
94.2
94.2
0.8

Coverage
CC MS MI
90.6 95.1 95.3
76.7 93.8 93.4
91.9 95.9 95.4
88.3 95.6 95.5
90.1 95.4 95.4
71.5 95.1 95.3
88.1 95.0 94.7
86.1 94.5 95.2
92.7 95.4 94.7
81.8 95.3 95.2
90.3 95.0 95.1
88.6 94.5 94.6
90.0 94.8 94.5
77.1 94.0 94.1
90.8 94.2 94.3
91.2 94.3 94.5
1.4 0.8 0.8

SM
95.8
94.6
96.6
97.1
97.0
96.9
97.2
98.2
95.3
95.3
95.4
95.3
95.0
94.6
95.0
94.8
0.7

34

Figure 1: QUATRO trial: sensitivity analysis for the estimated intervention
effect on the MCS (with 95% confidence interval) over a range of departures
from MAR.

-2
-4
-10

-5

0

0

.5

1

MCS score > 40: log odds ratio

-.5

Estimated intervention effect

0

2

MCS score: mean difference

-6

-4

-2

Sensitivity parameter δ
Sensitivity analysis method applied to
Intervention arm
Both arms
Control arm

35

0

Figure 2: QUATRO data: effective sample size in sensitivity analysis.

-5

0

MCS score > 40: log odds ratio

400

420

366

-10

380

Effective sample size

368

370

MCS score: mean difference

-6

-4

-2

Sensitivity parameter δ
Sensitivity analysis method applied to
Intervention arm
Both arms
Control arm

36

0

Supplementary materials
A

Details of mean score variance

The sandwich variance (6) is computed from B and C, where B has components
BSS = −dUS /dβ S =

X

h0 (β̂ S xSi )xSi xTSi

i

BSP

= −dUS /dβ P

=−

X

(1 − ri )h0 (β̂ P xP i + ∆i )xSi xTP i

i

BP S = −dUP /dβ S = 0
BP P = −dUP /dβ P =

X

ri h0 (β̂ P xP i )xP i xTP i

i

and C has components
CSS =

X

CSP = CTP S =

X

CP P =

X

e2Si xSi xTSi

i

eSi eP i xSi xTP i

i

e2P i xP i xTP i

i
T

T

where eSi = yi∗ (β̂ P ) − h(β̂ S xSi ) and eP i = ri {yi − h(β̂ P xP i )}.

B

Modifications for clustered data

If data are clustered, as in a cluster-randomised trial, we need to modify the
variance calculations in Sections 2.1 and 2.2 and the small-sample corrections
in Sections 3.1 and 3.2. Let m be the total number of clusters, mobs be
the number of clusters with at least one observed outcome, and mmis =
1

m − mobs be the number of clusters with no observed outcome. Let the data
be subscripted by cluster membership c = 1, . . . , m as well as individual i.
For the full sandwich variance method of Section 2.1, we only need to redefine the matrix C =

P

c

Uc (β̂)Uc (β̂)T where Uc (β̂) =

P

i

Uci (β̂) (Rogers,

1993).
For the two linear regressions method of Section 2.2, we similarly take








var β̂ P and var β̂ S − β̂ P as clustered sandwich variances.
For the small-sample methods of Section 3, we assume the standard methods use a small-sample correction factor f =

n−1 m
,
n−p∗ m−1

and use m − 1 degrees

of freedom for linear regression (StataCorp, 2011). We replace n and m by
neff and meff , calculated by the two methods explained below.
For the full sandwich variance method, we compute neff as in Section 3.1,
and compute meff = mobs + (Imis /Imis∗ )mmis .
For the two linear regressions method, the variance with small-sample






ar β̂ S − β̂ P
correction is (as before) vd
ar β̂ P + vd



sponding variance without small-sample correction is
n−p m−1
vd
ar
n−1 m



β̂ S − β̂ P



= Vsmall . The correnobs −p mobs −1
vd
ar
nobs −1 mobs

= Vlarge . The heuristic Vsmall ≈

leads to the equation |Vsmall | =



p

neff
neff −p





β̂ P +

neff −1 meff
V
neff −p meff −1 large

|Vlarge |. However, we have two un-

knowns neff and meff , so we take a second equation representing the variance
with small-sample correction only for the number of clusters:

2

mobs −1
vd
ar
mobs





β̂ P +

m−1 d
var
m





β̂ S − β̂ P = Vlargen say, with the heuristic Vsmall ≈

the second equation |Vsmall | =



meff
meff −1

p

and

|Vlargen |. We solve the second equa-

tion for meff and then the first equation for neff .

3

meff
V
meff −1 largen

