# A joint Bayesian framework for missing data and measurement error using
integrated nested Laplace approximations

Emma Skarstein  Department of Mathematical Sciences, Norwegian University of
Science and Technology, Norway  Sara Martino  Department of Mathematical
Sciences, Norwegian University of Science and Technology, Norway  Stefanie
Muff  Department of Mathematical Sciences, Norwegian University of Science and
Technology, Norway  Centre for Biodiversity Dynamics, Norwegian University of
Science and Technology, Norway

###### Abstract

Measurement error (ME) and missing values in covariates are often unavoidable
in disciplines that deal with data, and both problems have separately received
considerable attention during the past decades. However, while most
researchers are familiar with methods for treating missing data, accounting
for ME in covariates of regression models is less common. In addition, ME and
missing data are typically treated as two separate problems, despite practical
and theoretical similarities. Here, we exploit the fact that missing data in a
continuous covariate is an extreme case of classical ME, allowing us to use
existing methodology that accounts for ME via a Bayesian framework that
employs integrated nested Laplace approximations (INLA), and thus to
simultaneously account for both ME and missing data in the same covariate. As
a useful by-product, we present an approach to handle missing data in INLA,
since this corresponds to the special case when no ME is present. In addition,
we show how to account for Berkson ME in the same framework. In its broadest
generality, the proposed joint Bayesian framework can thus account for Berkson
ME, classical ME, and missing data, or for any combination of these in the
same or different continuous covariates of the family of regression models
that are feasible with INLA. The approach is exemplified using both simulated
and real data. We provide extensive and fully reproducible Supplementary
Material with thoroughly documented examples using R-INLA and inlabru.

Keywords: Bayesian joint model, Berkson measurement error, classical
measurement error, integrated nested Laplace approximation (INLA), missing
data

##  1 Introduction

Missing data and measurement error (ME) are both cases of information loss.
This loss means that at the very least, the level of uncertainty in our
inferences increases, and, depending on the underlying mechanism, the error
and missing data might also be biasing our conclusions in ways that are not
immediately obvious (Fuller,, 1987; Little and Rubin,, 1987; Carroll et al.,,
2006). Both problems are unavoidable in any field that deals with
observational data. A common example is missing responses in surveys, where ME
may also be present due to imprecise or untrue answers. Other examples can be
found in exposure assessment in epidemiology (Heid et al.,, 2004; Goldman et
al.,, 2011), phenotypic ME in quantitative genetics (Ponzi et al.,, 2018), or
errors in expert coded data in political science (Marquardt,, 2020). It is
therefore important to understand how to evaluate the impact of ME and missing
data on statistical inference, as well as how to adjust for both to ensure
valid inference.

Missing data and its implications for statistical analyses have been studied
extensively (_e.g._ , Little and Rubin,, 1987; van Buuren,, 2018; Carpenter
and Smuk,, 2021). The missingness of values in a covariate can be classified
based on three main mechanisms that generate the missingness. The first type
is referred to as _missing completely at random_ (MCAR), arising when the
missingness occurs entirely at random, not depending on any variables,
observed or unobserved. The second type, referred to as _missing at random_
(MAR), occurs when the missingness depends on some other variable(s) which are
themselves entirely observed. The last type is referred to as _missing not at
random_ (MNAR), and occurs when the missingness depends on the values of the
variables that are themselves unobserved (Little and Rubin,, 1987). In a
separate branch of research, ME has also been studied extensively (_e.g._ ,
Fuller,, 1987; Gustafson,, 2003; Carroll et al.,, 2006; Buonaccorsi,, 2010;
Yi,, 2017). Apart from increasing uncertainty, ME typically introduces bias
into parameter estimates, whereas the type and direction of bias heavily
depends on the actual error mechanism. As a major distinction, we can classify
ME into classical and Berkson error, which influence parameter estimates in
fundamentally different ways and must therefore be treated specifically (see,
_e.g._ , Berkson,, 1950; Carroll et al.,, 2006). There may be cases where the
ME does not influence the conclusions, but this cannot be known without
inspection. It is therefore recommended to model the ME or carry out a
simulation study when one suspects ME in the data in order to get more insight
into the effect on estimates and their uncertainty (van Smeden et al.,, 2019;
Keogh et al.,, 2020; Innes et al.,, 2021).

When dealing with ME or missing data, Bayesian methods enable us to specifiy
models that mirror the data generating process, and we can easily include
models that explain these errors in the data jointly along with the actual
models we are interested in. An imputation model, for example, can be fit as
part of a hierarchical Bayesian model (Erler et al.,, 2016; Ma and Chen,,
2018), and thus propagates the added uncertainty appropriately through the
levels. Similarly, if there is ME, we may specify an error model describing
the relation between the correct variable and the observed one, and
incorporate any information we have about the ME through priors (Richardson
and Gilks,, 1993; Dellaportas and Stephens,, 1995; Muff et al.,, 2015).

While most of the literature on missing data and ME does not overlap, the
similarities between the two problems have recently been pointed out. As an
example, Cole et al., (2006) present multiple imputation for ME correction. In
another step forward, Blackwell et al., (2017) treat ME and missing data in
the same framework by using the observation that missing data is a limiting
case of ME, corresponding to the case when the error variance tends to
infinity, that is, the complete absence of information (_i.e._ , missingness).
Soon after, Goldstein et al., (2018) described a hierarchical Bayesian model
with a classical ME and imputation model that allows the covariate with
missing values to depend on other variables, and Keogh and Bartlett, (2021)
considered ME as a missing data problem by demonstrating multiple ways to
utilize this link. Noghrehchi et al., (2020) combined multiple imputation with
functional methods for error correction in a joint framework, and most
recently van Smeden et al., (2021) compared four different methods that
address both ME and missing data, among them a Bayesian model. The integrated
nested Laplace approximation framework enables Bayesian inference for complex
hierarchical models and has over the past decade become a popular tool for a
wide variety of statistical models (Rue et al.,, 2009; Martino and Riebler,,
2020; Gómez-Rubio,, 2020). In the context of INLA, both missing data and ME in
covariates are a challenge, as the covariates are part of the latent Gaussian
Markov random field (GMRF), and INLA requires the values in this latent field
to be fully observed (Gómez-Rubio et al.,, 2022). No joint treatment of
missing data and ME in INLA has previously been proposed, although methodology
for missing data and ME has been suggested separately. In Berild et al.,
(2022), importance sampling is combined with INLA, which then allows for
missing data imputation. However, in scenarios where many observations are
missing this becomes unfeasible. On the other hand, Gómez-Rubio et al., (2022)
present a general approach by defining a covariate imputation model as a
latent effect with a GMRF structure. In the case of covariate ME, (Muff et
al.,, 2015) show how classical and Berkson ME in continuous covariates can be
incorporated into the INLA framework.

