# A weighting method for simultaneous adjustment for confounding and joint
exposure-outcome misclassifications

Bas B.L. Penning de Vries1,* Maarten van Smeden1 and Rolf H.H. Groenwold1,2
1Department of Clinical Epidemiology<EMAIL_ADDRESS>2Department of Biomedical Data Sciences
Leiden University Medical Center  PO Box 9600  2300 RC  The Netherlands

###### Abstract

Joint misclassification of exposure and outcome variables can lead to
considerable bias in epidemiological studies of causal exposure-outcome
effects. In this paper, we present a new maximum likelihood based estimator
for the marginal causal odd-ratio that simulaneously adjusts for confounding
and several forms of joint misclassification of the exposure and outcome
variables. The proposed method relies on validation data for the construction
of weights that account for both sources of bias. The weighting estimator,
which is an extension of the exposure misclassification weighting estimator
proposed by Gravel and Platt (Statistics in Medicine, 2018), is applied to
reinfarction data. Simulation studies were carried out to study its finite
sample properties and compare it with methods that do not account for
confounding or misclassification. The new estimator showed favourable large
sample properties in the simulations. Further research is needed to study the
sensitivity of the proposed method and that of alternatives to violations of
their assumptions. The implementation of the estimator is facilitated by a new
R function in an existing R package.

###### keywords:

Causal inference; confounding; propensity scores, inverse probability
weighting; joint exposure and outcome misclassification; validation data.

##  1 Introduction

In epidemiological research on causal associations between a particular
exposure and a certain outcome, erroneous information on either or both of
these variables poses a serious methodological obstacle in making valid
inferences. In particular, joint misclassification of exposure and outcome can
lead to considerable bias of standard causal effect estimators, with direction
and magnitude depending on various factors, including the misclassification
mechanism and the direction and magnitude of the true effect (Kristensen,
1992; Brenner et al., 1993; Vogel et al., 2005; Jurek et al., 2008;
VanderWeele and Hernán, 2012; Brooks et al., 2018).

Exposure and outcome misclassification is typically categorised according to
two separate properties: whether or not the misclassification is differential
and whether or not it is dependent relative to some covariate vector $L$
containing patient characteristics (Kristensen, 1992; VanderWeele and Hernán,
2012). Joint misclassification of exposure and outcome is said to be
_nondifferential_ if (1) the sensitivity and specificity of exposure
classification are constant across all categories of the (true) outcome given
$L$ and (2) the sensitivity and specificity of outcome classification are
constant across all categories of the (true) exposure given $L$; otherwise it
is _differential_. Misclassification is said to be _independent_ if the joint
probability of any exposure and outcome classification given any true exposure
and outcome categories and $L$ can be factored into the product of the
corresponding probabilities for exposure and outcome separately; otherwise, it
is _dependent_. In Dawid’s notation (1979), that is, if true exposure level
$A$ and true outcome $Y$ are (potentially mis)classified as $B$ and $Z$,
respectively, misclassification is nondifferential if and only if
$B\mathrel{\perp\\!\\!\\!\perp}Y|A,L$ and
$Z\mathrel{\perp\\!\\!\\!\perp}A|Y,L$ and independent if and only if
$Z\mathrel{\perp\\!\\!\\!\perp}B|Y,A,L$.

Epidemiological research hampered by joint misclassification of some type is
likely voluminous (Brooks et al., 2018). Examples of studies affected by
exposure and outcome misclassification can be found, for example, in the
literature on the causal effects of drug use, which is largely based on
routinely collected data, where exposures are typically operationalised on the
basis of prescription records and where outcomes are often self-reported
(Marcum et al., 2013; Culver et al., 2012; Leong et al., 2013; Ni et al.,
2017). Departures from differentiality are likely, particularly when it
concerns prescription-only medication. In applied epidemiological research,
misclassification or some of its potential consequences are often ignored
(Jurek et al., 2006; Brakenhoff et al., 2018). The assertion often made in the
discussion of study results that observed measures of association are biased
toward the null under nondifferentiality, for example, is not generally true
unless additional conditions are presupposed (Brenner et al., 1993; Brooks et
al., 2018).

Methods to adjust for misclassification rely on additional information that
can be used to estimate or correct for bias. One potential source of
information is validation data obtained through supposedly infallible
measurement. Recently, Gravel and Platt (2018) proposed an inverse probability
weighting (IPW) method to simultaneously address confounding and outcome
misclassification by means of internal validation data. In what follows, we
propose an extension of this method to allow for confounding adjustment and
joint exposure and outcome misclassification. This flexible estimator allows
for the misclassifications to be dependent, differential or both. In Section
2, inverse probability weights for confounding and joint misclassification are
introduced through a hypothetical study based on the illustrative example of
Gravel and Platt (2018). Section 3 details methods for estimation of the
various components of the proposed weights. In Section 4, we describe a series
of Monte Carlo simulations that were used to study properties of the proposed
method in finite samples. We conclude with a summary and discussion of our
findings in context of the existing literature.

##  2 Data distribution for illustration and development of weighting method

We first consider the data and setting described by Gravel and Platt (2018)
and suppose that Table 1 represents a simple random (i.i.d.) sample from (or
that its cell counts are proportional to the respective densities in) the
population of interest. This illustration is based on a cohort study on the
association between post-myocardial infarction statin use ($A$) and the 1-year
risk of reinfarction ($Y$). In what follows, we will refer to this example as
the ‘reinfarction example’.

Throughout we take the counterfactual framework for causal inference, formal
accounts of which are given for example by Neyman et al. (1935), Rubin (1974),
Holland (1986), Holland (1988) and Pearl (2009). The interest, we suppose,
lies in estimating a function of counterfactuals, in particular the causal
marginal odds ratio OR,

OR
$\displaystyle=\frac{\mathbb{E}[Y(1)]/(1-\mathbb{E}[Y(1)])}{\mathbb{E}[Y(0)]/(1-\mathbb{E}[Y(0)])},$
(1)

where $Y(0),Y(1)$ denote the counterfactual outcomes for hypothetical
interventions setting $A$ to $0$ and $1$, respectively.

Table 1: Cross-classification of the reinfarction data for 33,007 individuals as given by Gravel and Platt (2018). | $L=0$ |  | $L=1$
---|---|---|---
| $A=0$ | $A=1$ |  | $A=0$ | $A=1$
$Y=0$ | $11602$ | $13116$ |  | $1302$ | $5363$
$Y=1$ | $890$ | $589$ |  | $49$ | $96$

###  2.1 No misclassification

Under conditional exchangeability given $L$ (i.e.,
$(Y(0),Y(1))\mathrel{\perp\\!\\!\\!\perp}A|L$), consistency ($Y(a)=Y$ if
$A=a$) and positivity ($\Pr(A=a|L=l)>0$ for $a=0,1$ and all $l$ in the support
of $L$), the odds ratio of (1) can be expressed in terms of ‘observables’
(meaning, here, variables that would be observed had there been no measurement
error) as follows:

OR
$\displaystyle=\frac{\mathbb{E}[WY|A=1]/(1-\mathbb{E}[WY|A=1])}{\mathbb{E}[WY|A=0]/(1-\mathbb{E}[WY|A=0])}$
(2)

with weights $W$ defined as the inverse probability of the allocated exposure
level $A$ given $L$ (i.e., the inverse propensity score) multiplied by the
prevalence of the allocated exposure level $A$ (i.e., $W=\Pr(A)/\Pr(A|L)$;
Appendix I).

Replacing components of the right-hand side of (2) with sample analogues, we
obtain the following estimator for the setting where $L$ is binary:

$\displaystyle\widehat{\text{OR}}$
$\displaystyle:=\frac{\widehat{\mathbb{E}}[\widehat{W}Y|A=1]/(1-\widehat{\mathbb{E}}[\widehat{W}Y|A=1])}{\widehat{\mathbb{E}}[\widehat{W}Y|A=0]/(1-\widehat{\mathbb{E}}[\widehat{W}Y|A=0])}$
$\displaystyle=\frac{(\widehat{W}_{10}n_{110}+\widehat{W}_{11}n_{111})/(n_{110}+n_{111}+n_{010}+n_{011}-\widehat{W}_{10}n_{110}-\widehat{W}_{11}n_{111})}{(\widehat{W}_{00}n_{100}+\widehat{W}_{01}n_{101})/(n_{100}+n_{101}+n_{000}+n_{001}-\widehat{W}_{00}n_{100}-\widehat{W}_{01}n_{101})},$
(3)

where $n_{yal}$ denotes the number of subjects with $Y=y$, $A=a$, $L=l$ and
where $\widehat{W}_{al}$ is the product of the proportion of subjects in the
sample with $A=a$ and the inverse of the proportion of subjects with $A=a$
among those with $L=l$. For the data in Table 1, we obtain
$\widehat{\text{OR}}\approx 0.573$. The corresponding crude odds ratio (i.e.,
with $\widehat{W}=1$) is $0.509$.

###  2.2 Joint misclassification

Suppose that rather than observing $Y$ and $A$ we observe $Z$ and $B$, the
misclassified versions of $Y$ and $A$, respectively. The relation between $Z$
and $B$ on the one hand and $Y$, $A$ and $L$ on the other can be expressed as
follows:

$\displaystyle\Pr(Z=z,B=b|Y=y,A=a,L=l)$
$\displaystyle=(\pi_{byal})^{z}(1-\pi_{byal})^{1-z}(\lambda_{yal})^{b}(1-\lambda_{yal})^{1-b}$

for $z,b\in\\{0,1\\}$ and all possible realisations $y,a,l$ of $Y,A,L$, and
where $\pi_{byal}=\Pr(Z=1|B=b,Y=y,A=a,L=l)$ and
$\lambda_{yal}=\Pr(B=1|Y=y,A=a,L=l)$.

To simulate (dependent differential) misclassification in the reinfarction
dataset, we use the true positive and false positive rates given in Table 2.
The expected cell counts for these rates are given in Tables 3 and 4.

Table 2: True and false positive rates for reinfarction example. For $b,y,a,l\in\\{0,1\\}$, $\lambda_{yal}=\Pr(B=1|Y=y,A=a,L=l)$ and $\pi_{byal}=\Pr(Z=z|B=b,Y=y,A=a,L=l)$.  $\pi_{0000}=0.050$ | $\pi_{0001}=0.020$ | $\lambda_{000}=0.010$
---|---|---
$\pi_{1000}=0.060$ | $\pi_{1001}=0.108$ | $\lambda_{100}=0.181$
$\pi_{0100}=0.930$ | $\pi_{0101}=0.806$ | $\lambda_{010}=0.880$
$\pi_{1100}=0.938$ | $\pi_{1101}=0.692$ | $\lambda_{110}=0.910$
$\pi_{0010}=0.030$ | $\pi_{0011}=0.109$ | $\lambda_{001}=0.100$
$\pi_{1010}=0.060$ | $\pi_{1011}=0.050$ | $\lambda_{101}=0.265$
$\pi_{0110}=0.906$ | $\pi_{0111}=0.765$ | $\lambda_{011}=0.930$
$\pi_{1110}=0.950$ | $\pi_{1111}=0.861$ | $\lambda_{111}=0.823$
Table 3: Expected cell counts (rounded to integers) for reinfarction example after misclassification was introduced. Because of rounding, the sum of all cell entries is 33,006 rather than 33,007, the size of the reinfarction dataset. | $Z=0$ |  | $Z=1$
---|---|---|---
| $B=0$ | $B=1$ |  | $B=0$ | $B=1$
$Y=0,~{}A=0,~{}L=0$ | $10912$ | $109$ |  | $574$ | $7$
$Y=1,~{}A=0,~{}L=0$ | $51$ | $10$ |  | $678$ | $151$
$Y=0,~{}A=1,~{}L=0$ | $1527$ | $10850$ |  | $47$ | $693$
$Y=1,~{}A=1,~{}L=0$ | $5$ | $27$ |  | $48$ | $509$
$Y=0,~{}A=0,~{}L=1$ | $1148$ | $116$ |  | $23$ | $14$
$Y=1,~{}A=0,~{}L=1$ | $7$ | $4$ |  | $29$ | $9$
$Y=0,~{}A=1,~{}L=1$ | $334$ | $4738$ |  | $41$ | $249$
$Y=1,~{}A=1,~{}L=1$ | $4$ | $11$ |  | $13$ | $68$
Table 4: Expected cell counts (rounded to integers) for illustrative study setting after misclassification, collapsed over $Y$ and $A$. | $L=0$ |  | $L=1$
---|---|---|---
| $B=0$ | $B=1$ |  | $B=0$ | $B=1$
$Z=0$ | $12495$ | $10996$ |  | $1347$ | $1360$
$Z=1$ | $1493$ | $4869$ |  | $106$ | $340$