Here, we propose to take advantage of the ME and missing data connection to
provide joint ME and missing data modeling using INLA. To this end, we use the
interpretation of missingness as an extreme case of classical ME in order to
employ existing ME models in INLA for missing data imputation. In terms of
missing data, our approach is less general than the one by Gómez-Rubio et al.,
(2022), since their approach also allows for modelling of the missingness
mechanism, which is necessary when the missingness is MNAR. However, in many
practical applications it is reasonable to make the assumption that the
missingness is MAR or even MCAR, in which case the missingness mechanism is
ignorable. The model presented in this paper is therefore well suited for
those cases, and has the advantage of being able to account for ME
simultaneously. The proposed Bayesian model can be fit in any framework for
Bayesian modeling, but we especially highlight the necessary adjustments that
must be made in INLA. This paper is organized as follows. In Section 2, we
describe the Bayesian model that will be used to model ME and missing data in
covariates. In Section 3 we explain how the model is implemented in INLA. In
Section 4, we illustrate the model and its implementation with three examples.
Lastly, Section 5 provides a discussion of the method, results, limitations
and extensions.

##  2 Model specification

###  2.1 Classical measurement error

In the case of classical ME, we are interested in $n$ observations of a
variable $\bm{x}=(x_{1},x_{2},\dots,x_{n})^{T}$, but we actually observe (in
vector notation) $\bm{w}=\bm{x}+\bm{u}_{c}$, with an additive noise term
$\bm{u}_{c}\sim\mathcal{N}(\bm{0},\tau_{u_{c}}\bm{D}_{u_{c}})$ that is
independent of $\bm{x}$. The subscript $c$ is used here to indicate that it is
the error term in a _classical_ error model,
$\tau_{u_{c}}=1/\sigma_{u_{c}}^{2}$ denotes the precision of the error term,
and $\bm{D}_{u_{c}}$ may capture heteroscedasticity or dependence structures,
but in the simplest case $\bm{D}_{u_{c}}$ is the identity matrix, implying
constant ME precision across all observations. A characteristic of classical
ME is that the variance of the observed variable $\bm{w}$ can be additively
split into $\sigma^{2}_{w}=\sigma_{x}^{2}+\sigma_{u_{c}}^{2}$, thus the
variance of $\bm{x}$ will be smaller than the variance of the misobserved
version $\bm{w}$.

When $\bm{x}$ is a covariate in a regression model, a Bayesian hierarchical
model to account for classical ME in the variable $\bm{x}$ can be formulated
as (see _e.g._ , Richardson and Gilks,, 1993; Dellaportas and Stephens,, 1995;
Muff et al.,, 2015)

$\displaystyle\bm{\eta}$
$\displaystyle=\beta_{0}\bm{1}+\beta_{x}\bm{x}+\bm{Z}\bm{\beta}_{z}\ ,$ (1)
$\displaystyle\bm{w}$ $\displaystyle=\bm{x}+\bm{u}_{c}\ ,$
$\displaystyle\bm{u}_{c}\sim\mathcal{N}(\bm{0},\tau_{u_{c}}\bm{D}_{u_{c}})\ ,$
(2) $\displaystyle\bm{x}$
$\displaystyle=\alpha_{0}\bm{1}+\widetilde{\bm{Z}}\bm{\alpha}_{z}+\bm{\varepsilon}_{x}\
,$
$\displaystyle\bm{\varepsilon}_{x}\sim\mathcal{N}(\bm{0},\tau_{\varepsilon_{x}}\bm{D}_{\varepsilon_{x}})\
,$ (3)

where the components are:

Regression model (1) :

$g(\bm{\mu})=\bm{\eta}$ is the linear predictor in a generalized linear model
(GLM), given the correct covariate values for $\bm{x}$, as well as other
covariates observed without error, stored in the matrix $\bm{Z}$. Here,
$g(\cdot)$ is the link function and $\bm{\mu}$ is the mean vector. The
coefficients $\beta_{0}$, $\beta_{x}$ and $\bm{\beta}_{z}$ are the intercept
and the slopes for $\bm{x}$ and $\bm{Z}$, respectively, and these may be given
Gaussian priors centered at $0$ with small precisions.

Classical error model (2) :

The error model describes how the covariate with ME is related to the correct
variable, where $\bm{u}_{c}$ is the error in the observed variable $\bm{w}$,
and $\bm{x}$ and $\bm{u}_{c}$ are assumed to be independent. Importantly, the
prior of the precision $\tau_{u_{c}}$, typically a gamma prior, must be chosen
based on expert knowledge, validation data or repeated measurements to ensure
identifiability (Gustafson,, 2005, 2010).

Imputation model (3) :

Describes how the covariate with error might depend on other fully observed
covariates, stored in the matrix $\widetilde{\bm{Z}}$ (which may or may not be
equal to $\bm{Z}$). Here, $\bm{\varepsilon}_{x}$ is the residual term in the
model for the covariate $\bm{x}$, and its precision $\tau_{\varepsilon_{x}}$
should be given an carefully chosen prior. The intercept $\alpha_{0}$ and
slopes $\bm{\alpha}_{z}$ for the covariates in $\widetilde{\bm{Z}}$ may be
given Gaussian priors centered at $0$ with small precisions. In the ME
literature, Model (3) is generally referred to as the _exposure_ model.

###  2.2 Missing data in relation to classical measurement error

A convenient feature of the ME model in Section 2.1 is that it can at the same
time be used to account for missing values. Imagine that the classical ME
variance $\sigma_{u_{c}}^{2}$ exists on a continuum where we have an
observation with no error when the variance is zero ($w_{i}=x_{i}$), and as
the variance increases our ME gradually increases, to the point where we have
no information about the observation, which is analogous to missing data
(Figure 1). We can thus think of missing data as an extreme case of ME. Using
this, we can reflect the degree of error through the prior of the classical
error variance: either there is no error ($\sigma_{u_{c}}^{2}=0$), some finite
error ($0<\sigma_{u_{c}}^{2}<\infty$), or an observation is completely missing
($\sigma_{u_{c}}^{2}\rightarrow\infty$). We cover the technical aspects of
this in Section 3.2. Note that Figure 1 shows the ME variance
$\sigma_{u_{c}}^{2}$ for illustration purposes, whereas we use the precision
$\tau_{u_{c}}=1/\sigma_{u_{c}}^{2}$ in the rest of the paper.

Figure 1: The scale of classical measurement error variance. The figure is
adapted with permission from Blackwell et al., (2017) (copyright 2015 by the
authors).

###  2.3 Berkson measurement error

The Berkson ME (Berkson,, 1950) is fundamentally different from the classical
ME. While the classical ME can be thought of as random noise, the Berkson ME
arises for example in situations where all members of a population are
assigned the same value for some property, when in reality the value varies
between the members. Berkson error is quite typical in experimental setups,
where a variable is intended to be fixed at a given value, such as a
concentration or a specific time point, when in reality there is some level of
variation due to practical deviations from the target value (see _e.g._ ,
Buonaccorsi and Lin,, 2002). The Berkson error model can be formulated as
$\bm{x}=\bm{w}+\bm{u}_{b}$, where $\bm{x}$ is as before the correct value,
$\bm{w}$ is the observed value and
$\bm{u}_{b}\sim\mathcal{N}(\bm{0},\tau_{u_{b}}\bm{D}_{u_{b}})$, with Berkson
error precision $\tau_{u_{b}}=1/\sigma_{u_{b}}^{2}$ that scales the precision
matrix $\bm{D}_{u_{b}}$. Assuming $\bm{w}$ and $\bm{u}_{b}$ to be independent,
the variance of the correct obervations is additively split like
$\sigma_{x}^{2}=\sigma_{w}^{2}+\sigma_{u_{b}}^{2}$. A property of the Berkson
ME is thus that the variance of the observed values $\bm{w}$ will be smaller
than the variance of the correct value $\bm{x}$, as opposed to the classical
error situation. As with the classical ME, we can set up a Bayesian
hierarchical model that can be used to account for Berkson ME in the variable
$\bm{x}$ (Muff et al.,, 2015),

$\displaystyle\bm{\eta}$
$\displaystyle=\beta_{0}\bm{1}+\beta_{x}\bm{x}+\bm{Z}\bm{\beta}_{z}\ ,$ (4)
$\displaystyle\bm{x}$ $\displaystyle=\bm{w}+\bm{u}_{b}\ ,$
$\displaystyle\bm{u}_{b}\sim\mathcal{N}(\bm{0},\tau_{u_{b}}\bm{D}_{u_{b}})\ ,$
(5)

where Model (4) is the regression model of interest (and corresponds to Model
(1) above), and Model (5) is the Berkson ME model with independent error term
$\bm{u}_{b}$. Note that we do not need an imputation model, as $\bm{x}$ is
defined conditionally on the observations $\bm{w}$ in the error model (5).

###  2.4 A generalized measurement error model

In some cases both error types occur together, for instance in exposure
assessment in epidemiology (see _e.g._ , Reeves et al.,, 1998; Heid et al.,,
2004; Deffner et al.,, 2018). By combining the hierarchical classical error
model given by components (1) – (3) with the Berkson ME model (4) – (5) it is
possible to handle classical and Berkson ME in the same covariate. We again
denote the correct variable as $\bm{x}$ and the observed version by $\bm{w}$,
but we need to introduce an additional latent variable $\bm{r}$ (Muff et al.,,
2017), where the relation between $\bm{r}$ and the correct covariate $\bm{x}$
is the Berkson error model $\bm{x}=\bm{r}+\bm{u}_{b}$. This leads to the
combined model formulated as

$\displaystyle\bm{\eta}$
$\displaystyle=\beta_{0}\bm{1}+\beta_{x}\bm{x}+\bm{Z}\bm{\beta}_{z}\,,$ (6)
$\displaystyle\bm{x}$ $\displaystyle=\bm{r}+\bm{u}_{b}\,,$
$\displaystyle\bm{u}_{b}\sim\mathcal{N}(\bm{0},\tau_{u_{b}}\bm{D}_{u_{b}})\ ,$
(7) $\displaystyle\bm{w}$ $\displaystyle=\bm{r}+\bm{u}_{c}\,,$
$\displaystyle\bm{u}_{c}\sim\mathcal{N}(\bm{0},\tau_{u_{c}}\bm{D}_{u_{c}})\ ,$
(8) $\displaystyle\bm{r}$
$\displaystyle=\alpha_{0}\bm{1}+\widetilde{\bm{Z}}\bm{\alpha}_{z}+\bm{\varepsilon}_{x}\,,$
$\displaystyle\bm{\varepsilon}_{r}\sim\mathcal{N}(\bm{0},\tau_{\varepsilon_{r}}\bm{D}_{\varepsilon_{r}})\
.$ (9)

Note that the imputation model is defined for $\bm{r}$, as this is the latent
variable entering the classical ME model, whereas the model for the latent
”correct” variable $\bm{x}$ is given through the Berkson ME model. To better
understand the role of the latent variable $\bm{r}$, it may help to imagine
what would happen if one of the error terms was ”turned off” in the model. If
we have $\bm{u}_{b}=\bm{0}$, then (7) simply says $\bm{x}=\bm{r}$, and we are
left with the classical error model we introduced in Section 2.1. On the other
hand, if there is no classical error ($\bm{u}_{c}=\bm{0}$) then
$\bm{w}=\bm{r}$, and we are only left with the Berkson error model. In that
case, the coefficients of the imputation model will be estimated, but it will
not feed back into the model, so in practice it does not affect $\bm{x}$. If
we have Berkson error and missing data, we need to keep the classical error
model in order to to impute the missing values, but we scale the classical ME
precision to a large value for the values that are observed to indicate that
there is no classical ME. The handling of the missingness is then done just as
described in Section 2.2.

###  2.5 Possible extensions

The model described in Section 2.4 can in the most general case account for
classical ME, Berkson ME, and observations missing at random in the same
variable, with the regression model of interest being any generalized linear
model that is compatible with the assumptions INLA takes. There are some
straightforward extensions to this – if there is ME or missingness in multiple
covariates, then separate error and imputation models need to be formulated
for each such covariate. If the error in one covariate is heteroscedastic,
then the elements of $\bm{D}_{u_{c}}$ can correspondingly be adjusted to
reflect the differing precisions for the different observations. Biased errors
can be specified in the error model by introducing an offset. Due to the
modular structure of Bayesian hierarchical models, the components may be
combined with any other model formulations that are possible within a Bayesian
framework, for example random effects or spatial effects (see _e.g._ ,
Gryparis et al.,, 2008).

##  3 Technical considerations for missing and mismeasured covariates in INLA

A requirement to use INLA is that the model can be expressed as a latent
Gaussian Markov random field, which makes INLA slightly more restrictive than
standard sampling approaches. There are also some practical issues that arise
when fitting the particular model described in Section 2.4, which we will now
address.

###  3.1 How are missing values treated in INLA?

The technical issue of how missing values are handled in INLA is not addressed
in the original INLA publication (Rue et al.,, 2009), but the details were
later discussed, for example on the FAQ website https://www.r-inla.org/faq,
question 7: ”How does inla() deal with ’NA’s in the data-argument?”. First of
all, it is crucial to discriminate between missing values in the response
versus missingness in the covariates. If a value is missing in the response,
then that data point gives no contribution to the likelihood, while the
predictive distributions for these missing values in the response will be
computed automatically (Gómez-Rubio,, 2020, Chapter 12.3). If a value in a
covariate is missing, on the other hand, that entry is re-coded to be $0$
internally. This way of handling missing values in the covariates may lead to
biased posteriors, and it is the case we address here.