We redefine the weights in (2) as a function of $B$ and $L$ (according to
Appendix I) such that

$\displaystyle W$
$\displaystyle=\frac{p(B)\varepsilon_{BL}}{\sum_{y}\sum_{a}\pi_{ByaL}(\lambda_{yaL})^{B}(1-\lambda_{yaL})^{1-B}(\varepsilon_{aL})^{y}(1-\varepsilon_{aL})^{1-{y}}(\delta_{L})^{a}(1-\delta_{L})^{1-a}},$
(4)

where $p(B)$ is the prevalence of level $B$ of the potentially misclassified
version of the exposure variable and where $\varepsilon_{al}=\Pr(Y=1|A=a,L=l)$
and $\delta_{l}=\Pr(A=1|L=l)$ for all possible realisations $a$ and $l$ of $A$
and $L$, respectively. In Appendix I, it is shown that

OR
$\displaystyle=\frac{\mathbb{E}[WZ|B=1]/(1-\mathbb{E}[WZ|B=1])}{\mathbb{E}[WZ|B=0]/(1-\mathbb{E}[WZ|B=0])},$
(5)

which suggests the plug-in estimator

$\displaystyle\widehat{\text{OR}}$
$\displaystyle:=\frac{\widehat{\mathbb{E}}[\widehat{W}Z|B=1]/(1-\widehat{\mathbb{E}}[\widehat{W}Z|B=1])}{\widehat{\mathbb{E}}[\widehat{W}Z|B=0]/(1-\widehat{\mathbb{E}}[\widehat{W}Z|B=0])},$
(6)

where $\widehat{\mathbb{E}}$ denotes the sample mean operator and
$\widehat{W}$ the sample analogue (i.e., consistent estimator) of $W$ in (4).

In the absence of exposure misclassification, (4) reduces to

$\displaystyle
W=\Bigg{(}\frac{(\delta_{L})^{A}(1-\delta_{L})^{1-A}}{p(A)}\Bigg{[}\pi_{A0AL}\frac{1-\varepsilon_{AL}}{\varepsilon_{AL}}+\pi_{A1AL}\Bigg{]}\Bigg{)}^{-1}.$
(7)

The first term within the round brackets corrects for confounding and
represents the propensity score divided by prevalence of exposure level $A$.
The term within square brackets is a factor that corrects for
misclassification in the outcome variable. This correction factor is similar
to that proposed by Gravel and Platt (2018). The only difference is that where
in (7) it does not depend on the fallible measurement $Z$ of $Y$, Gravel and
Platt define different weights for subjects with $Z=0$. Note, however, that
the choice of weights for subjects with $Z=0$ does not affect the population
quantity in (5) or the estimator defined by (6).

As for the reinfarction example, the odds ratio estimate for the exposure-
outcome effect based on inverse probability weighting that assumes absence of
exposure or outcome misclassification is $1.120$, while the corresponding
misclassification naive crude odds ratio is $1.031$. Estimation of the
population weights $W$ from observables using validation data is discussed in
the next section. As shown below, weighting using the proposed weights that
account for confounding and outcome and exposure misclassification results in
an odds ratio of $\text{OR}=\widehat{\text{OR}}\approx 0.573$. Inference based
on (7) rather than (4), i.e., ignoring misclassification in the exposure but
correcting for outcome misclassification, yields an odds ratio estimate of
0.934.

###  2.3 Parameterisation based on positive and negative predictive values

In the foregoing discussion, the proposed weights were expressed in terms of
sensitivity and specificity parameters. The sensitivity and specificity of $Z$
with respect to $Y$, given $(B,A,L)$, are $\pi_{{B}1AL}$ and $1-\pi_{{B}0AL}$,
respectively. Similarly, $\lambda_{Y1L}$ and $1-\lambda_{Y0L}$ reflect the
sensitivity and specificity, respectively, with respect to $A$, conditional on
$Y$ and $L$.

As discussed below, it may be more convenient to choose a parameterisation
that is based on (positive and negative) predictive values. Define
$\delta^{\ast}_{l}=\Pr(B=1|L=l)$, $\varepsilon^{\ast}_{bl}=\Pr(Z=1|B=b,L=l)$,
$\lambda^{\ast}_{zbl}=\Pr(A=1|Z=z,B=b,L=l)$ and
$\pi^{\ast}_{azbl}=\Pr(Y=1|A=a,Z=z,B=b,L=l)$. The weights in (4) can be
rewritten as

$\displaystyle W$
$\displaystyle=\frac{\sum_{y}\sum_{a}\pi^{\ast}_{ByaL}(\lambda^{\ast}_{yaL})^{B}(1-\lambda^{\ast}_{yaL})^{1-B}(\varepsilon^{\ast}_{aL})^{y}(1-\varepsilon^{\ast}_{aL})^{1-{y}}(\delta^{\ast}_{L})^{a}(1-\delta^{\ast}_{L})^{1-a}}{\sum_{y}\sum_{a}(\lambda^{\ast}_{yaL})^{B}(1-\lambda^{\ast}_{yaL})^{1-B}(\varepsilon^{\ast}_{aL})^{y}(1-\varepsilon^{\ast}_{aL})^{1-{y}}(\delta^{\ast}_{L})^{a}(1-\delta^{\ast}_{L})^{1-a}}$
$\displaystyle\qquad\times\frac{p(B)}{\varepsilon^{\ast}_{BL}(\delta^{\ast}_{L})^{B}(1-\delta^{\ast}_{L})^{1-B}}.$
(8)

In the absence of exposure misclassification, these weights simplify to

$\displaystyle
W=\frac{p(A)}{(\delta_{L})^{A}(1-\delta_{L})^{1-A}}\frac{\varepsilon_{AL}}{\varepsilon^{\ast}_{AL}}.$

##  3 Estimation of weights

Estimation of the proposed weights can be done using a number of approaches
and we will here consider a maximum likelihood approach that assumes the
availability of internal validation data, i.e., that some study participants
have their observed exposure or outcome measured by an ‘infallible’ or ‘gold
standard’ (100% accurate) classifier.

###  3.1 Validation subset inclusion mechanism

Let $R_{Y}$ be the indicator variable that takes the value of 1 if the outcome
is observed (i.e., measured by an infallible classifier) and 0 otherwise.
Similarly, define $R_{A}$ to be the indicator variable that takes the value of
1 if the exposure variable is observed and 0 otherwise. $R_{Y}$ and $R_{A}$
reflect which subjects have validation data available on $Y$ and $A$,
respectively. The subset of subjects with validation data on $Y$ need not
fully overlap with the subset with validation data on $A$.

The validation subsets can be approached from the missing data framework of
Rubin (1976). Provided that $Z,B,L$ are free of missing values, Rubin’s
missing at random (MAR) condition is met whenever the vector $(R_{Y},R_{A})$
is conditionally independent of $(Y,A)$ given $(Z,B,L)$.

###  3.2 Full likelihood approach based on parameterisation in terms of
sensitivities and specificities

Simultaneous estimation of the whole vector of $\delta$, $\varepsilon$,
$\lambda$ and $\pi$ parameters can be done via maximum likelihood estimation
as follows. Assuming i.i.d. observations $(Z_{i},B_{i},Y_{i},A_{i},L_{i})$ and
ignorable missingness in the sense of Rubin (1976) (MAR and distinctness), for
valid likelihood-based inference it is appropriate to maximise the following
log-likelihood over the parameter space of $\theta$, the vector of $\delta$,
$\varepsilon$, $\lambda$ and $\pi$ parameters:

$\displaystyle\ell(\theta)$ $\displaystyle=\sum_{i:R_{Yi}=R_{Ai}=1}\log
f(\theta;Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle\qquad+\sum_{i:R_{Yi}=1\wedge
R_{Ai}=0}\log\sum_{A_{i}}f(\theta;Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle\qquad+\sum_{i:R_{Yi}=0\wedge
R_{Ai}=1}\log\sum_{Y_{i}}f(\theta;Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle\qquad+\sum_{i:R_{Yi}=R_{Ai}=0}\log\sum_{Y_{i}}\sum_{A_{i}}f(\theta;Z_{i},B_{i},Y_{i},A_{i},L_{i}),$

where

$\displaystyle f(\theta;Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle=(\pi_{B_{i}Y_{i}A_{i}L_{i}})^{Z_{i}}(1-\pi_{B_{i}Y_{i}A_{i}L_{i}})^{1-Z_{i}}(\lambda_{Y_{i}A_{i}L_{i}})^{B_{i}}(1-\lambda_{Y_{i}A_{i}L_{i}})^{1-B_{i}}$
$\displaystyle\qquad\times(\varepsilon_{A_{i}L_{i}})^{Y_{i}}(1-\varepsilon_{A_{i}L_{i}})^{1-{Y_{i}}}(\delta_{L_{i}})^{A_{i}}(1-\delta_{L_{i}})^{1-A_{i}}.$

Evaluating this log-likelihood involves marginalising over unobserved
quantities in the last three terms of $\ell(\theta)$. The log-likelihood
equations may become considerably more tractable if we choose a
parameterisation of the likelihood that is based on predictive values rather
than sensitivities and specificities.

###  3.3 Full likelihood approach based on parameterisation in terms of
predictive values

Inference may alternatively be based on a log-likelihood that is parameterised
in terms of the vector $\theta^{\ast}$ of the $\delta^{\ast}$,
$\varepsilon^{\ast}$, $\lambda^{\ast}$ and $\pi^{\ast}$ parameters, i.e.,

$\displaystyle\ell^{\ast}(\theta^{\ast})$
$\displaystyle=\sum_{i:R_{Yi}=R_{Ai}=1}\log
g(\theta^{\ast};Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle\qquad+\sum_{i:R_{Yi}=1\wedge
R_{Ai}=0}\log\sum_{A_{i}}g(\theta^{\ast};Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle\qquad+\sum_{i:R_{Yi}=0\wedge
R_{Ai}=1}\log\sum_{Y_{i}}g(\theta^{\ast};Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle\qquad+\sum_{i:R_{Yi}=R_{Ai}=0}\log\sum_{Y_{i}}\sum_{A_{i}}g(\theta^{\ast};Z_{i},B_{i},Y_{i},A_{i},L_{i}),$

where

$\displaystyle g(\theta^{\ast};Z_{i},B_{i},Y_{i},A_{i},L_{i})$
$\displaystyle=(\pi^{\ast}_{A_{i}Z_{i}B_{i}L_{i}})^{Y_{i}}(1-\pi^{\ast}_{A_{i}Z_{i}B_{i}L_{i}})^{1-Y_{i}}(\lambda^{\ast}_{Z_{i}B_{i}L_{i}})^{A_{i}}(1-\lambda^{\ast}_{Z_{i}B_{i}L_{i}})^{1-A_{i}}$
$\displaystyle\qquad\times(\varepsilon^{\ast}_{B_{i}L_{i}})^{Z_{i}}(1-\varepsilon^{\ast}_{B_{i}L_{i}})^{1-{Z_{i}}}(\delta^{\ast}_{L_{i}})^{B_{i}}(1-\delta^{\ast}_{L_{i}})^{1-B_{i}}.$

If validation data is available on $Y$ if and only if it is available on $A$,
the complete data log-likelihood ignoring the missing data mechanism can be
conveniently expressed as follows:

$\displaystyle\ell^{\ast}(\theta^{\ast})$
$\displaystyle=\ell_{1}^{\ast}(\theta^{\ast})+\ell_{2}^{\ast}(\theta^{\ast})+\ell_{3}^{\ast}(\theta^{\ast})+\ell_{4}^{\ast}(\theta^{\ast}),$
(9)

with $\theta^{\ast}$ denoting the vector of $\delta^{\ast}$,
$\varepsilon^{\ast}$, $\lambda^{\ast}$ and $\pi^{\ast}$ parameters and where

$\displaystyle\ell_{1}^{\ast}(\theta^{\ast})$
$\displaystyle=\sum_{i:R_{Yi}=R_{Ai}=1}Y_{i}\log(\pi^{\ast}_{A_{i}Z_{i}B_{i}L_{i}})+(1-Y_{i})\log(1-\pi^{\ast}_{A_{i}Z_{i}B_{i}L_{i}})$
$\displaystyle\ell_{2}^{\ast}(\theta^{\ast})$
$\displaystyle=\sum_{i:R_{Yi}=R_{Ai}=1}A_{i}\log(\lambda^{\ast}_{Z_{i}B_{i}L_{i}})+(1-A_{i})\log(1-\lambda^{\ast}_{Z_{i}B_{i}L_{i}})$
$\displaystyle\ell_{3}^{\ast}(\theta^{\ast})$
$\displaystyle=\sum_{i}Z_{i}\log(\varepsilon^{\ast}_{B_{i}L_{i}})+(1-Z_{i})\log(1-\varepsilon^{\ast}_{B_{i}L_{i}})$
$\displaystyle\ell_{4}^{\ast}(\theta^{\ast})$
$\displaystyle=\sum_{i}B_{i}\log(\delta^{\ast}_{L_{i}})+(1-B_{i})\log(1-\delta^{\ast}_{L_{i}}).$

Now, assuming distinct parameter spaces for the vectors of $\pi^{\ast}$,
$\lambda^{\ast}$, $\varepsilon^{\ast}$, and $\delta^{\ast}$ parameters, the
parameter values that maximise $\ell^{\ast}(\theta^{\ast})$ can be found by
separately maximising $\ell_{1}^{\ast}(\theta^{\ast})$ and
$\ell_{2}^{\ast}(\theta^{\ast})$ in the validation subset with respect to the
$\pi^{\ast}$ and $\lambda^{\ast}$ parameters, respectively, and
$\ell_{3}^{\ast}(\theta^{\ast})$ and $\ell_{4}^{\ast}(\theta^{\ast})$ in the
entire dataset with respect to $\varepsilon^{\ast}$ and $\delta^{\ast}$.
Following Gravel and Platt (2018) and Tang et al. (2013), the sum of the first
and last two terms are therefore suitably labelled the internal validation and
main study log-likelihood, respectively. With this parameterisation, finding
the maximum likelihood estimates is readily achieved by taking advantage of
standard statistical software.

###  3.4 Equivalence of likelihood approaches based on different
parameterisations

Without restrictions imposed on

$\displaystyle\theta_{l}$
$\displaystyle:=(\pi_{000l},\pi_{100l},\pi_{010l},\pi_{110l},\pi_{001l},\pi_{101l},\pi_{011l},\pi_{111l},\lambda_{00l},\lambda_{10l},\lambda_{01l},\lambda_{11l},\varepsilon_{0l},\varepsilon_{1l},\delta_{l})\text{~{}~{}and}$
$\displaystyle\theta^{\ast}_{l}$
$\displaystyle:=(\pi^{\ast}_{000l},\pi^{\ast}_{100l},\pi^{\ast}_{010l},\pi^{\ast}_{110l},\pi^{\ast}_{001l},\pi^{\ast}_{101l},\pi^{\ast}_{011l},\pi^{\ast}_{111l},\lambda^{\ast}_{00l},\lambda^{\ast}_{10l},\lambda^{\ast}_{01l},\lambda^{\ast}_{11l},\varepsilon^{\ast}_{0l},\varepsilon^{\ast}_{1l},\delta^{\ast}_{l}),$

other than that $\theta_{l},\theta_{l}^{\ast}\in(0,1)^{15}$, it can be shown
that the maximum likelihood estimator based on the internal validation design
is invariant to its parameterisation (sensitivities/specificities versus
positive and negative predictive values). This is because there exists a
function mapping every $\theta_{l}\in(0,1)^{15}$ to a unique
$\theta^{\ast}_{l}\in(0,1)^{15}$ and vice versa. Maximising $\ell(\theta)$
with respect to $\theta$ is then equivalent to maximising
$\ell(\sigma(\theta^{\ast}))$ ($=\ell^{\ast}(\theta^{\ast})$) with respect to
$\theta^{\ast}$ for some bijection $\sigma$ such that
$\theta=\sigma(\theta^{\ast})$; that is,

$\displaystyle\underset{\theta}{\arg\max}~{}\ell(\theta)=\sigma\left(\underset{\theta^{\ast}}{\arg\max}~{}\ell(\sigma(\theta^{\ast}))\right).$

If more restrictions are imposed on $\theta$ or $\theta$, e.g., if we assume
non-saturated logistic models for the components of $\theta$ and
$\theta^{\ast}$, this equivalence no longer holds and the resulting weight
estimates may differ depending on the parameterisation.

###  3.5 Application

For the re-infarction data example, we assume validation data are available
according to a MAR mechanism characterised by

$\displaystyle\Pr(R_{Y}=1|R_{A}=s,Z=z,B=b,Y=y,A=a,L=l)$ $\displaystyle=s,$
$\displaystyle\Pr(R_{A}=1|Z=z,B=b,Y=y,A=a,L=l)$ $\displaystyle=0.25+0.10b.$

This mechanism assigns validation data to an individual on either both $Y$ and
$A$ (30% of all individuals) or neither depending on their realisation of $B$,
the misclassified version of the exposure variable $A$ (Table 5). Tables S.1
and S.2 (see Supplementary Web Appendix) give the likelihood contributions for
the parameterisation based on predictive values and the closed form maximum
likelihood expressions, respectively. Maximum likelihood estimates can also be
found by fitting to the data the saturated logistic regression models of $B$
and $Z$ on $L$ and $(B,L)$, respectively, and to the validation subset the
fully saturated logistic regression models of $A$ and $Y$ on $(Z,B,L)$ and
$(A,Z,B,L)$, respectively. Estimated weights are then obtained by plugging in
the maximum likelihood estimates into (8). As in the complete data setting
where we assumed the weights to be known, evaluating (6) then yields an odds
ratio of $\widehat{\text{OR}}=\text{OR}\approx 0.573$.

Table 5: Expected cell counts (rounded to integers) for illustrative study setting after misclassification and formation of validation subsets. |  |  |  |  | $B=0$ |  | $B=1$
---|---|---|---|---|---|---|---
$R_{Y}$ | $R_{A}$ | $Y$ | $A$ | $L$ | $Z=0$ | $Z=1$ |  | $Z=0$ | $Z=1$
$0$ | $0$ |  |  | $0$ | $m_{1}=9371$ | $m_{2}=7147$ |  | $m_{3}=1011$ | $m_{4}=\phantom{0}884$
$0$ | $0$ |  |  | $1$ | $m_{5}=1120$ | $m_{6}=3165$ |  | $m_{7}=\phantom{0}\phantom{0}80$ | $m_{8}=\phantom{0}221$
$0$ | $1$ |  |  |  | $m_{9}=\phantom{0}\phantom{0}\phantom{0}0$ | $m_{10}=\phantom{0}\phantom{0}\phantom{0}0$ |  | $m_{11}=\phantom{0}\phantom{0}\phantom{0}0$ | $m_{12}=\phantom{0}\phantom{0}\phantom{0}0$
$1$ | $0$ |  |  |  | $m_{13}=\phantom{0}\phantom{0}\phantom{0}0$ | $m_{14}=\phantom{0}\phantom{0}\phantom{0}0$ |  | $m_{15}=\phantom{0}\phantom{0}\phantom{0}0$ | $m_{16}=\phantom{0}\phantom{0}\phantom{0}0$
$1$ | $1$ | $0$ | $0$ | $0$ | $m_{17}=2728$ | $m_{18}=\phantom{0}\phantom{0}38$ |  | $m_{19}=\phantom{0}144$ | $m_{20}=\phantom{0}\phantom{0}\phantom{0}2$
$1$ | $1$ | $1$ | $0$ | $0$ | $m_{21}=\phantom{0}\phantom{0}13$ | $m_{22}=\phantom{0}\phantom{0}\phantom{0}3$ |  | $m_{23}=\phantom{0}169$ | $m_{24}=\phantom{0}\phantom{0}53$
$1$ | $1$ | $0$ | $1$ | $0$ | $m_{25}=\phantom{0}382$ | $m_{26}=3797$ |  | $m_{27}=\phantom{0}\phantom{0}12$ | $m_{28}=\phantom{0}242$
$1$ | $1$ | $1$ | $1$ | $0$ | $m_{29}=\phantom{0}\phantom{0}\phantom{0}1$ | $m_{30}=\phantom{0}\phantom{0}\phantom{0}9$ |  | $m_{31}=\phantom{0}\phantom{0}12$ | $m_{32}=\phantom{0}178$
$1$ | $1$ | $0$ | $0$ | $1$ | $m_{33}=\phantom{0}287$ | $m_{34}=\phantom{0}\phantom{0}41$ |  | $m_{35}=\phantom{0}\phantom{0}\phantom{0}6$ | $m_{36}=\phantom{0}\phantom{0}\phantom{0}5$
$1$ | $1$ | $1$ | $0$ | $1$ | $m_{37}=\phantom{0}\phantom{0}\phantom{0}2$ | $m_{38}=\phantom{0}\phantom{0}\phantom{0}1$ |  | $m_{39}=\phantom{0}\phantom{0}\phantom{0}7$ | $m_{40}=\phantom{0}\phantom{0}\phantom{0}3$
$1$ | $1$ | $0$ | $1$ | $1$ | $m_{41}=\phantom{0}\phantom{0}84$ | $m_{42}=1658$ |  | $m_{43}=\phantom{0}\phantom{0}10$ | $m_{44}=\phantom{0}\phantom{0}87$
$1$ | $1$ | $1$ | $1$ | $1$ | $m_{45}=\phantom{0}\phantom{0}\phantom{0}1$ | $m_{46}=\phantom{0}\phantom{0}\phantom{0}4$ |  | $m_{47}=\phantom{0}\phantom{0}\phantom{0}3$ | $m_{48}=\phantom{0}\phantom{0}24$

##  4 Simulations

We performed a series of Monte Carlo simulation experiments to illustrate the
implementation of the proposed method, to study its finite sample properties
and to compare the method to estimators that ignore the presence of
confounding or joint exposure and outcome misclassification. All simulations
were conducted using R-3.5.0 (R Core Team, 2018) on x86_64-pc-linux-gnu
platforms of the high performance computer cluster of Leiden University
Medical Center.

###  4.1 Methods

For all 36 simulation experiments, we generated $n_{\mathrm{sim}}=1000$
samples of size $n$ according to the data generating mechanisms depicted in
the directed acyclic graphs of Figure 1. This multi-step data generating
process included generating values on measurement error-free variables,
introducing misclassification and allocating individuals validation data. We
applied various estimators to each of the simulation samples to yield, for
each scenario, an empirical distribution of each point estimator and
corresponding precision estimators. These distributions were then summarised
into various performance metrics. These metrics include the empirical bias of
the estimator on the log-scale (i.e., the mean estimated log-OR minus the
target log-OR across the $n_{\mathrm{sim}}$ samples), the empirical standard
error (SE) of the estimator on the log-scale (i.e., the square root of the
mean squared deviation of the estimated log-OR from the mean log-OR), the
empirical mean squared error (MSE) (i.e., the sum of the squared SE and the
squared bias), the square root of the mean estimated variance (SSE, sample
standard error) and the empirical coverage probability (CP) (i.e., the
fraction of simulation runs per scenario where the 95% confidence interval
(95%CI) contained the target quantity).

####  4.1.1 Distribution of measurement error-free variables

Following Gravel and Platt (2018), we consider a setting based on that of
“Scenario A” in Setoguchi et al. (2008) with slight modifications to the
propensity score and outcome models. We consider a fully observed covariate
vector $L=(L_{0},...,L_{10})$ whose distribution coincides with that of
$h(V)$, where $V=(V_{1},...,V_{10})$ has the multivariate normal distribution
with zero means, unit variances and correlations equal to zero except for the
correlations between $W_{1}$ and $V_{5}$, $V_{2}$ and $V_{6}$, $V_{3}$ and
$V_{8}$, and $V_{4}$ and $V_{9}$, which were set to $0.2$, $0.9$, $0.2$, and
$0.9$, respectively. Function $h$ was defined such that

$h(V)=(I(V_{1}>0),V_{2},I(V_{3}>0),V_{4},I(V_{5}>0),I(V_{6}>0),V_{7},I(V_{8}>0),I(V_{9}>0),V_{10}).$

Thus, sampling from the distribution of $L$ is equivalent to sampling from the
multivariate normal distribution with the given parameter values and
dichotomising the 1st, 3rd, 5th, 6th, 8th and 9th elements.

Next, let $U_{1}$ and $U_{2}$ be binary variables distributed according to the
following logistic models:

$\displaystyle\mathrm{logit}\Pr(U_{1}=1|L)$ $\displaystyle=\eta_{0},$ (10)
$\displaystyle\mathrm{logit}\Pr(U_{2}=1|L,U_{1})$ $\displaystyle=\mu_{0}.$
(11)

The distribution of the binary exposure variable $A$ was defined according to
the model

$\displaystyle\mathrm{logit}\Pr(A=1|L,U_{1},U_{2})=\alpha_{0}+\textstyle\sum_{j=1}^{10}\alpha_{j}L_{j}+\alpha_{11}U_{1}.$
(12)

Letting $U_{3}$ be a scalar random variable that is independent of
$(A,L_{1},...,L_{10},U_{1},U_{2})$ and uniformly distributed over the interval
$[0,1]$, we defined the counterfactual outcome $Y(a)$, under the intervention
setting $A$ to $a$, as

$\displaystyle Y(a)=I\Big{(}U_{3}<\mathrm{expit}\Big{\\{}\beta_{0}+\gamma
a+\textstyle\sum_{j=1}^{10}\beta_{j}L_{j}+\beta_{11}U_{2}\Big{\\}}\Big{)}.$
(13)

With $Y:=Y(A)$, the above implies consistency, conditional exchangeability
given $L$ and structural positivity.

####  4.1.2 Misclassification mechanism

For scenarios with joint misclassification, we defined $B=U_{1}$ and
$Z=U_{2}$, so that the predictive values take a standard logistic form:

$\displaystyle\mathrm{logit}\Pr(Y=1|A,B,L,Z)$ $\displaystyle=\beta_{0}+\gamma
A+\textstyle\sum_{j=1}^{10}\beta_{j}L_{j}+\beta_{11}Z$ (14)
$\displaystyle\mathrm{logit}\Pr(A=1|B,L,Z)$
$\displaystyle=\alpha_{0}+\textstyle\sum_{j=1}^{10}\alpha_{j}L_{j}+\alpha_{11}B.$
(15)

For scenarios without exposure misclassification, we set $\alpha_{11}=0$ and
defined $B=A$ and $Z=U_{2}$, so that

$\displaystyle\mathrm{logit}\Pr(Y=1|A,B,L,Z)$ $\displaystyle=\beta_{0}+\gamma
A+\textstyle\sum_{j=1}^{10}\beta_{j}L_{j}+\beta_{11}Z$ (16)
$\displaystyle\mathrm{logit}\Pr(B=1|L,Z)$
$\displaystyle=\alpha_{0}+\textstyle\sum_{j=1}^{10}\alpha_{j}L_{j}.$ (17)

For simplicity, we removed any marginal dependence of $Z$ on the covariates
$L$ and $U_{1}$ as well as any marginal dependence of $U_{1}$ on $L$ (cf.
equations (10) and (11)). Although models (10) through (15) take a standard
logistic form, they do not imply that the corresponding sensitivities and
specificities can be written in the same form. We chose the predictive values
rather than the sensitivities and specificities to take a standard logistic
form so as to ensure correct model specification in the estimation of the
weights in the simulation experiments, in which a likelihood approach based on
predictive values was adopted (cf. (9)).

Figure 1: Data structure for scenarios with misclassification on the outcome
only (left) or on both the exposure and outcome (right). Bullet arrowheads
represent deterministic relationships.

####  4.1.3 Missing data mechanism

For these simulations, we stipulated $L$, $B$ and $Z$ to be observed for all
subjects. We consider scenarios where the dataset can be partitioned into a
subset with validation data on all misclassified variables (denoted $R=1$) and
a dataset with validation data on neither ($R=0$). That is, we simulated data
such that subjects have validation data on both $A$ and $Y$ or neither on $A$
nor on $Y$. Values for the response indicator $R$ were generated according to
the following (MAR) model:

$\displaystyle\mathrm{logit}\Pr(R=1|Z,B,Y,A,L)$
$\displaystyle=\mathrm{logit}\Pr(R=1|Z,B,L)$
$\displaystyle=\xi_{0}+\xi_{1}Z+\xi_{2}B+\xi_{3}ZB.$

####  4.1.4 Scenarios

We initially fixed most parameters of models (12) and (13) at the respective
values of “Scenario A” of Setoguchi et al. (2008): $\alpha_{1}=0.8$,
$\alpha_{2}=-0.25$, $\alpha_{3}=0.6$, $\alpha_{4}=-0.4$, $\alpha_{5}=-0.8$,
$\alpha_{6}=-0.5$, $\alpha_{7}=0.7$, $\alpha_{8}=0$, $\alpha_{9}=0$,
$\alpha_{10}=0$, $\beta_{0}=-3.85$, $\beta_{1}=0.3$, $\beta_{2}=-0.36$,
$\beta_{3}=-0.73$, $\beta_{4}=-0.2$, $\beta_{5}=0$, $\beta_{6}=0$,
$\beta_{7}=0$, $\beta_{8}=0.71$, $\beta_{9}=-0.19$ and $\beta_{10}=0.26$.
Parameters $\eta_{0}$ and $\alpha_{0}$ were fixed at zero and $\xi_{1}$,
$\xi_{2}$ and $\xi_{3}$ at $2$, $1$ and $-1$, respectively. The remaining
parameters and $\beta_{0}$ were allowed to vary across scenarios as per Table
6.

Scenarios differ by sample size $n$, the presence of outcome
misclassification, potentially misclassified outcome prevalence (via
$\mu_{0}$), the associations between the exposure and outcome on the one hand
and the respective misclassified versions on the other (via $\alpha_{11}$ and
$\beta_{11}$), outcome model intercept $\beta_{0}$, the conditional log-OR
$\gamma$, or the size of the validation subset (via $\xi_{0}$). Based on an
iterative Monte Carlo integration approach (Austin and Stafford, 2008), we
specified $\gamma$ so as to keep the target marginal log odds ratio at $-0.4$.

Table 6: Simulation parameter values used in the Monte Carlo studies. | Exposure |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---
Scenario | misclassification | $n$ | $\mu_{0}$ | $\alpha_{11}$ | $\beta_{0}$ | $\beta_{11}$ | $\gamma$ | $\xi_{0}$
1 | Absent | $5000$ | $-2$ | $0$ | $-3.85$ | $2$ | $-0.431$ | $-1.5$
2 | Absent | $5000$ | $-3$ | $0$ | $-3.85$ | $2$ | $-0.417$ | $-1.5$
3 | Absent | $5000$ | $-2$ | $0$ | $-3.85$ | $4$ | $-0.624$ | $-1.5$
4 | Absent | $5000$ | $-2$ | $0$ | $-3.85$ | $2$ | $-0.431$ | $-2.5$
5 | Present | $5000$ | $-2$ | $2$ | $-3.85$ | $2$ | $-0.431$ | $-1.5$
6 | Present | $5000$ | $-3$ | $2$ | $-3.85$ | $2$ | $-0.417$ | $-1.5$
7 | Present | $5000$ | $-2$ | $4$ | $-3.85$ | $2$ | $-0.431$ | $-1.5$
8 | Present | $5000$ | $-2$ | $2$ | $-3.85$ | $4$ | $-0.624$ | $-1.5$
9 | Present | $5000$ | $-2$ | $2$ | $-3.85$ | $2$ | $-0.431$ | $-2.5$
10 | Absent | $10000$ | $-2$ | $0$ | $-3.85$ | $2$ | $-0.431$ | $-1.5$
11 | Absent | $10000$ | $-3$ | $0$ | $-3.85$ | $2$ | $-0.417$ | $-1.5$
12 | Absent | $10000$ | $-2$ | $0$ | $-3.85$ | $4$ | $-0.624$ | $-1.5$
13 | Absent | $10000$ | $-2$ | $0$ | $-3.85$ | $2$ | $-0.431$ | $-2.5$
14 | Present | $10000$ | $-2$ | $2$ | $-3.85$ | $2$ | $-0.431$ | $-1.5$
15 | Present | $10000$ | $-3$ | $2$ | $-3.85$ | $2$ | $-0.417$ | $-1.5$
16 | Present | $10000$ | $-2$ | $4$ | $-3.85$ | $2$ | $-0.431$ | $-1.5$
17 | Present | $10000$ | $-2$ | $2$ | $-3.85$ | $4$ | $-0.624$ | $-1.5$
18 | Present | $10000$ | $-2$ | $2$ | $-3.85$ | $2$ | $-0.431$ | $-2.5$
19 | Absent | $5000$ | $-2$ | $0$ | $-2$ | $2$ | $-0.470$ | $-1.5$
20 | Absent | $5000$ | $-3$ | $0$ | $-2$ | $2$ | $-0.445$ | $-1.5$
21 | Absent | $5000$ | $-2$ | $0$ | $-2$ | $4$ | $-0.641$ | $-1.5$
22 | Absent | $5000$ | $-2$ | $0$ | $-2$ | $2$ | $-0.470$ | $-2.5$
23 | Present | $5000$ | $-2$ | $2$ | $-2$ | $2$ | $-0.470$ | $-1.5$
24 | Present | $5000$ | $-3$ | $2$ | $-2$ | $2$ | $-0.445$ | $-1.5$
25 | Present | $5000$ | $-2$ | $4$ | $-2$ | $2$ | $-0.470$ | $-1.5$
26 | Present | $5000$ | $-2$ | $2$ | $-2$ | $4$ | $-0.641$ | $-1.5$
27 | Present | $5000$ | $-2$ | $2$ | $-2$ | $2$ | $-0.470$ | $-2.5$
28 | Absent | $10000$ | $-2$ | $0$ | $-2$ | $2$ | $-0.470$ | $-1.5$
29 | Absent | $10000$ | $-3$ | $0$ | $-2$ | $2$ | $-0.445$ | $-1.5$
30 | Absent | $10000$ | $-2$ | $0$ | $-2$ | $4$ | $-0.641$ | $-1.5$
31 | Absent | $10000$ | $-2$ | $0$ | $-2$ | $2$ | $-0.470$ | $-2.5$
32 | Present | $10000$ | $-2$ | $2$ | $-2$ | $2$ | $-0.470$ | $-1.5$
33 | Present | $10000$ | $-3$ | $2$ | $-2$ | $2$ | $-0.445$ | $-1.5$
34 | Present | $10000$ | $-2$ | $4$ | $-2$ | $2$ | $-0.470$ | $-1.5$
35 | Present | $10000$ | $-2$ | $2$ | $-2$ | $4$ | $-0.641$ | $-1.5$
36 | Present | $10000$ | $-2$ | $2$ | $-2$ | $2$ | $-0.470$ | $-2.5$

####  4.1.5 Estimators

We considered five estimators of the OR for the marginal exposure-outcome
effect: a crude estimator (labeled Crude) that ignores both confounding and
misclassication of any variable, a misclassification naive estimator (labeled
PS) that addresses confounding through IPW, complete cases analysis (CCA) in
which IPW is applied only to the subset of subjects with validation data, the
Gravel and Platt estimator (GP) that ignores exposure misclassification, and
the method proposed in this article (labeled IPWM). Both GP and IPWM are
implemented using the R function mecor::ipwm (Nab, 2019; Nab et al., 2018),
which in the simulation settings considered uses iteratively reweighted least
squares via the stats::glm function for maximum likelihood estimation. Unlike
Gravel and Platt (Gravel and Platt, 2018), we used a non-parametric rather
than a semi-parametric bootstrap procedure for estimating standard errors and
constructing confidence intervals. Semi-parametrically generating response
indicators would preferably require modelling of (or making additional
assumptions about) the missing data mechanism. For all methods and each
original dataset, we drew 1000 bootstrap samples for variance estimation and
the construction of percentile confidence intervals.

All estimators are based on a function of the estimated outcome probability
$P_{1}$ in the exposed group and the estimated outcome probability $P_{0}$ in
the unexposed group. However, since $P_{1}$ and $P_{0}$ may take a value of 0
or 1, the crude odds ratio $[P_{1}/(1-P_{1})]/[P_{0}/(1-P_{0})]$ need not
exist. In contrast to what is often (implicitly) done in simulation
studies—i.e., studying the properties of the estimators after conditioning on
datasets where $[P_{1}/(1-P_{1})]/[P_{0}/(1-P_{0})]$ is defined—we first
define $P_{1}^{\ast}=(P_{1}s+1)/(s+2)$ and $P_{0}^{\ast}=(P_{0}s+1)/(s+2)$ for
a large positive number $s$ (here set to $10^{6}$) and then regard
$[P_{1}^{\ast}/(1-P_{1}^{\ast})]/[P_{0}^{\ast}/(1-P_{0}^{\ast})]$ as the
estimator of the OR for the exposure-outcome association. This ensures the
estimator is always defined and effectively shrinks the outcome probabilities
towards 0.5 and the OR towards 1 (Appendix II).

For PS and CCA, we used a logistic regression of $B$ and $A$, respectively, on
covariates $L_{1}$ through $L_{10}$ as main effects to estimate the propensity
scores. Taking the crude OR for the association between $B$ and $Z$ (PS) or
$A$ and $Y$ (CCA) over the data weighted by the reciprocal of the propensity
scores provided an estimate of target OR. R code for the methods GP and IPWM
is given in Appendix III.

###  4.2 Results

The treatment assignment mechanism detailed above resulted in average exposure
rates ranging from 17% to 51%, whereas average outcome rates ranged from 3% to
22%. Across all simulation studies, the average outcome, exposure and joint
misclassification rates ranged from 6% to 18%, from 0 to 33% and from 0% to
6%, respectively. Approximately 16% to 32% of subjects were allocated
validation data.

The results on the performance of the various methods in simulations studies
1-9 are provided in Table 7 (see Supplementary Table S.3 for the results on
all scenarios).

As expected, Crude, PS and CCA clearly showed bias with respect to the target
log OR of $-0.4$. The bias associated with restricting the analysis to records
with validation data is likely brought on to a large extent by collider
stratification, with $R$ acting as the collider here (cf. Figure 1). Both
Crude and PS indicated a null effect, as one would anticipate in view of the
marginal and $L$-conditional independence of $B$ and $Z$ implied by the
simulation set-up. The empirical coverage probabilities were, although low for
both estimators, similar to substantially larger for PS as compared with
Crude. Paralleling this is that Crude, whose (implicit) propensity score model
is inherently at least as parsimonious, yielded similar to smaller empirical
and sample standard errors as compared with PS. With the average fraction of
subjects with validation data being as low as 16% (in scenarios with low
$\xi_{0}$) to 32%, it is not unsurprising that Crude was subject to the
largest degree of variability.

The results for the IPWM approach are generally favourable and in line with
its theoretical (large sample) properties. Note that the results for GP and
IPWM are identical for scenarios 1-4, 10-13, 19-22 and 28-31, since the
methods are equivalent in the absence of exposure misclassification. In all
other scenarios, i.e., scenarios for which GP was not developed, GP performed
substantially worse than IPWM. The non-zero, albeit relatively small,
systematic deviations of the IPWM point estimates from the target $-0.4$,
notably the estimated bias of $-0.097$ (scenario 2), may be attributable in
part to the outcome being rare (with prevalence ranging from 3% to 8% across
scenarios 1-9). This is indicated by the superior performance of IPWM in
scenarios where the outcome is more prevalent (scenarios 10-36, prevalence up
to 22%). A similar observation was made by Gravel and Platt (2018). The
standard errors for GP and IPWM were noticeably higher than those of Crude and
PS, which is unsurprising in view of the discrepancies in the number of
estimated parameters. As expected, increasing the sample size, the true
outcome rate (via $\beta_{0}$) or both led to a decrease in the variability of
IPWM (cf. Table 6 and Supplementary Table S.3). Throughout the empirical
coverage probabilities of IPWM were close to the nominal level of 0.95.

Table 7: Results for simulation studies 1-9 on the performance of different causal estimators in various scenarios of confounding and misclassification in exposure and outcome. Abbreviations: PS, propensity score method ignoring misclassification; CCA, complete case analysis; GP, Gravel and Platt estimator ignoring exposure misclassification; IPWM, inverse probability weighting method for confounding and joint exposure and outcome misclassification; BSE, estimated standard error for the bias due to Monte Carlo error; SE, empirical standard error; SSE, sample standard error; CP, empirical coverage probability. In all scenarios, the true marginal log OR (estimand) was $-0.4$. | Crude
---|---
Scenario | Bias | BSE | MSE | SE | SSE | CP
1 | $\phantom{-}0.394$ | $0.004$ | $0.275$ | $0.119$ | $0.118$ | $0.122$
2 | $\phantom{-}0.382$ | $0.006$ | $0.329$ | $0.183$ | $0.184$ | $0.492$
3 | $\phantom{-}0.394$ | $0.004$ | $0.272$ | $0.117$ | $0.118$ | $0.116$
4 | $\phantom{-}0.401$ | $0.004$ | $0.278$ | $0.117$ | $0.118$ | $0.102$
5 | $\phantom{-}0.401$ | $0.003$ | $0.252$ | $0.090$ | $0.088$ | $0.007$
6 | $\phantom{-}0.407$ | $0.004$ | $0.298$ | $0.132$ | $0.134$ | $0.133$
7 | $\phantom{-}0.396$ | $0.003$ | $0.243$ | $0.086$ | $0.088$ | $0.009$
8 | $\phantom{-}0.395$ | $0.003$ | $0.242$ | $0.086$ | $0.088$ | $0.005$
9 | $\phantom{-}0.398$ | $0.003$ | $0.247$ | $0.089$ | $0.088$ | $0.005$
| PS
Scenario | Bias | BSE | MSE | SE | SSE | CP
1 | $\phantom{-}0.392$ | $0.005$ | $0.321$ | $0.168$ | $0.169$ | $0.382$
2 | $\phantom{-}0.379$ | $0.008$ | $0.407$ | $0.264$ | $0.258$ | $0.738$
3 | $\phantom{-}0.389$ | $0.006$ | $0.327$ | $0.175$ | $0.169$ | $0.402$
4 | $\phantom{-}0.389$ | $0.006$ | $0.327$ | $0.176$ | $0.168$ | $0.392$
5 | $\phantom{-}0.402$ | $0.003$ | $0.252$ | $0.090$ | $0.088$ | $0.010$
6 | $\phantom{-}0.407$ | $0.004$ | $0.297$ | $0.131$ | $0.135$ | $0.136$
7 | $\phantom{-}0.396$ | $0.003$ | $0.243$ | $0.086$ | $0.088$ | $0.009$
8 | $\phantom{-}0.395$ | $0.003$ | $0.242$ | $0.086$ | $0.088$ | $0.004$
9 | $\phantom{-}0.398$ | $0.003$ | $0.247$ | $0.089$ | $0.088$ | $0.005$
| CCA
Scenario | Bias | BSE | MSE | SE | SSE | CP
1 | $-0.078$ | $0.015$ | $0.476$ | $0.469$ | $0.491$ | $0.899$
2 | $-0.117$ | $0.019$ | $0.615$ | $0.601$ | $0.900$ | $0.887$
3 | $-0.020$ | $0.010$ | $0.301$ | $0.301$ | $0.300$ | $0.919$
4 | $-0.093$ | $0.020$ | $0.640$ | $0.631$ | $1.158$ | $0.899$
5 | $-0.145$ | $0.009$ | $0.307$ | $0.286$ | $0.286$ | $0.903$
6 | $-0.109$ | $0.011$ | $0.357$ | $0.345$ | $0.362$ | $0.930$
7 | $-0.213$ | $0.007$ | $0.282$ | $0.237$ | $0.250$ | $0.865$
8 | $-0.209$ | $0.006$ | $0.231$ | $0.187$ | $0.186$ | $0.775$
9 | $-0.175$ | $0.012$ | $0.422$ | $0.392$ | $0.411$ | $0.902$

(Continued.)  GP Scenario Bias BSE MSE SE SSE CP 1 $-0.036$ $0.011$ $0.360$
$0.359$ $0.428$ $0.958$ 2 $-0.097$ $0.016$ $0.515$ $0.505$ $0.861$ $0.938$ 3
$-0.019$ $0.007$ $0.234$ $0.233$ $0.240$ $0.939$ 4 $-0.045$ $0.016$ $0.503$
$0.501$ $1.087$ $0.944$ 5 $\phantom{-}0.269$ $0.008$ $0.316$ $0.244$ $0.244$
$0.799$ 6 $\phantom{-}0.280$ $0.010$ $0.392$ $0.314$ $0.339$ $0.862$ 7
$\phantom{-}0.134$ $0.008$ $0.259$ $0.241$ $0.252$ $0.926$ 8
$\phantom{-}0.259$ $0.004$ $0.207$ $0.140$ $0.144$ $0.570$ 9
$\phantom{-}0.263$ $0.010$ $0.394$ $0.325$ $0.339$ $0.883$ IPWM Scenario Bias
BSE MSE SE SSE CP 1 $-0.036$ $0.011$ $0.360$ $0.359$ $0.428$ $0.958$ 2
$-0.097$ $0.016$ $0.515$ $0.505$ $0.861$ $0.938$ 3 $-0.019$ $0.007$ $0.234$
$0.233$ $0.240$ $0.939$ 4 $-0.045$ $0.016$ $0.503$ $0.501$ $1.087$ $0.944$ 5
$-0.017$ $0.009$ $0.286$ $0.286$ $0.284$ $0.942$ 6 $-0.014$ $0.011$ $0.359$
$0.359$ $0.386$ $0.958$ 7 $\phantom{-}0.004$ $0.008$ $0.243$ $0.243$ $0.261$
$0.969$ 8 $-0.004$ $0.006$ $0.180$ $0.180$ $0.181$ $0.958$ 9 $-0.025$ $0.012$
$0.375$ $0.374$ $0.415$ $0.956$

##  5 Discussion

The analysis of epidemiologic data is often complicated by the presence of
confounding and misclassifications in exposure and outcome variables. In this
paper we proposed a new estimator for estimating a marginal odds-ratio in the
presence of confouding and joint misclassification of the exposure and outcome
variables. In simulation studies, this weighting estimator showed promising
finite sample performance, reducing bias and mean squared error as compared
with simpler methods.

The proposed IPWM estimator is an extension of the inverse probability
weighting estimator recently proposed by Gravel and Platt (2018) (GP) which
only addresses the misclassification in the outcome. IPWM and GP are
(mathematically) equivalent when the exposure is (assumed to be) measured
without misclassification error.

Like the Gravel and Platt approach, IPWM relies on estimates of sensitivity
and specificity or positive and negative predictive values for the
misclassified variables. In this paper, we used an internal approach where a
portion of subjects would receive error-free (‘gold standard’) measurements on
either or both the outcome and exposure.

We anticipate that in some settings the likelihood may not be fully
identifiable from the data at hand. In these settings, it may be possible to
incorporate external rather than internal information on the misclassification
rates, possibly through a Bayesian approach using prior assumptions about
misclassification probabilities. When validation data is external, however, it
may be necessary to assume misclassification to be independent of covariates
$L$, because external studies seldom consider the same covariates as the main
study (Lyles et al., 2011). External validation approaches also require the
assumption that the misclassification parameters targeted in the validation
sample are transportable to the main study.

An important advantage of the IPWM approach is that the subjects with
validation data need not form a completely random subset. The proposed method
was developed under the assumption that validation data allocation occurs in
an ‘ignorable’ fashion (Rubin, 1976). In practice, it may be that the
researchers have limited control over the validation data allocation
mechanism. For instance, it is conceivable that individuals with specific
indications (e.g., with a realisation of $L$, $B$ or $Z$) are practically
ineligible to be assigned a double measurement of the exposure ($A$ and $B$)
and outcome ($Y$ and $Z$). Further, the estimator also allows for validation
subjects to receive either the double exposure or double outcome measurement.
We simulated data such that subjects have validation data on both the exposure
and outcome variables or on neither. Although this may greatly simplify
analysis and enhance efficiency, in practice it is not necessary to assume
that this condition holds. An interesting scenario is where subjects have
validation data on at most one variable, i.e., on the exposure variable or the
outcome variable but not both. In this case, valid estimation would require
additional modelling assumptions; for example, the error-free outcome variable
cannot then be regressed on the error-free exposure variable.

To accommodate settings where validation data allocation is not completely at
random, we deviated from the semi-parametric bootstrap procedure for variance
estimation proposed by Gravel and Platt. Instead, the non-parametric procedure
we used requires less assumptions regarding the validation subset sampling
procedure. The non-parametric procedure showed good performance in our
simulations.

Whilst we have discussed under what conditions the proposed method
consistently estimates or at least identifies the target quantity, the
assumptions may be untenable in particular settings. Particularly, an
infallible measurement tool for the exposure and outcome that can be performed
on a subset of the data need not always exist. The robustness to deviations of
infallibility is an interesting and important direction for further research.
This is especially relevant where there exists considerable uncertainty about
the tenability of the assumptions that is difficult to incorporate in the
analysis. An obvious and flexible alternative to IPWM is to multiply impute
missing values including absent measurement error-free variables before
implementing IPW (MI+IPW). Although MI+IPW and IPWM may be comparable in terms
of their assumptions, it is yet unclear how they behave under assumption
violations such as misspecification of the outcome model.

An advantageous property of MI+IPW is that it can easily accommodate missing
covariate values. Other alternatives that can accommodate missing covariates
were recently developed by Shu and Yi (2018). Their proposed weighting
estimators simultaneously addresses confounding, misclassification of the
outcome (but not of the exposure) and measurement error on the covariates
under a classical additive measurement error model. The methods can be
implemented using validation data or repeated measurements and use a simple
misclassification model (in which the outcome surrogate is independent of
exposure or covariates given the target outcome) that is suitable for
performing sensitivity analyses.

Another interesting area for further research is where the researchers do have
control over who is referred for further testing by the assumed infallible
measurement tool(s). An obvious choice is to adopt a completely at random
strategy (simple random sampling). However, other referral (sampling)
strategies exist and it is not clear what strategy leads to the most
favourable estimator properties for the given setting.

In summary, we have developed an extension to an existing method, to allow for
valid estimation of a marginal causal OR in the presence of confounding and a
commonly ignored and misunderstood source of bias—joint exposure and outcome
misclassification. The R function mecor::ipwm has been made available to
facilitate implementation (Nab, 2019; Nab et al., 2018).

## Acknowledgements

RHHG was funded by the Netherlands Organization for Scientific Research (NWO-
Vidi project 917.16.430). The views expressed in this article are those of the
authors and not necessarily any funding body.

## References

  * Austin and Stafford (2008) Austin, P. C. and Stafford, J. (2008).  The performance of two data-generation processes for data with specified marginal treatment odds ratios.  Communications in Statistics - Simulation and Computation 37, 1039–1051.
  * Brakenhoff et al. (2018) Brakenhoff, T. B., Mitroiu, M., Keogh, R. H., Moons, K. G., Groenwold, R. H., and van Smeden, M. (2018).  Measurement error is often neglected in medical literature: a systematic review.  Journal of Clinical Epidemiology 98, 89–97.
  * Brenner et al. (1993) Brenner, H., Savitz, D. A., and Gefeller, O. (1993).  The effects of joint misclassification of exposure and disease on epidemiologic measures of association exposure and disease on epidemiologic measures of association.  Journal of Clinical Epidemiology 46, 1195–1202.
  * Brooks et al. (2018) Brooks, D. R., Getz, K. D., Brennan, A. T., Pollack, A. Z., and Fox, M. P. (2018).  The impact of joint misclassification of exposures and outcomes on the results of epidemiologic research.  Current Epidemiology Reports 5, 166–174.
  * Culver et al. (2012) Culver, A. L., Ockene, I. S., Balasubramanian, R., Olendzki, B. C., Sepavich, D. M., Wactawski-Wende, J., Manson, J. E., Qiao, Y., Liu, S., Merriam, P. A., et al. (2012).  Statin use and risk of diabetes mellitus in postmenopausal women in the women’s health initiative.  Archives of internal medicine 172, 144–152.
  * Dawid (1979) Dawid, A. (1979).  Conditional independence in statistical theory.  Journal of the Royal Statistical Society, Series B (Methodological) pages 1–31.
  * Gravel and Platt (2018) Gravel, C. A. and Platt, R. W. (2018).  Weighted estimation for confounded binary outcomes subject to misclassification.  Statistics in Medicine 37, 425–436.
  * Holland (1986) Holland, P. (1986).  Statistics in causal inference.  Journal of the American Statistical Association 81, 945–960.
  * Holland (1988) Holland, P. (1988).  Causal inference, path analysis, and recursive structural equations models.  Sociological Methodology 18, 449–484.
  * Jurek et al. (2008) Jurek, A. M., Greenland, S., and Maldonado, G. (2008).  Brief report: how far from non-differential does exposure or disease misclassification have to be to bias measures of association away from the null?  International Journal of Epidemiology 37, 382–385.
  * Jurek et al. (2006) Jurek, A. M., Maldonado, G., Greenland, S., and Church, T. R. (2006).  Exposure-measurement error is frequently ignored when interpreting epidemiologic study results.  European journal of epidemiology 21, 871–876.
  * Kristensen (1992) Kristensen, P. (1992).  Bias from nondifferential but dependent misclassification of exposure and outcome.  Epidemiology pages 210–215.
  * Leong et al. (2013) Leong, A., Dasgupta, K., Bernatsky, S., Lacaille, D., Avina-Zubieta, A., and Rahme, E. (2013).  Systematic review and meta-analysis of validation studies on a diabetes case definition from health administrative records.  PloS one 8, e75256.
  * Lyles et al. (2011) Lyles, R. H., Tang, L., Superak, H. M., King, C. C., Celentano, D. D., Lo, Y., and Sobel, J. D. (2011).  Validation data-based adjustments for outcome misclassification in logistic regression: an illustration.  Epidemiology 22, 589\.
  * Marcum et al. (2013) Marcum, Z. A., Sevick, M. A., and Handler, S. M. (2013).  Medication nonadherence: a diagnosable and treatable medical condition.  Jama 309, 2105–2106.
  * Nab (2019) Nab, L. (2019).  mecor: Measurement Error Corrections.  R package version 0.1.0. Available from: https://github.com/LindaNab/mecor.git.
  * Nab et al. (2018) Nab, L., Groenwold, R. H., Welsing, P. M., and van Smeden, M. (2018).  Measurement error in continuous endpoints in randomised trials: problems and solutions.  arXiv preprint arXiv:1809.07068 .
  * Neyman et al. (1935) Neyman, J., Iwaszkiewicz, K., and St. Kolodziejczyk (1935).  Statistical problems in agricultural experimentation.  Supplement to the Journal of the Royal Statistical Society 2, 107–180.
  * Ni et al. (2017) Ni, J., Leong, A., Dasgupta, K., and Rahme, E. (2017).  Correcting hazard ratio estimates for outcome misclassification using multiple imputation with internal validation data.  Pharmacoepidemiology and drug safety 26, 925–934.
  * Pearl (2009) Pearl, J. (2009).  Causality: Models, Reasoning and Inference.  Cambridge University Press, New York.
  * R Core Team (2018) R Core Team (2018).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rubin (1974) Rubin, D. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology 66, 688–701.
  * Rubin (1976) Rubin, D. (1976).  Inference and missing data.  Biometrika 63, 581–592.
  * Setoguchi et al. (2008) Setoguchi, S., Schneeweiss, S., MA, M. B., Glynn, R., and Cook, E. (2008).  Evaluating uses of data mining techniques in propensity score estimation: a simulation study.  Pharmacoepidemiology and Drug Safety 17, 546–555.
  * Shu and Yi (2018) Shu, D. and Yi, G. Y. (2018).  Weighted causal inference methods with mismeasured covariates and misclassified outcomes.  Statistics in Medicine .
  * Tang et al. (2013) Tang, L., Lyles, R. H., Ye, Y., Lo, Y., and King, C. C. (2013).  Extended matrix and inverse matrix methods utilizing internal validation data when both disease and exposure status are misclassified.  Epidemiologic Methods 2, 49–66.
  * VanderWeele and Hernán (2012) VanderWeele, T. J. and Hernán, M. A. (2012).  Results on differential and dependent measurement error of the exposure and the outcome using signed directed acyclic graphs.  American journal of epidemiology 175, 1303–1310.
  * Vogel et al. (2005) Vogel, C., Brenner, H., Pfahlberg, A., and Gefeller, O. (2005).  The effects of joint misclassification of exposure and disease on the attributable risk.  Statistics in Medicine 24, 1881–1896.

##  Appendix A I

###### Theorem A.1

For any $a,l$, let

$\displaystyle\varphi(a,l)$
$\displaystyle=\frac{\varphi^{\ast}(a,l)}{\mathbb{E}[\varphi^{\ast}(A,L)|A=a]}~{}~{}\text{and}~{}~{}\varphi^{\ast}(a,l)=\frac{1}{\Pr(A=a|L=l)}.$

If $Y(A)=Y$ (consistency), $(Y(0),Y(1))\mathrel{\perp\\!\\!\\!\perp}A|L=l$
(conditional exchangeability), $\Pr(A=a)>0$ and $\Pr(A=a|L=l)>0$ (positivity)
for all $a$ and every $l$ in the support of $L$, then

$\displaystyle\mathbb{E}[Y(a)]$
$\displaystyle=\mathbb{E}[\varphi(A,L)I(Y=1)|A=a].$

###### Proof A.2

We begin by considering $\mathbb{E}[\varphi^{\ast}(A,L)|A=a]$. By the law of
the unconscious statistician and Bayes’ theorem, we have

$\displaystyle\mathbb{E}[\varphi^{\ast}(A,L)|A=a]$
$\displaystyle=\sum_{l}\frac{\Pr(L=l|A=a)}{\Pr(A=a|L=l)}$
$\displaystyle=\sum_{l}\frac{\Pr(A=a|L=l)\Pr(L=l)}{\Pr(A=a)\Pr(A=a|L=l)}$
$\displaystyle=\frac{1}{\Pr(A=a)}\sum_{l}\Pr(L=l)$
$\displaystyle=\frac{1}{\Pr(A=a)}.$

Hence, for all $a,y$, we have

$\displaystyle\sum_{l}\varphi(a,l)\Pr(Y=y,L=l|A=a)$
$\displaystyle=\sum_{l}\frac{\Pr(Y=y,L=l|A=a)\Pr(A=a)}{\Pr(A=a|L=l)}$
$\displaystyle=\sum_{l}\frac{\Pr(Y=y|A=a,L=l)\Pr(A=a|L=l)\Pr(L=l)}{\Pr(A=a|L=l)}$
$\displaystyle=\sum_{l}\Pr(Y=y|A=a,L=l)\Pr(L=l)$
$\displaystyle=\sum_{l}\Pr(Y(a)=y|A=a,L=l)\Pr(L=l)$ (A.2.1)
$\displaystyle=\sum_{l}\Pr(Y(a)=y|L=l)\Pr(L=l)$ (A.2.2)
$\displaystyle=\Pr(Y(a)=y),$

where (A.2.1) and (A.2.2) hold under consistency and conditional
exchangeability given $L$, respectively. Positivity ensures the weights are
defined/exist. Hence,
$\mathbb{E}[\varphi(A,L)I(Y=1)|A=a]=\sum_{l}\varphi(a,l)\Pr(Y=1,L=l|A=a)=\mathbb{E}[Y(a)]$,
as desired.

###### Corollary A.1

For any $y,a,l$, let

$\displaystyle\varphi(a,l)$
$\displaystyle=\frac{\varphi^{\ast}(a,l)}{\mathbb{E}[\varphi^{\ast}(A,L)|A=a]},~{}~{}\varphi^{\ast}(a,l)=\frac{1}{\Pr(A=a|L=l)},~{}~{}\text{and}$
$\displaystyle\phi(a,l)$
$\displaystyle=\frac{\Pr(Y=1,L=l|A=a)}{\Pr(Z=1,L=l|B=a)}.$

If $Y(A)=Y$, $(Y(0),Y(1))\mathrel{\perp\\!\\!\\!\perp}A|L$ and positivity
holds, then

$\displaystyle\mathbb{E}[Y(a)]$
$\displaystyle=\sum_{l}\varphi(a,l)\Pr(Y=1,L=l|A=a)$
$\displaystyle=\sum_{l}\varphi(a,l)\phi(a,l)\Pr(Z=1,L=l|B=a)$
$\displaystyle=\mathbb{E}[\varphi(B,L)\phi(B,L)Z|B=a].$

##  Appendix B II

###### Theorem B.1

Fix some $s>0$ and let $P^{\ast}=(Ps+1)/(s+2)$ for all $P\in[0,1]$. If
$(P_{0},P_{1})\in(0,1)\times(0,1)$, then

$\displaystyle
1<\frac{P_{1}^{\ast}/(1-P_{1}^{\ast})}{P_{0}^{\ast}/(1-P_{0}^{\ast})}<\frac{P_{1}/(1-P_{1})}{P_{0}/(1-P_{0})}~{}~{}$
$\displaystyle\text{if }P_{1}>P_{0},$ $\displaystyle
1=\frac{P_{1}^{\ast}/(1-P_{1}^{\ast})}{P_{0}^{\ast}/(1-P_{0}^{\ast})}=\frac{P_{1}/(1-P_{1})}{P_{0}/(1-P_{0})}~{}~{}$
$\displaystyle\text{if }P_{1}=P_{0},\text{ and}$ $\displaystyle
1>\frac{P_{1}^{\ast}/(1-P_{1}^{\ast})}{P_{0}^{\ast}/(1-P_{0}^{\ast})}>\frac{P_{1}/(1-P_{1})}{P_{0}/(1-P_{0})}~{}~{}$
$\displaystyle\text{if }P_{1}<P_{0}$

###### Proof B.2

Suppose $(P_{0},P_{1})\in(0,1)\times(0,1)$. If and only if

$\displaystyle\frac{P_{1}^{\ast}/(1-P_{1}^{\ast})}{P_{0}^{\ast}/(1-P_{0}^{\ast})}<\frac{P_{1}/(1-P_{1})}{P_{0}/(1-P_{0})},$
(B.2.3)

then

$\displaystyle\frac{P_{1}s+1}{s+1-P_{1}s}\frac{s+1-P_{0}s}{P_{0}s+1}$
$\displaystyle<\frac{P_{1}}{1-P_{1}}\frac{1-P_{0}}{P_{0}},$
$\displaystyle\frac{P_{1}s+1}{s+1-P_{1}s}\frac{1-P_{1}}{P_{1}}$
$\displaystyle<\frac{P_{0}s+1}{s+1-P_{0}s}\frac{1-P_{0}}{P_{0}}.$

Now, since

$\displaystyle\frac{\partial}{\partial
P}\bigg{\\{}\frac{Ps+1}{s+1-Ps}\frac{1-P}{P}\bigg{\\}}=\frac{(-2P^{2}+2P-1)S-1}{P^{2}(1-(P-1)S)^{2}}<0$

over the interval $(0,1)$ for $P$, it follows that inequality (B.2.3) holds if
$P_{1}>P_{0}$. Also, if $P_{1}>P_{0}$, then, since $\partial/(\partial
P)\\{(Ps+1)/(s+1-Ps)\\}>0$ if $P\in(0,1)$, we have

$1<\frac{P_{1}^{\ast}/(1-P_{1}^{\ast})}{P_{0}^{\ast}/(1-P_{0}^{\ast})}.$

Similar arguments establish the assertion for the case where $P_{1}<P_{0}$. It
is easily verified that if $P_{1}=P_{0}$, then

$\displaystyle\frac{P_{1}^{\ast}/(1-P_{1}^{\ast})}{P_{0}^{\ast}/(1-P_{0}^{\ast})}$
$\displaystyle=\frac{P_{1}s+1}{s+1-P_{1}s}\frac{s+1-P_{0}s}{P_{0}s+1}=1$
$\displaystyle=\frac{P_{1}}{1-P_{1}}\frac{1-P_{0}}{P_{0}}=\frac{P_{1}/(1-P_{1})}{P_{0}/(1-P_{0})},$

as desired.

##  Appendix C III

GP and IPWM were applied to every dataset data in R using the function
mecor::ipwm and the following code:

# GP:
formulasGP <\- list(
Y~Z+B+L1+L2+L3+L4+L5+L6+L7+L8+L9+L10,
B~Z+L1+L2+L3+L4+L5+L6+L7+L8+L9+L10,
Z~L1+L2+L3+L4+L5+L6+L7+L8+L9+L10
)
mecor::ipwm(
formulas=formulasGP, data=data, outcome_true=‘‘Y’’,
outcome_mis=‘‘Z’’, exposure_true=‘‘B’’, exposure_mis=NULL, sp=1e6
)

# IPWM:
formulasIPWM <\- list(
Y~A+Z+B+L1+L2+L3+L4+L5+L6+L7+L8+L9+L10,
A~Z+B+L1+L2+L3+L4+L5+L6+L7+L8+L9+L10,
Z~B+L1+L2+L3+L4+L5+L6+L7+L8+L9+L10,
B~L1+L2+L3+L4+L5+L6+L7+L8+L9+L10
)
mecor::ipwm(
formulas=formulasIPWM, data=data, outcome_true=‘‘Y’’,
outcome_mis=‘‘Z’’, exposure_true=‘‘A’’, exposure_mis=‘‘B’’, sp=1e6
)


# Modeling Institutional Credit Risk with Financial News

Tam Tran-The
MassMutual Data Science
470 Atlantic Ave
Boston, Massachusetts 02210
<EMAIL_ADDRESS>
###### Abstract

Credit risk management, the practice of mitigating losses by understanding the
adequacy of a borrower’s capital and loan loss reserves, has long been
imperative to any financial institution’s long-term sustainability and growth.
MassMutual is no exception. The company is keen on effectively monitoring
downgrade risk, or the risk associated with the event when credit rating of a
company deteriorates. Current work in downgrade risk modeling depends on
multiple variations of quantitative measures provided by third-party rating
agencies and risk management consultancy companies. As these structured
numerical data become increasingly commoditized among institutional investors,
there has been a wide push into using alternative sources of data, such as
financial news, earnings call transcripts, or social media content, to
possibly gain a competitive edge in the industry. The volume of qualitative
information or unstructured text data has exploded in the past decades and is
now available for due diligence to supplement quantitative measures of credit
risk. This paper proposes a predictive downgrade model using solely news data
represented by neural network embeddings. The model standalone achieves an
Area Under the Receiver Operating Characteristic Curve (AUC) of more than 80%.
The output probability from this news model, as an additional feature,
improves the performance of our benchmark model using only quantitative
measures by more than 5% in terms of both AUC and recall rate. A qualitative
evaluation also indicates that news articles related to our predicted
downgrade events are specially relevant and high-quality in our business
context.

##  1 Introduction

Credit risk refers to the possibility of loss resulting from a borrower’s or a
bond issuer’s failure to repay a loan or meet contractual obligations. One
subcategory of credit risk is downgrade risk, which occurs when third-party
rating agencies, such as Moody’s and S&P, lower their ratings on a bond or a
company. For example, a change by S&P from a B to a CCC rating is considered a
downgrade event. Such rating information is used extensively by regulatory
organization, specifically the National Association of Insurance Commissioners
(NAIC), to ensure financial solvency of insurance companies. More precisely,
NAIC requires a company with a higher exposure to risk to hold a higher amount
of capital in reserve. Accurately monitoring rating classes and potential
deterioration events is therefore critical to MassMutual. Receiving enough
notice before any possible impending downgrade event would help the company to
manage its investment portfolio more efficiently by preparing necessary
capital and/or switching holding positions.

In this paper, we design 3 natural language processing (NLP) frameworks to
recognize credit relevant patterns (i.e. future 1-year downgrade events) in
news articles about more than 2.2K companies of interest. The study
demonstrates that with an appropriate methodology, we can achieve an AUC of
more than 80% using solely news data to predict downgrade events and yield a
more than 5% performance gain when this adverse credit signal in media
coverage is added to the quantitative downgrade risk model.

The structure of this paper is as follows: Section 2 provides background on
credit risk modeling and the NLP techniques used to convert text into
meaningful numerical data; Section 3 describes the dataset, benchmark model
and main model development; Section 4 presents performance results and
stability test details of the final model; and Section 5 discusses limitations
of this work and future directions.

##  2 Background

###  2.1 Credit Risk Modeling

The motivation to develop credit risk models stems from the need to construct
quantitative estimates of the amount of economic capital needed to support a
financial institution’s risk-taking activities. Specifically, minimum capital
in reserve is often set in proportion to the risk exposure of a company’s
portfolio. Typical credit risk models take as input the conditions of the
general economy and those of the specific firms in question, and generate as
output a credit quality measure.

###  2.2 Natural Language Processing (NLP) Techniques: Topic Modeling,
Sentiment Lexicons and Neural Net Embeddings

When it comes to using text in a machine learning model, one of the main
challenges is how to represent texts as numerical inputs so that we can feed
them into the model. This project aims to experiment with different
methodologies to represent news articles and evaluate the algorithms’
performances based on the downgrade prediction task. Following are some NLP
techniques we use in the study:

Latent Dirichlet Allocation (LDA). As it is believed that every article we
have is composed of major themes, we use topic modeling to extract the hidden
thematic structure in the text. Topic modeling is a type of statistical model
used to discover the abstract topics that occur in a collection of documents.
Latent Dirichlet allocation (LDA), a generative probabilistic model, is a
common algorithm of topic model. The basic idea is documents are represented
as random mixtures over latent topics, where each topic is characterized by a
distribution over words. Mathematically, this can be formulated as:

$P(\theta,\textbf{z},\textbf{w}|\alpha,\beta)=P(\theta|\alpha)\prod_{n=1}^{N}P(z_{n}|\theta)P(w_{n}|z_{n},\beta)$

Given $\alpha$, a parameter vector on the per-document topic distributions and
$\beta$, a parameter vector on the per-topic word distributions, we are
finding the joint distribution of a topic mixture $\theta$, a set of N topics
z, and a set of N words w (?). In the context of our study, each news article
is represented as a set of topic probabilities and the number of topics is a
tunable parameter. To evaluate an LDA model, we use coherence measure, which
gives an estimate of how well each topic can be represented as a composition
of parts that can be combined (?).

Sentiment Lexicons. There is a growing body of sentiment lexicons, or
affective word lists, to examine the tone and sentiment of textual data. Some
lexicons we experiment with are:

  * •

Loughran and McDonald, designed to particularly reflect tone in financial
texts (?)

  * •

VADER, specifically attuned to sentiment in microblog-like contexts (?)

  * •

AFINN, also constructed based on micro-blog posts including Internet slang
words (?)

  * •

SentiWordNet and OPINION, for general sentiment analysis

Neural Net Embeddings. Based on the underlying idea that “a word is
characterized by the company it keeps,” each word is represented by an
embedding, or a vector of continuous numbers so that words of similar semantic
(relative to the task) are closer to one another in the vector space. By using
neural network on a supervised/unsupervised task, the resulting
weights/parameters that have been adjusted to minimize loss on the task are
the embeddings.

  * •

Doc2Vec Document Embeddings. Doc2Vec is an unsupervised framework that learns
fixed-length feature representations from variable-length pieces of texts. In
the framework, a vector representation is trained by stochastic gradient
descent and backpropagation to be useful for predicting the next word in the
sentence. The Doc2Vec model we use in this study is a distributed memory one,
where each paragraph vector acts as a memory that remembers what is missing
from the current context (?). After being trained, the paragraph vectors can
be used as features for the paragraph and be fed directly into conventional
machine learning algorithms such as logistic regression.

  * •

fastText Word Embeddings. Inspired by the same hypothesis as Doc2Vec, fastText
word representations are trained to predict words that appear in their
contexts. What distinguishes fastText model is it takes into account
morphology. More precisely, the model represents each word by a sum of its
character n-grams, which allows us to compute representations for words that
did not appear in the training data and proves to be helpful while working
with morphologically rich language (?).

##  3 Modeling Downgrade Risk

Figure 1:  Overview of modeling pipeline. For details of benchmark model, see
Section 3.2. For details of news model, see Section 3.3. Note that features
for the final model include the outputted probability from news model and all
quantitative features from benchmark model.

###  3.1 Data

News. For each company at a point in time (daily), we have all news articles
corresponding to the record and a label indicating whether the company is
going to downgrade within future 1-year period.

  * •

Data source. We have access to news data source from Thomson Reuters via
CreditEdge’s API provided by Moody’s. This source brings us daily article
pieces at company level. Reuters articles in the dataset could take various
forms: news stories about company’s activities (e.g. merger & acquisition,
board assignment, etc.), market snapshots (e.g. individual stock movements,
etc.), or quarterly earnings summaries. There are approximately 21K
observations identified at company and date level, covering more than 2.2K
Moody’s unique permanent company identifiers (PID) and spanning from mid
November 2017 to early October 2019 (although 99% of these articles are in
2019).

  * •

Data preprocessing. Since we are working in a very narrow domain (financial
news of a company of interest) and the data is particularly noisy, we take
multiple cleaning and preprocessing steps. First, we remove small machine-
generated articles, such as NYSE/NASDAQ/AMEX order imbalance, articles that
only include video links or articles whose accuracy hasn’t been verified by
Reuters. Second, we remove unnecessary headers, footers, HTML metadata, and
lowercase all texts. Third, we manually read through a handful of articles and
flag which format tends to mention multiple companies in the text and which
tends to talk about a single company. Last but not least, for pieces that
cover information of multiple companies, we use a combination of named entity
recognition and fuzzy matching techniques, as well as rule-based exact
matching to extract only sentences about the company of relevance. A few more
preprocessing layers, such as tokenization, stemming and stop words removal,
are added depending on which modeling approach we take. By performing a
thorough text preprocessing, we hope to increase the signal to noise ratio
which suggests news is more likely to have a material impact on our downgrade
prediction. To ensure there is only one row of data per company per day, for
companies that have multiple articles on a single day, we concatenate these
articles together.

Ratings. Credit ratings for each company are available from both Moody’s and
S&P. We take the worst combination between the two rating systems for a
company at a point in time (daily) to derive current rating and the worst
rating within the next future 1 year, from which we determine whether that
company has an impending downgrade event or not.

Table 1: Downgrade events by year downgrade | year | count
---|---|---
0 | 2017 | 42
2018 | 152
2019 | 20,519
1 | 2019 | 324

For our benchmark model, all features are variations of quantitative credit
measures provided by third-party rating agencies and risk management
consultancy companies - this is further discussed in Section 3.2.

###  3.2 Benchmark Model

The benchmark downgrade model of this study is a logistic regression. For each
training data point, we have a vector of features $X$ and an observed class
$Y$, where $X$ provides quantitative metrics about a company on a daily basis
and $Y=1$ indicates that the company is downgrading within 1 year from the
date of the data observation. This setup ensures that both our benchmark and
news models have the same level of data frequency and can be compared and
combined later. Assuming that $P(Y=1|X=x)=p(x;\theta)$ for some function $p$
parameterized by $\theta$, what we want to model is:

$P(Y=1|X=x,\theta)=\frac{1}{1+e^{-\theta x}}$

The model is trained on 9 features that are variations in terms of either term
structures (i.e. 1-year, 5-year) or transformation (i.e. lag, diff) of
following quantitative measures purchased from risk management providers:

  * •

The probability that a firm will default over 1 year based on company-specific
attributes, industry-related measures and relevant macro-economic factors

  * •

The probability that a firm will downgrade within 1 year based on market
implied ratings and rating outlooks

  * •

Historical credit rating for any firm at a point in time

Since our data is highly imbalanced, where the occurrence of downgrade events
accounts for 1.5% of the dataset, we employ SMOTE algorithm to replicate
observations from the minority class. Overall, this benchmark model achieves
an AUC of 82.7% and a recall rate of 69.5%.

###  3.3 Model Development

Figure 1 is an overview of our modeling pipeline. Throughout the training
process, we deliberately employ logistic regression due to the algorithm’s
great interpretability which is a focus of our business users who work in a
highly regulated industry. This section details 3 different approaches we take
to translate unstructured texts into meaningful numerical data that can be fed
into a logistic regression model. We evaluate these NLP methodologies on our
downstream task at hand, which is downgrade prediction, to select the best
news model that can be incorporated into and improve the existing benchmark
model.

Approach 1: Lexicon-based Sentiment and Topic Scores. We first experiment with
one of the simplest sentiment analysis approaches, which is to compare words
of an article against a labeled word list, where each word has been scored for
valence. Based on 5 different lexicons (Loughran & McDonald, VADER, AFINN,
SentiWordNet, and OPINION), we compute 5 sentiment scores for each news
article. Specifically, we check whether a word in the text exists in each
positive or negative word list, then count the frequency of that word. For
positive words, we do not count those that have a negation in one of three
places preceding it. The sentiment score is calculated as
$\frac{\\#positives-\\#negatives}{\\#positives+\\#negatives}$. Additionally,
we train an LDA model with 10 different topics where each topic is a
combination of keywords and each keyword contributes a certain weight to the
topic. We decide on 10 as the optimal number of topics by running multiple LDA
models with number of topics ranging from 5 to 40 and picking the one that
gives the highest coherence value, as shown in Figure 2. Two features that
indicate the number of articles for each company at a point in time and
whether an article contains any variation of the word “downgrade” are also
created. In the end, our news model is trained on 17 features, including 5
sentiment scores, 10 topic probabilities, and 2 count indicators. This news
model standalone has an unimpressive AUC of 59.8%. Although sentiment and
topic scores are helpful in giving us a basic analysis of our news data at
hand, these features are far from sufficient for accurate downgrade
prediction.

Figure 2:  Coherence values of LDA model with different topic numbers. The
optimal number is 10 topics with coherence value of 0.704.

Approach 2: Doc2Vec Embeddings. In this approach, we train a distributed-
memory Doc2Vec model based on the vocabulary of our training dataset over 20
epochs and generate a document embedding for each news article. We choose to
use Doc2Vec paragraph vectors since they are learned from unlabeled data and
theoretically can work well for tasks that have small amount of labeled data
(?), which is our case in this study. These embeddings are taken as input into
a logistic regression model, which results in an AUC of 71.6%.

Approach 3: fastText Embeddings. Taking advantage of transfer learning, we use
1M word vectors pre-trained by a team of Facebook researchers on Wikipedia
2017, UMBC webbase corpus and statmt.org news dataset (?). Since the fastText
model was trained on massive data sources, their representations perform very
well at transferring to other NLP problems and improve the generalization of
models learned on limited amount of data. Unlike Doc2Vec which produces
representations at document level, fastText generates a 300-dimensional vector
for each word. To create an embedding for the whole document, we take the
average of embeddings for words contained in it. With these final document
embeddings fed into a logistic regression, our news model achieves an AUC of
80.9%, which is the best performance for this standalone model so far.

Table 2: Performance of standalone news model approach | AUC (%)
---|---
sentiment scores and LDA topic scores | 59.8
Doc2Vec embeddings | 71.6
fastText embeddings (averaging out) | 80.9

Based on the performance results, we decide to employ fastText embeddings as
the main feature of our news model. The downgrade probability outputted from
this is then combined with other quantitative measures in the benchmark,
coupled with SMOTE oversampling technique, to build a final augmented model.
Results of this final model are further explored in Section 4.

##  4 Final Model Results and Validation

###  4.1 Quantitative Evaluation

The following results are based on the evaluation on a 20% hold-out test set
that are untouched, unseen throughout the training process. Both the benchmark
and final models are tested on this same set to ensure fair comparison.

AUC. One of the main metrics we use to compare performances of different
models is AUC, which tells how good a model is at distinguishing classes (the
higher the number, the better). Figure 3 indicates that a downgrade model
using news alone can achieve a test AUC of 80.9% and there is a 5% gain when
adding this news probability on top of the existing model (from 82.7% for
benchmark to 87.7% for final model).

Figure 3:  AUC comparison between benchmark, standalone news, and final
models. Final model achieves the highest value, 87.7%, resulting in a 5% gain
compared with benchmark model.

Recall. Since the cost of a false negative in our business context (an
impending downgrade event that the model is not able to capture) is very
steep, we put an emphasis on optimizing the model’s recall rate. Although a
standalone news model produces a humble recall rate of 54.2%, incorporating
news into benchmark model improves recall rate from 69.5% to 74.6% on the test
set, resulting in a 5.1% gain in recall.

Cumulative Gains. Given limited human and time resources that can be dedicated
to the task of analyzing companies in the list of anticipated downgrades, we
would like to be able to capture as many downgrades as possible using as few
test cases as possible. Figure 4 indicates that with news information added,
analyzing the top 10% companies with the highest predicted downgrade
probabilities can achieve a recall rate of 74.6%, which is a 27.1% improvement
compared with the benchmark model. Similarly, analyzing the top 20% companies
using the final model results in a recall rate of 83.1%, an increase of 13.6%
compared with benchmark.

Figure 4:  Cumulative gains chart. Given limited resources, the final model is
able to provide a more optimal list of companies to analyze.

Overall, our model with the additional downgrade probability of news not only
has higher accuracy rate (in terms of AUC and recall) but is also more
efficient by providing a more optimal list of companies to analyze given
restrained resources.

###  4.2 Qualitative Evaluation

To enhance the pragmatic sense of our model, we carefully examine all news
articles associated with the true positives (actual downgrade events that the
model is able to capture) and false negatives (actual downgrade events that
the model is unable to capture) in our test set. Investigation shows that the
model can pick up high-quality and very relevant article pieces to infer a
downgrade event of a company.

True Positive Samples. There are 44 true positives corresponding with 16
unique companies in our test data. Following is a representative list of
articles appearing in our true positive cases. Each example is shown at
company level, contains only keywords due to space limit, and gives a good
summary of the entire article piece.

  * •

Company A: close factories; cut down 12,000 jobs; has been among the hardest
to be hit by the trade war so far

  * •

Company B: lose clients and get sued following misconduct revelations

  * •

Company C: challenged by [an investor] for stretching itself financially to
buy rival oil driller

  * •

Company D: struggle with a host of issues; cut its dividend and report a
wider-than-expected loss in its main engineering and construction unit

  * •

Company E: hire restructuring firms and may choose to seek bankruptcy
protection

  * •

Company F: plan to wind down its dress-barn retail operations, resulting in
the closure of about 650 stores

  * •

Company G: anticipate having discussions on constructive basis relating to its
underperformance

  * •

Company H: make sophisticated missiles that use rare earth metals in their
guidance systems, and sensors

  * •

Company I: face lawsuit over art fraud; had been the willing auction house
that knowingly and intentionally made the fraud possible

Although none of the articles explicitly mention that a company is being
downgraded, the majority of them provide valuable insights into a company’s
financial health or how the company is perceived in the market. These signals
could potentially be considered as leading indicators toward a downgrade. The
model seems to be able to pick up subtle and indirect signs of an impending
downgrade event as well. For example, in the case of Company H, the article
refers to the company’s heavy dependence on rare earth metals. This does not
appear to relate to downgrade at first glance; however, more examination
suggests that rare earth metal has become a political spotlight recently due
to the trade tension between the U.S, who considers this mineral critical to
the country’s economic and national security, and China, who is the largest
producer and manufacturer of this element in the world. Company H indeed was
downgraded on 07-31-2019.

False Negative Samples. The test set includes 15 false negatives associated
with 6 distinct companies. After review, we see these cases fall under either
of two following categories:

  * •

The news articles follow unpopular format structures compared with others in
the dataset, which leads to their noisy quality even after preprocessing.

  * •

The news articles appear in this false negative list are about companies that
also exist in the true positive list. However, these pieces have a
positive/neutral tone, were published shortly earlier, and might revolve
around a different activity mentioned in the articles of the true positives.
For example:

    * –

Company I’s shareholders approve proposed acquisition; clients are pleased
with the company’s help to file their first lawsuit against the government
(This piece was published on 06-21-2019. The one in true positive list was
published on 06-25-2019. Company I’s downgrade event happened on 09-17-2019).

    * –

Company C could divest most or all of [an investment vehicle] after buying out
[another company]. (This piece was published on 06-14-2019. The one in true
positive list was published on 06-27-2019. Company C’s downgrade event
happened on 08-01-2019).

Our takeaways here are that bad news can “outweigh” good/neutral news when it
comes to contributing to the prediction of downgrade and companies’ ratings
could take downturn very shortly after emergence of some negative sentiments
about them in the market.

###  4.3 Model Robustness

To evaluate the robustness of our final model’s performance gain, we run 100
experiments using 100 different random seeds that are involved in the
training/splitting, cross-validation and SMOTE oversampling processes. The
performance gain in AUC from adding news to the benchmark model is positive
100 times out of 100 experiments. The mean value of these gains is 6.3% with a
standard error of 1.6%, as shown in Figure 5.

Figure 5:  Histogram of the performance gain in 100 experiment runs. The gain
is positive 100 times, with mean value 6.3% and standard error 1.6%.

This suggests that our positive performance gain is reliable and the 5%
increase in AUC mentioned in Section 4.1 is still on the lower side. As we
have access to a larger set of data, the standard error will theoretically be
pushed closer to 0.

##  5 Conclusions and Future Work

In the hope of obtaining unique insights into companies’ financial health that
are not available in traditional quantitative credit measures, we train a
downgrade risk model on solely news information. We demonstrate that news
coverage, if represented appropriately in the data, can help detect adverse
credit signal and considerably improve the performance of existing model that
is trained on conventional credit measures.

There are multiple avenues for research opened from this work. Throughout our
modeling development, we notice how we preprocess the text data could
introduce a material change in our final model performance. Extracting only
relevant information about the company of interest plays a key role in
increasing the gain. Thus, we would like to further explore more advanced
techniques regarding information extraction, such as coreference resolution to
find sentences that do not exactly include a company’s name but still refer to
the same entity, or sentence segmentation to pick out sentences when the
boundary is ambiguous. In addition, because the news model is built solely
based on online articles, it is inherently subject to media bias. Fact-
checking news or tackling different forms of bias within mass media is outside
the scope of this study but offers an exciting opportunity for future
research.

##  6 Acknowledgement

The author is grateful for helpful discussions about model development made by
Zizhen Wu and Jasmine Geng, for knowledge sharing about the raw datasets made
by Yi Li and Yi Wang, and for paper review by Nailong Zhang, Jasmine Geng,
Adam Fox and Sears Merritt.

## References

  * [Blei et al. 2003] Blei, D. M.; Ng, A. Y.; Jordan, M. I.; and Lafferty, J.  2003\.  Latent dirichlet allocation.  Journal of Machine Learning Research 3:2003.
  * [Bojanowski et al. 2016] Bojanowski, P.; Grave, E.; Joulin, A.; and Mikolov, T.  2016\.  Enriching word vectors with subword information.  CoRR abs/1607.04606.
  * [Hutto and Gilbert 2014] Hutto, C. J., and Gilbert, E.  2014\.  Vader: A parsimonious rule-based model for sentiment analysis of social media text.  In Adar, E.; Resnick, P.; Choudhury, M. D.; Hogan, B.; and Oh, A. H., eds., ICWSM.  The AAAI Press.
  * [Le and Mikolov 2014] Le, Q., and Mikolov, T.  2014\.  Distributed representations of sentences and documents.  In Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32, ICML’14, II–1188–II–1196.  JMLR.org.
  * [Loughran and McDonald 2011] Loughran, T., and McDonald, B.  2011\.  When is a liability not a liability? textual analysis, dictionaries, and 10‐ks.  Journal of Finance 66(1):35–65.
  * [Mikolov et al. 2018] Mikolov, T.; Grave, E.; Bojanowski, P.; Puhrsch, C.; and Joulin, A.  2018\.  Advances in pre-training distributed word representations.  In Proceedings of the International Conference on Language Resources and Evaluation (LREC 2018).
  * [Nielsen 2011] Nielsen, F. Å.  2011\.  A new ANEW: evaluation of a word list for sentiment analysis in microblogs.  CoRR abs/1103.2903.
  * [Röder, Both, and Hinneburg 2015] Röder, M.; Both, A.; and Hinneburg, A.  2015\.  Exploring the space of topic coherence measures.  In Proceedings of the Eighth ACM International Conference on Web Search and Data Mining, WSDM ’15, 399–408.  New York, NY, USA: ACM.