As mentioned in the introduction, a few publications have tackled missing
values in covariates in INLA previously. Gómez-Rubio and Rue, (2018) combine
Markov chain Monte Carlo with INLA to this end, and Berild et al., (2022) show
how to combine importance sampling with INLA, and both use missing covariate
imputation to illustrate their approaches. However, in both cases the missing
observations are treated as additional parameters of the model, and thus the
model becomes very computationally expensive when the number of missing
observations is large, which makes it of limited practical use. The method
proposed in Gómez-Rubio et al., (2022), on the other hand, is based purely on
INLA, and is much more widely applicable. They also include a logistic model
to describe the actual missingness mechanism, which means that it is possible
to do a sensitivity analysis to examine the influence of the missingness. This
is advisable when one suspects that the missingness mechanism may be MNAR,
that is, when the variable is missing depending on the value of the variable
itself. However, in the cases where we have reason to believe that missingness
is MCAR or MAR it is not necessary to include a model for the missingness
mechanism itself (van Buuren,, 2018), and we can employ a simpler model than
Gómez-Rubio et al., (2022). In that case, we can take advantage of the
interpretation of missing data as an extreme form of ME used here.

###  3.2 Scaling the classical error precision for the case of missing values

In Section 2.2 we described the interpretation of missingness as an extreme
form of ME. Conceptually, scaling the classical error precision to a large
value (_e.g._ , $10^{12}$) would indicate that we have a small error, while
scaling to a very small value (_e.g._ , $10^{-12}$) would correspond to a
missing observation, since we have very large uncertainty. Values for the
precision in between these two extremes would indicate a regular classical ME.
The respective information can thus directly be incorporated into the
corresponding element of the (scaled) diagonal matrix
$\tau_{u_{c}}\bm{D}_{u_{c}}$. However, in practice it is even not necessary to
set the precision to a small value if a covariate entry is missing, because
the values that are completely missing in $\bm{w}$ do not contribute to the
likelihood in the error model component (8) thanks to the way missing values
are handled in INLA (see Section 3.1). We thus only need to scale the
precision according to whether a given observation has ME or not. In the case
when the value is missing, the imputation model will impute the missing value
for $w_{i}$ based on the other covariates through the imputation and error
models in (9) and (8).

###  3.3 Reformulating sub-models that have latent variables as responses

In both missing data and ME problems, the actual covariate $\bm{x}$ is
unobserved and therefore part of the latent field. In a standard Bayesian
hierarchical model, the imputation model imputes the unobserved values based
on other available covariates. However, imputation models like the one
presented in Equation (9) cannot be formulated directly in INLA, since the
framework does not allow latent variables in the response. For the case of ME,
a solution to this problem has however already been suggested by Muff et al.,
(2015): the imputation model in Equation (9) is reformulated with pseudo-
observations $\bm{0}$ as the response, which circumvents the problem as the
imputation model can now be handled as part of the observation model. The same
workaround needs to be used for the specification of a Berkson ME model, where
the unobserved variable $\bm{x}$ would otherwise stand as the response. Our
full model as formulated in INLA then becomes

$\displaystyle\bm{\eta}$
$\displaystyle=\beta_{0}\bm{1}+\beta_{x}\bm{x}+\bm{Z}\bm{\beta}_{z}\ ,$ (10)
$\displaystyle\bm{0}$ $\displaystyle=-\bm{x}+\bm{r}+\bm{u}_{b}\,,$
$\displaystyle\bm{u}_{b}\sim\mathcal{N}(\bm{0},\tau_{u_{b}}\bm{D}_{u_{b}})\ ,$
(11) $\displaystyle\bm{w}$ $\displaystyle=\bm{r}+\bm{u}_{c}\,,$
$\displaystyle\bm{u}_{c}\sim\mathcal{N}(\bm{0},\tau_{u_{c}}\bm{D}_{u_{c}})\ ,$
(12) $\displaystyle\bm{0}$
$\displaystyle=-\bm{r}+\alpha_{0}\bm{1}+\widetilde{\bm{Z}}\bm{\alpha}_{z}+\bm{\varepsilon}_{x}\
,$
$\displaystyle\bm{\varepsilon}_{x}\sim\mathcal{N}(\bm{0},\tau_{\varepsilon_{x}}\bm{D}_{\varepsilon_{x}})\
.$ (13)

The same approach is directly applicable to the case with missing data.
Another issue that arises is the estimation of the coefficient $\beta_{x}$,
which is not straightforward since $\bm{x}$ is unobserved. In fact, although
both $\bm{x}$ and $\beta_{x}$ are assumed to be Gaussian distributed, their
product is not, while INLA requires that the latent field is Gaussian. Again,
a solution to this was presented in Muff et al., (2015) for the ME case, and
thus the approach is directly applicable to the case with missing data as
well. In brief, instead of estimating $\beta_{x}$ as a part of the latent
field, we let it be a hyperparameter, which ensures that the latent field is
still Gaussian, conditioned on $\beta_{x}$.

##  4 Applications

In this section we examine two real-world data sets and carry out a simulation
study. Code for all three examples is available in the Supplementary Material,
where each example is presented in separate documents with extensive
explanations.

###  4.1 Missing covariate imputation: a linear model for cholesterol

In order to illustrate how to model missing covariate imputation in INLA using
the framework of a classical ME model, we consider the common case where only
missing data is present, but no ME. We consider the `nhanes2` (National Health
and Nutrition Examination Survey) dataset from the `mice` R package (van
Buuren and Groothuis-Oudshoorn,, 2011), which has often been used to evaluate
missing data imputation methods, (see, _e.g._ , van Buuren and Groothuis-
Oudshoorn,, 2011; Gómez-Rubio and Rue,, 2018; Gómez-Rubio et al.,, 2022;
Berild et al.,, 2022). The data contains 25 observations and four variables:
age group (which is categorical with three levels), BMI, hypertension and
cholesterol. Note that this data set is not large enough to make inference.
The example is provided as an illustration for how to implement the method for
missing covariate data, since the data set is commonly used for this purpose.
Here, we fit a linear regression model with cholesterol (chl) as the response,
and age group and BMI as covariates. The response cholesterol and the
covariate BMI have 10 and 9 missing observations, respectively, but since the
missingness in the response is automatically handled by INLA (see Section
3.1), the problem reduces to imputing the missing values in the BMI variable.
We specify the joint Bayesian model with a regression model for cholesterol
(14), the error model (15) and the imputation model (16), where the latter
describes how BMI depends on age category:

$\displaystyle\mathit{chl}_{i}$
$\displaystyle=\beta_{0}+\beta_{1}\mathit{age}^{(40-59)}_{i}+\beta_{2}\mathit{age}^{(60-99)}_{i}+\beta_{3}\mathit{bmi}_{i}+\varepsilon_{i}\
,\quad$ $\displaystyle\varepsilon_{i}\sim\mathcal{N}(0,\tau_{y})\ ,$ (14)
$\displaystyle\mathit{bmi}_{i}^{obs}$ $\displaystyle=\mathit{bmi}_{i}+u_{i}\
,\quad$ $\displaystyle u_{i}\sim\mathcal{N}(0,\tau_{u}\bm{D}_{ii})\ ,$ (15)
$\displaystyle\mathit{bmi}_{i}$
$\displaystyle=\alpha_{0}+\alpha_{1}\mathit{age}^{(40-59)}_{i}+\alpha_{2}\mathit{age}^{(60-99)}_{i}+\varepsilon^{(x)}_{i}\
,\quad$ $\displaystyle\varepsilon^{(x)}_{i}\sim\mathcal{N}(0,\tau_{x})\ .$
(16)

Here, $\mathit{bmi}^{obs}_{1},\dots\mathit{bmi}^{obs}_{n}$ are our
observations, which contain missing values, whereas
$\mathit{bmi}_{1},\dots\mathit{bmi}_{n}$ are the latent variables. Age is
dummy-coded for the three age-groups. The coefficients in Model (14) and Model
(16) are all given $\mathcal{N}(0,10^{-6})$ priors. The priors for the
precisions of the residuals $\bm{\varepsilon}$ and $\bm{\varepsilon}^{(x)}$
are $\tau_{y}\sim\mathsf{Gamma}(2,846.8)$ and
$\tau_{x}\sim\mathsf{Gamma}(2,16.8)$, chosen based on information from the
data (this is described in detail in the Supplementary Material), while the
precision for the error term, $\tau_{u}$, is fixed to 1. Since we have no ME
in this example, and since the missing values for $\mathit{bmi}^{obs}_{i}$
(which is the response of the ME model (15)) do not contribute to the
likelihood, we can actually set _all_ the diagonal elements of $\bm{D}$ to
$10^{12}$.

For comparison, we fit two additional models, one complete case model that
uses only the observed BMI values, as well as a simpler imputation model,
where the imputation is independent of age (Figure 2). In this simpler model,
the missing body mass values have a simple Gaussian prior centered at the
average of the observed values and variance four times the variance of the
observed values, as in Gómez-Rubio and Rue, (2018). For completeness, we also
report the exact results from Gómez-Rubio and Rue, (2018). Across all
coefficients, the posterior means for the complete case model are consistently
larger than when an imputation model is included. It is to be expected that
uncertainties decrease when we are able to include more observations, however
the very small credible intervals in this case are more likely an
inconsistency due to the very small sample size. Overall, the example
illustrates that our approach implemented in INLA is a valid alternative to
existing tools for missing data imputation.

Figure 2: Posterior means and 95% credible intervals for the parameters in
Section 4.1. The complete case analysis (red) is a model fit using only the
complete cases ($n=16$). The simple imputation model (yellow) uses a Gaussian
prior for the missing observations, and is identical to the model used in
Gómez-Rubio and Rue, (2018) (orange) ($n=25$). The full imputation model
(blue) also depends on the age of the patient for the imputation of the
missing BMI values ($n=25$).

###  4.2 Imprecise and missing systolic blood pressure measurements

We now consider a case of ME and missing values in the covariates of a
survival model, using the framework described in Section 2. The data is
originally from the Third National Health and Nutrition Examination Survey
(NHANES III), but we use the data as pre-processed and provided by Bartlett
and Keogh, (2018). They linked the NHANES III data to data from the US
National Death Index, with information about the mortality status of each
participant in 2011. Following Keogh and Bartlett, (2021), we consider a
Weibull survival model with the continuous covariates systolic blood pressure
(SBP) and age, and binary covariates diabetes status, sex and smoking status.
The response is time until death by cardiovascular disease. Deaths by other
causes are treated as censorings. Measurements of SBP are known to vary
substantially within the same patient, and SBP is therefore measured twice for
some participants, enabling us to estimate the variance of the error. However,
for some participants there is only one measurement, and for others the SBP is
completely missing. Note that the smoking status is also missing for around
half of the participants, but since we only consider continuous ME, all
observations with missingness in the (binary) smoking status are removed for
this illustrative example, resulting in $n=3433$ observations.

To model the survival time we use a Weibull survival model with hazard
$h(t_{i})=rt^{r-1}\lambda_{i}^{r}$, where $r$ is the shape parameter of the
Weibull distribution and $t_{i}$ is the time until death by cardiovascular
disease for individual $i$. The linear predictor $\eta_{i}$ for individual $i$
is linked to the rate parameter by $\log(\lambda_{i})=\eta_{i}$. The joint
model including the ME model for SBP is then

$\displaystyle\log(\lambda_{i})$
$\displaystyle=\beta_{0}+\beta_{sbp}\mathit{SBP}_{i}+\beta_{sex}\mathit{sex}_{i}+\beta_{age}\mathit{age}_{i}+\beta_{smoke}\mathit{smoke}_{i}+\beta_{diabetes}\mathit{diabetes}_{i}\,$
(17) $\displaystyle\mathit{SBP}^{(1)}_{i}$
$\displaystyle=\mathit{SBP}_{i}+u^{(1)}_{i}\,,$ (18)
$\displaystyle\mathit{SBP}^{(2)}_{i}$
$\displaystyle=\mathit{SBP}_{i}+u^{(2)}_{i}\,,$ (19)
$\displaystyle\mathit{SBP}_{i}$
$\displaystyle=\alpha_{0}+\alpha_{sex}\mathit{sex}_{i}+\alpha_{age}\mathit{age}_{i}+\alpha_{smoke}\mathit{smoke}_{i}+\alpha_{diabetes}\mathit{diabetes}_{i}+\varepsilon^{(x)}_{i}\,,$
(20)

for $i=1,\ldots,n$, where
$\bm{SBP}^{(1)}=(\mathit{SBP}^{(1)}_{1},\dots,\mathit{SBP}^{(1)}_{n})^{T}$ and
$\bm{SBP}^{(2)}=(\mathit{SBP}^{(2)}_{1},\dots,\mathit{SBP}^{(2)}_{n})^{T}$ are
the first and second blood pressure measurements, respectively, with
$\bm{SBP}$ being the latent error-free version of the blood pressure, and the
remaining four variables that are assumed to be without error: sex, age,
smoking status and diabetes status. The ME terms $\bm{u}^{(1)}$ and
$\bm{u}^{(2)}$ and the residual term for the imputation model are assumed to
be independently distributed as $N(\bm{0},\tau_{u}\bm{I})$ and
$\bm{\varepsilon}^{(x)}\sim\mathcal{N}(\bm{0},\tau_{x}\bm{I})$, respectively.
For the analysis using INLA to be comparable to that of Keogh and Bartlett,
(2021), we used the same priors: the coefficients of the model of interest
(the betas) are given $\mathcal{N}(0,10^{-6})$ priors, while the coefficients
of the imputation model (the alphas) are given $\mathcal{N}(0,10^{-4})$
priors, $\tau_{u}\sim\mathsf{Gamma}(0.5,0.5)$, and
$\tau_{x}\sim\mathsf{Gamma}(0.5,0.5)$.

Figure 3: Posterior means and 95% credible intervals for the parameters in
Section 4.2. The naive complete case analysis (red) is a model fit using only
the complete cases ($n=2667$), and with no ME adjustment. The ME adjusted
complete case analysis (yellow) is fit on the same complete case data, but
accounting for the ME in the systolic blood pressure. Finally, the ME adjusted
and imputation model (blue) uses all the observations ($n=3433$) and accounts
for ME and missing data. A dotted vertical line is drawn from the posterior
means of the full imputation and ME model to ease comparison.

The model that accounts for ME and missing values ($n=3433$) is compared to a
naive complete case analysis ($n=2667$) (Figure 3), as well as to the same ME
model but with the complete case data, in order to study the effect of
accounting for the ME but not the missing observations. The naive complete
case analysis is only applied to the individuals that had no missing values in
$\bm{SBP}^{(1)}$, and treats $\bm{SBP}^{(1)}$ as the correct value for
$\bm{SBP}$. We see that the estimated posterior mean for $\beta_{sbp}$
increases in the ME models, which indicates that there may be some attenuation
due to the ME. Comparing the ME adjusted complete case model with the model
that also imputes the missing values, the coefficients for age and diabetes
increase a bit, while the one for SBP shows only a very small increase. It is
notable that in this case the missingness seems to be influencing other
coefficient estimates than the coefficient of the covariate that actually has
the missingness. This may often be the case for both ME and missing data, when
the covariate with the ME and missingness is correlated with other covariates.
Finally, we compared the runtime of the model fit by Keogh and Bartlett,
(2021) in JAGS to the runtime of the model implemented in R-INLA. The model
fit in JAGS took 31 minutes (10000 iterations were used, and the first 5000
were discarded), while INLA took 20 seconds (both were run on a 2 GHz Quad-
Core Intel Core i5 processor). This improvement shows the significant
practical advantage of being able to fit these complex models in R-INLA.

###  4.3 Simulation study: Berkson and classical measurement error alongside
missing data

For this example, we simulated a linear regression model with a mismeasured
covariate $\bm{x}$, observed as $\bm{w}$, as well as an error-free covariate
$\bm{z}$. For the simulation, we first generated the error free covariate as

$\bm{z}\sim\mathcal{N}(\bm{0},\bm{I})\ .$ (21)

Keeping with the notation used in Section 2.4, the latent variable with
Berkson error was then generated as

$\bm{r}=\bm{1}+2\bm{z}+\bm{\varepsilon}_{r}\
,\qquad\bm{\varepsilon}_{r}\sim\mathcal{N}(\bm{0},\bm{I})\ ,$ (22)

and the correct value for the variable was generated by adding further
variation to $\bm{r}$ by

$\bm{x}=\bm{r}+\bm{u}_{b}\ ,\qquad\bm{u}_{b}\sim\mathcal{N}(\bm{0},\bm{I})\ .$
(23)

The classical ME was then added to the variable with Berkson error to obtain
the variable that is actually observed, $\bm{w}$:

$\bm{w}=\bm{r}+\bm{u}_{c}\ ,\qquad\bm{u}_{c}\sim\mathcal{N}(\bm{0},\bm{I})\ .$
(24)

Next, the probability for an observation to be missing was assumed to be MAR,
depending only on $\bm{z}$ as

$\text{P}(w_{i}=\texttt{NA})=\frac{\exp(-1.5+0.5z_{i})}{1+\exp(-1.5+0.5z_{i})}\
,$ (25)

and the response $\bm{y}$ is generated from the correct covariate $\bm{x}$ as

$\bm{y}=\beta_{0}\bm{1}+\beta_{x}\bm{x}+\beta_{z}\bm{z}+\bm{\varepsilon}\
,\quad\bm{\varepsilon}\sim\mathcal{N}(\bm{0},\tau_{y}\bm{I})\ ,$ (26)

with $(\beta_{0},\beta_{x},\beta_{z})=(1,2,2)$. The simulation was repeated
100 times, with $n=1000$ observations per iteration. For each data set we fit
three different models: A complete case model without error correction using
$\bm{w}$ as covariate, a model that accounts for both ME (classical and
Berkson) and missing data (and thus reflects the data generating process), and
a best-case model that is a regression of $\bm{y}$ on $\bm{x}$ and $\bm{z}$,
which shows how this basic regression model would perform if we did not have
any ME or missing data. As for the other examples, complete code to reproduce
this simulation is available in the Supplementary Material. We also provide an
implementation using inlabru rather than R-INLA, since some may find the
syntax of inlabru more user-friendly.

The posterior means of the coefficients in the ME model are very close to
those actually used for the data generation (Figure 4), and clearly improved
compared to the naive model, for which the posterior mean for $\beta_{x}$ is
over-estimated, while the posterior mean for $\beta_{z}$ is under-estimated.
Although $\beta_{z}$ is observed without error, the errors in $\bm{x}$ will
also greatly bias the estimate for $\beta_{z}$, since $\bm{z}$ and $\bm{x}$
are correlated.

Figure 4: Sample means and 95% quantile intervals of the posterior means from
fitting three different models to 100 different simulated data sets. The dots
are the posterior means for each of the runs, and are jittered along the
$y$-axis to reduce overplotting. The first model (red) is the naive complete
case analysis, simply fitting a regression model using $\bm{w}$ instead of
$\bm{x}$, the second model (yellow) adjusts for ME and missing data, and the
third model (blue) is the model using the correct value of $\bm{x}$. A dotted
vertical line is drawn from the posterior means of the best case model to ease
comparison.

##  5 Discussion

In this paper we have addressed how existing ME modelling approaches for INLA
can be used directly to impute missing covariates. By using a model motivated
by a ME problem, we are able to work around the practical issue with covariate
imputation in INLA. The resulting model formulation provides a way to deal
with missingness in a covariate, but it also gives a conceptually simple way
to handle cases where ME and missingness both occur in the same variable. In
addition to classical ME and missing data, we show how the model can be
generalized to account for Berkson ME at the same time. We show that the
computational efficiency of INLA allows us to fit the models significantly
faster than previously possible in sampling-based software for Bayesian
inference.

To reach our goal, we are combining two previously disconnected lines of
research in a synergistic way: the first line of research is treating missing
data as an extreme case of ME, and the second is doing missing data imputation
in INLA. A few publications cited here have used a ME model for missing data
imputation in sampling-based inference, but the idea has until now not been
incorporated in INLA. Simultaneously, there have been a few different
approaches to missing data imputation in INLA, but these have typically not
been feasible if the number of missing observations is large. By using the ME
model to enable imputation in INLA, we have addressed the computational
challenge that was previously a limiting factor. Moreover, even though we here
only considered examples with a single continuous covariate, the model can
easily be expanded to accommodate multiple covariates with missing or
mismeasured variables. On the other hand, it is not straightforward to account
for missingness or ME in categorical covariates, since in INLA the latent
variables must be assumed to be Gaussian. Modelling the misclassification and
missingness of a categorical variable might be achieved by combining this
method with MCMC within INLA.

In summary, we have shown how a joint Bayesian model for ME and missingness
allows us to account for both classical and Berkson error in continuous
covariates of regression models, as well as providing a new way to do
covariate imputation in INLA. A Bayesian analysis gives us the flexibility we
need to reflect the error-generating process directly in the model, as well as
the option to include prior knowledge about the error variance. Our approach
to handle missing data and ME in INLA enables significantly faster inference
compared to sampling-based methods, which makes these complex Bayesian models
for ME and missing data feasible for a wide range of practical applications.

## Data availability statement

The data used in this paper is available in the Supplementary Material, along
with the R-code to reproduce the results. The Supplementary Material can be
found at https://emmaskarstein.github.io/Missing-data-and-measurement-error/.

## Acknowledgements

The authors thank Malgorzata Roos, Ekaterina Poliakova and Christoffer Wold
Andersen for insightful comments on an early version of the manuscript.

## Conflict of interest statement

The authors have declared no conflict of interest.

## References

  * Bartlett and Keogh, (2018) Bartlett, J. W. and Keogh, R. H. (2018).  Bayesian correction for covariate measurement error: A frequentist evaluation and comparison with regression calibration.  Statistical Methods in Medical Research, 27(6):1695–1708.  PMID: 27647812.
  * Berild et al., (2022) Berild, M. O., Martino, S., Gómez-Rubio, V., and Rue, H. (2022).  Importance sampling with the integrated nested laplace approximation.  Journal of Computational and Graphical Statistics, 31(4):1225–1237.
  * Berkson, (1950) Berkson, J. (1950).  Are there two regressions?  Journal of the American Statistical Association, 45(250):164–180.
  * Blackwell et al., (2017) Blackwell, M., Honaker, J., and King, G. (2017).  A unified approach to measurement error and missing data: Overview and applications.  Sociological Methods and Research, 46(3):303–341.
  * Buonaccorsi, (2010) Buonaccorsi, J. (2010).  Measurement Error: Models, Methods, and Applications.  Chapman & Hall/CRC Interdisciplinary Statistics. CRC Press.
  * Buonaccorsi and Lin, (2002) Buonaccorsi, J. P. and Lin, C.-D. (2002).  Berkson measurement error in designed repeated measures studies with random coefficients.  Journal of Statistical Planning and Inference, 104(1):53–72.
  * Carpenter and Smuk, (2021) Carpenter, J. R. and Smuk, M. (2021).  Missing data: A statistical framework for practice.  Biometrical Journal, 63(5):915–947.
  * Carroll et al., (2006) Carroll, R., Ruppert, D., Stefanski, L., and Crainiceanu, C. (2006).  Measurement Error in Nonlinear Models: A Modern Perspective, Second Edition.  Chapman & Hall/CRC Monographs on Statistics & Applied Probability. CRC Press.
  * Cole et al., (2006) Cole, S. R., Chu, H., and Greenland, S. (2006).  Multiple-imputation for measurement-error correction.  International Journal of Epidemiology, 35(4):1074–1081.
  * Deffner et al., (2018) Deffner, V., Küchenhoff, H., Breitner, S., Schneider, A., Cyrys, J., and Peters, A. (2018).  Mixtures of Berkson and classical covariate measurement error in the linear mixed model: Bias analysis and application to a study on ultrafine particles.  Biometrical Journal, 60(3):480–497.
  * Dellaportas and Stephens, (1995) Dellaportas, P. and Stephens, D. A. (1995).  Bayesian analysis of errors-in-variables regression models.  Biometrics, 51(3):1085–1095.
  * Erler et al., (2016) Erler, N. S., Rizopoulos, D., van Rosmalen, J., Jaddoe, V. W. V., Franco, O. H., and Lesaffre, E. M. E. H. (2016).  Dealing with missing covariates in epidemiologic studies: a comparison between multiple imputation and a full Bayesian approach.  Statistics in Medicine, 35(17):2955–2974.
  * Fuller, (1987) Fuller, W. A. (1987).  Measurement Error Models.  John Wiley & Sons, Inc., USA.
  * Goldman et al., (2011) Goldman, G., Mulholland, J., Russell, A., Strickland, M., Klein, M., Waller, L., and Tolbert, P. (2011).  Impact of exposure measurement error in air pollution epidemiology: Effect of error type in time-series studies.  Environ Health, 10:61.
  * Goldstein et al., (2018) Goldstein, H., Browne, W. J., and Charlton, C. (2018).  A Bayesian model for measurement and misclassification errors alongside missing data, with an application to higher education participation in Australia.  Journal of Applied Statistics, 45(5):918–931.
  * Gómez-Rubio, (2020) Gómez-Rubio, V. (2020).  Bayesian Inference with INLA.  CRC Press.
  * Gómez-Rubio et al., (2022) Gómez-Rubio, V., Cameletti, M., and Blangiardo, M. (2022).  Missing data analysis and imputation via latent Gaussian Markov random fields.  SORT-Statistics and Operations Research Transactions, 46(2):217–244.
  * Gómez-Rubio and Rue, (2018) Gómez-Rubio, V. and Rue, H. (2018).  Markov chain Monte Carlo with the integrated nested Laplace approximation.  Statistics and Computing, 28:1033–1051.
  * Gryparis et al., (2008) Gryparis, A., Paciorek, C. J., Zeka, A., Schwartz, J., and Coull, B. A. (2008).  Measurement error caused by spatial misalignment in environmental epidemiology.  Biostatistics, 10(2):258–274.
  * Gustafson, (2003) Gustafson, P. (2003).  Measurement Error and Misclassification in Statistics and Epidemiology: Impacts and Bayesian Adjustments.  Chapman & Hall/CRC Interdisciplinary Statistics. CRC Press.
  * Gustafson, (2005) Gustafson, P. (2005).  On model expansion, model contraction, identifiability and prior information: Two illustrative scenarios involving mismeasured variables.  Statistical Science, 20(2):111–140.
  * Gustafson, (2010) Gustafson, P. (2010).  Bayesian inference for partially identified models.  The International Journal of Biostatistics, 6(2).
  * Heid et al., (2004) Heid, I., Küchenhoff, H., Miles, J., Kreienbrock, L., and Wichmann, H. (2004).  Two dimensions of measurement error: Classical and Berkson error in residential radon exposure assessment.  Journal of Exposure Analysis and Environmental Epidemiology, 14(5):365–377.
  * Innes et al., (2021) Innes, G. K., Bhondoekhan, F., Lau, B., Gross, A. L., Ng, D. K., and Abraham, A. G. (2021).  The measurement error elephant in the room: Challenges and solutions to measurement error in epidemiology.  Epidemiologic Reviews, 43(1):94–105.
  * Keogh and Bartlett, (2021) Keogh, R. H. and Bartlett, J. W. (2021).  Measurement error as a missing data problem.  In Handbook of Measurement Error Models, Chapman & Hall/CRC Handbooks of Modern Statistical Methods, pages 429–452. CRC Press.  Available as preprint, arXiv:1910.06443.
  * Keogh et al., (2020) Keogh, R. H., Shaw, P. A., Gustafson, P., Carroll, R. J., Deffner, V., Dodd, K. W., Küchenhoff, H., Tooze, J. A., Wallace, M. P., Kipnis, V., and Freedman, L. S. (2020).  Stratos guidance document on measurement error and misclassification of variables in observational epidemiology: Part 1 — Basic theory and simple methods of adjustment.  Statistics in Medicine, 39(16):2197–2231.
  * Little and Rubin, (1987) Little, R. and Rubin, D. (1987).  Statistical Analysis With Missing Data.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc.
  * Ma and Chen, (2018) Ma, Z. and Chen, G. (2018).  Bayesian methods for dealing with missing data problems.  Journal of the Korean Statistical Society, 47(3):297–313.
  * Marquardt, (2020) Marquardt, K. L. (2020).  How and how much does expert error matter? Implications for quantitative peace research.  Journal of Peace Research, 57(6):692–700.
  * Martino and Riebler, (2020) Martino, S. and Riebler, A. (2020).  Integrated Nested Laplace Approximations (INLA), pages 1–19.  John Wiley & Sons, Ltd.
  * Muff et al., (2017) Muff, S., Ott, M., Braun, J., and Held, L. (2017).  Bayesian two-component measurement error modelling for survival analysis using inla–A case study on cardiovascular disease mortality in Switzerland.  Computational Statistics & Data Analysis, 113:177–193.
  * Muff et al., (2015) Muff, S., Riebler, A., Held, L., Rue, H., and Saner, P. (2015).  Bayesian analysis of measurement error models using integrated nested Laplace approximations.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 64(2):231–252.
  * Noghrehchi et al., (2020) Noghrehchi, F., Stoklosa, J., and Penev, S. (2020).  Multiple imputation and functional methods in the presence of measurement error and missingness in explanatory variables.  Computational Statistics, 35(3):1291–1317.
  * Ponzi et al., (2018) Ponzi, E., Keller, L. F., Bonnet, T., and Muff, S. (2018).  Heritability, selection, and the response to selection in the presence of phenotypic measurement error: Effects, cures, and the role of repeated measurements.  Evolution, 72(10):1992–2004.
  * Reeves et al., (1998) Reeves, G. K., Cox, D. R., Darby, S. C., and Whitley, E. (1998).  Some aspects of measurement error in explanatory variables for continuous and binary regression models.  Statistics in Medicine, 17(19):2157–2177.
  * Richardson and Gilks, (1993) Richardson, S. and Gilks, W. R. (1993).  Conditional independence models for epidemiological studies with covariate measurement error.  Statistics in Medicine, 12(18):1703–1722.
  * Rue et al., (2009) Rue, H., Martino, S., and Chopin, N. (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2):319–392.
  * van Buuren, (2018) van Buuren, S. (2018).  Flexible Imputation of Missing Data. Second Edition. CRC Press, Boca Raton, FL.
  * van Buuren and Groothuis-Oudshoorn, (2011) van Buuren, S. and Groothuis-Oudshoorn, K. (2011).  mice: Multivariate imputation by chained equations in R.  Journal of Statistical Software, 45(3):1–67.
  * van Smeden et al., (2019) van Smeden, M., Lash, T. L., and Groenwold, R. H. H. (2019).  Reflection on modern methods: five myths about measurement error in epidemiological research.  International Journal of Epidemiology, 49:338–347.
  * van Smeden et al., (2021) van Smeden, M., Penning de Vries, B. B., Nab, L., and Groenwold, R. H. (2021).  Approaches to addressing missing values, measurement error, and confounding in epidemiologic studies.  Journal of Clinical Epidemiology, 131:89–100.
  * Yi, (2017) Yi, G. (2017).  Statistical Analysis with Measurement Error or Misclassification: Strategy, Method and Application.  Springer Series in Statistics. Springer New York.

###### List of Figures

  1. 1 The scale of classical measurement error variance. The figure is adapted with permission from Blackwell et al., (2017) (copyright 2015 by the authors).
  2. 2 Posterior means and 95% credible intervals for the parameters in Section 4.1. The complete case analysis (red) is a model fit using only the complete cases ($n=16$). The simple imputation model (yellow) uses a Gaussian prior for the missing observations, and is identical to the model used in Gómez-Rubio and Rue, (2018) (orange) ($n=25$). The full imputation model (blue) also depends on the age of the patient for the imputation of the missing BMI values ($n=25$).
  3. 3 Posterior means and 95% credible intervals for the parameters in Section 4.2. The naive complete case analysis (red) is a model fit using only the complete cases ($n=2667$), and with no ME adjustment. The ME adjusted complete case analysis (yellow) is fit on the same complete case data, but accounting for the ME in the systolic blood pressure. Finally, the ME adjusted and imputation model (blue) uses all the observations ($n=3433$) and accounts for ME and missing data. A dotted vertical line is drawn from the posterior means of the full imputation and ME model to ease comparison.
  4. 4 Sample means and 95% quantile intervals of the posterior means from fitting three different models to 100 different simulated data sets. The dots are the posterior means for each of the runs, and are jittered along the $y$-axis to reduce overplotting. The first model (red) is the naive complete case analysis, simply fitting a regression model using $\bm{w}$ instead of $\bm{x}$, the second model (yellow) adjusts for ME and missing data, and the third model (blue) is the model using the correct value of $\bm{x}$. A dotted vertical line is drawn from the posterior means of the best case model to ease comparison.