DP2020-06

Semiparametric Bayesian
Instrumental Variables
Estimation for Nonignorable
Missing Instruments
Ryo KATO
Takahiro HOSHINO

February 5, 2020

Semiparametric Bayesian Instrumental
Variables Estimation for Nonignorable
Missing Instruments
Ryo Kato
Research Institute for Economics and Business Administration, Kobe
University, 2-1 Rokkodai-cho, Nada-ku, Kobe, Japan
and
Takahiro Hoshinoy
Department of Economics, Keio University, 2-15-45 Mita, Minato-ku,
Tokyo, Japan / RIKEN Center for Advanced Intelligence Project, 1-4-1
Nihonbashi, Chuo-ku, Tokyo, Japan
Summary.This paper considers the case where instrumental variable (IV)are available to infer the e¤ect of interested variable to the outcome (orthe causal e¤ect), but some components of IV are missing with the missing
mechanism of not missing at random (NMAR). Although NMAR requires theanalysis to prespecify the missing mechanism, it is unknown for us and what
is worse, it is generally not identi…ed. We use the IV distribution of original
population as an auxiliary information, and show that missing mechanism
can be represented as identi…able nonparametric generalized additive model.
We also introduce MCMC algorithm that impute the missing values andsimultaneously estimate parameters of interested.
Key words: Instrumental variable; missing not at random; auxiliary information.
1. Introduction
When it is not feasible to conduct randomized controlled trials (RCT) or
quasi-randomized experiments, the IV approach can be a very useful tool toinfer the causal e¤ect if it is possible to …nd su¢ cient IVs. Since they canproperly eliminate the confoundings caused by unobserved factors, IV models
are developed and applied in many empirical economics researches, whereunobserved confoundings are ubiquitous. Therefore, introducing su¢ cientyemail: kato ryo@rieb.kobe-u.ac.jp
email: bayesian@jasmine.ocn.ne.jp
1

IVs can in itself be great invention, and many researchers are trying to …ndthem.
Despite the desperate e¤orts, IVs tend to be missing. For example, information of twin are often used as an instrument, but it is only observed
for sub-samples. Aaslund and Gronquist (2010) used twin birth as an instrument to survey the e¤ect of family size on the quality of children. In thiscase, a twin birth can be observed for families with twins. Of course, thecomplete case analysis seems to result in biased results, and as an alternative
approach, restricting the sample to families with more than two children losethe e¢ ciency in sample size. Variables on children, such as child BMI, arealso frequently employed as instruments. However, since child information
comes from di¤erent sources than the endogenous parents’information (e.g.
BMI), there is a tendency for the former to be missing. Real data analysis
section is an example wherein instruments are missing for many observations
since they are sourced from other surveys.
In addition to economics, missing IV is a common problem in other …elds.
Mendelian randomization uses genotype information as an instrumental variable to infer the causal e¤ect of a biomarker to a disease. Since the appropriate genetic variant is independent of the confounders of the intermediate
phenotype-outcome association and can a¤ect the outcome only through thecausal intermediate phenotype as long as it is related to the intermediate
phenotype, it has recently been applied in economics as well as in biostatistics. In general, as genetic variants explain only a small portion of the endogenous population, Mendelian randomization requires large sample sizes(Smith, 2006) to satisfy enough causal associations. However, Mendelianrandomization datasets are often missing (Palmer et al., 2012) and a large
enough sample size cannot be guaranteed.
These example of missing IV shows us that complete case analysis, namely,
only those samples where all the instruments are observed, are used in theanalysis, thus resulting in biased results and wrong decision-making. In thispaper, we develop a semiparametric method to impute the missing portionof IV and simultaneously infer the causal e¤ect. As the most general case,we consider the case with not missing at random (NMAR). Therefore, missingness of instruments remains associated with the missing instruments evenafter controlling for other observed variables. In the NMAR case, the interested regression model cannot be identi…ed without an additional assumption (Little and Rubin, 2002) . An example of such an assumption is strong
parametric assumption on the regression and missing mechanism (Kott andChang, 2010) . However, Miao et al. (2016) showed that probit speci…cation2

on the missing mechanism can identify normal and normal-mixture models
while logit speci…cation can less identify them.
Although there exists a lot of literature on IV, the model for missing instrument is scant (Kennedy and Small, 2017) except for (Burgess et al., 2011;
Mogstad and Wiswall, 2012; Chaudhuri and Guilkey, 2016). Burgess et al.
(2011) considered missing instruments for Mendelian randomization, and proposed a Bayesian multiple imputation method. However, they considered thecase of missing at random. Mogstad and Wiswall (2012) proposed an IV estimator for partially missing instrument, but they also assume missing at
random. Chaudhuri and Guilkey (2016) developed semiparametric e¢ cient
GMM method for missing IV, assuming IV missing at random. As anotherrelated work Ertefaie et al. (2017) considered the case of NMAR with observed confonders, but they assume IVs are completely observed. Therefore,
there exist no literature which considered missing IV with NMAR.
We take an assumption that the IV distribution of the original population
is available as an auxiliary information. In many cases, the population-level
information is available from other data sources. Government statistics issuch an example and some researches utilize this auxiliary information toestimate individual-level causality. Imbens and Lancaster (1994) and Hellerstein and Imbens (1999) incorporated population-level information as moment conditions to infer individual-level models using the generalized methodof moments (GMM). Another instance where population-level information isused is the empirical likelihood estimation (Qin, 2000; Chaudhuri et al.,2008). Such approaches are also applied to the missing data issues. Nevo
(2003) proposed the propensity score weighting method using the moment
conditions obtained from auxiliary population-level information. Igari andHoshino (2018) introduced the Bayesian method with population-level information that dealt with repeated durations under unobserved missing indicators.
Although prior works incorporating population-level information to deal
with missing variables use moment conditions, our proposed method uses
probability distribution of population as auxiliary information since the moment conditions have less information than original distribution. Under thecondition that the original population distribution of the missing IV is known,
followed by the theorem in Hirano et al. (2001), we show that the missing
mechanism is nonparametrically identi…ed with generalized additive model,
and the substantive IV regression models are also identi…ed. Figure 1 illustrates the model considered in this paper. In general, since fully nonparametric missing mechanism are not identi…ed, parametric missing mechanism are3

frequently assumed (Kott and Chang, 2010). However, misspeci…cation of
missing mechanisms results in severely biased estimates (Kim and Yu, 2011).
Kim and Yu (2011) developed a semiparametric missing mechanism approach
which incorporates nonparametric speci…cations on observed variables but
not on unobserved variables. However, their method, as well as other prior
works, cannot identify the nonparametric part of unobserved variables. On
the other hand, we assume the availability of the information of the original population distribution of missing IV so that our proposed method canspecify the fully nonparametric missing mechanisms on observed and unobserved variables. Furthermore, our missing mechanism can incorporate crossterms of observed and unobserved variables, which cannot be identi…ed by
the existing methods.

Fig.1 The model considered in this paper

2. The challenges of handling missing IV
2.1 Setup and notations
Before discussing the methodologies for handling datasets with missing
IVs, we introduce the notations and setup.
IV approach aims to overcome the unmeasured confounders problem in order to make causal inference. We consider potential outcome approach, thende…ne the causal e¤ect. Let Y be the outcome variable, X be a treatment(X )variable or endogenous variable. Let Yi
denote the potential outcomethat would be observed for unit i if the individual were to have the treatment level of X . For each unit, only one possible realization of Yiobs = Yi and
Xiobs = X is observed. In this potential outcome approach, although each4

unit has potential outocomes corresponding to the possible level of X, only
(X )one potential outcome Yi = Yi
can be observed. Let Zi be a instrumental
variable for unit i.
Let (Yi ; Xi ; Zi ) ; (i = 1; :::; N ) be an i.i.d. sample. We consider IV regression model with one structural equation or “second-stage” and one reduced
form equation or “…rst-stage”in the following form
Yi = M (Xi ; 0 ) + 1;i:
Xi = f (Zi ; ) + 2;i(1)where M (Xi ; 0 ) has only a …nite-dimensional unknown parameter 0 ; butf (Zi ) are an unknown functions. We consider the situation where 1;i and2;i are correlated, which causes the endogeneity problem. We assume Zi to
be independent of 2;i and E( 1;i Zi ) = 0. We call
We consider the case where X and are scalar, that is M (Xi ; 0 ) = 0 Xi :In this case, our parameter of interest is the causal e¤ect of increasing X by
(X +1)(X )one unit: 0 = Yi
Yi:
Although some recent literature consider the case that the function M is
unknown and they estimate M nonparametrically, we assume M (Xi ; 0 ) =0 Xi because nonparametric estimation of M require us to very strong identi…ability assumptions (Newey and Powell, 2003; Hall and Horowitz, 2005;
Darolles et al., 2011; Liao and Jiang, 2011; Kato, 2013).
We consider the missingness of instrumental variables Z. We assume that
Z for some observations are missing and the other variables, Y and X arecompletely observed. We denote the missing indicator as Ri which takes 1
when the corresponding element of the Zi is observed and 0 otherwise.
We assume the missing mechanism to be nonignorable (NMAR). We consider the case with the missing probability as depending on all the otherobserved variables p (R = 1jY; X; Z) (See also Figure 1).
2.2 Missing instruments and complete case analysis
The simplest and the most applied way to deal with missing IVs arethought to be complete case analysis. We consider the results obtained fromcomplete case analysis using estimating equation. The estimating equationof complete case analysis is
E [I (R = 1) 1 fY

f (Zi ;

0)g]

where 0 is the true value of on the …rst stage, and I (R = 1) = 1 forcomplete case and I (R = 1) = 0, otherwise.

5

If the probability of missing depends only on Z (NMAR), that is, p (R = 1jZ),then the estimates of obtained from complete case analysis results has consistency, as shown from the following estimating equation since E ( 1 jZ) is
assumed to be 0:
E [I (R = 1jZ) 1 fY f (Zi ; 0 ) g]
= EZ (EY [I (R = 1jZ) 1 fY f (Zi ; 0 ) g jZ])
Z
Z
p (R = 1jZ) 1
fY f (Zi ; 0 ) g f (Y jZ) dY f (Z) dZ
=
Z
Y
Z
Z
p (R = 1jZ) 1=
1 f ( 1 jZ) d 1 f (Z) dZ
Z

1

= 0

On the other hand, for example, if the probability of missing only depends
on X (MAR), that is, p (R = 1jX),
E [I (R = 1jX) 1 fY f (Zi ; 0 ) g]
= EX (EY [I (R = 1jX) 1 fY f (Zi ; 0 ) g jX])
Z
Z
p (R = 1jX) 1
fY f (Zi ; 0 ) g f (Y jX) dY f (X) dX
=
X
Y
Z
Z
p (R = 1jZ) 1=
1 f ( 1 jX) d 1 f (Z) dZ
z

1

6= 0

shows complete case analysis results in inconsistent result since E ( 1 jX) 6= 0:
In the setting, since E [ (R = 1jY; X; Z) 1 fY f (z; 0 ) g] 6= 0, thecomplete case analysis results in biased estimates. Therefore, we must specify
the missing mechanism to obtain the consistent estimates of the interested
parameters.
2.3 Missing instruments and existing imputation methods
As stated above, no method has proposed to impute missing instruments
when the pattern of missingness is NMAR. We review some imputation methods which can just imputed missing components, though they are not expected to obtain consistent estimates.
2.3.1 Multiple imputation by chained equation Multiple imputation bychained equation (MICE) (van Buuren, 2007) is one of the most applied imputation method since researchers can avoid the di¢ culty in specifying the6distributional assumptions. MICE speci…es a multivariate distribution by asequence of univariate regressions for each missing variable. Moreover, it caneasily implement using existing several software packages, such as the “mice”
package in R, “proc mi” with the FCS option in SAS, and “mi impute” in
STATA.
However, Liu et al. (2014) showed that the imputations with MICE do
not guarantee the asymptotic distributions to be consistent with the existing Bayesian joint model multiple imputation estimator when the family of
conditional models and their joint distributions are not “compatible”. The
violation of compatible is ubiquitous and MICE is shown to result in severely
biased estimated in many cases using simulation studies (Kato and Hoshino,2019). What is worse, since it do not specify the pattern of missingness whichis required to NMAR, MICE result in biased in the case. Then the situationconsidered here are not appropriate for MICE. In fact, the simulation study
conducted in this paper shows considerable biased results when missing IV
is imputed using MICE.
2.3.2 Imputation by machine learning method Recently, machine leaning
based method has been applied to impute the missing values. Especially,
Stekhoven and Buhlmann (2012) proposed the missForest algorithm, in which
the missing values are imputed by the predictors from the random forest.
Since missForest predicts missing values based on random forest, it accommodate interactions and nonlinearities. Therefore, we do not have to specify aparticular regression model for imputation. missForest is often implemented
to the real world dataset since they are reported to provide lower imputation
errors than the FCS method (Waljee et al., 2013; Liao et al., 2014).
However, it is reported that missForest approach results in biased estimates in some cases when the number of covariates are small (Kato andHoshino, 2019). Moreover, even if the algorithm correctly estimates themean, it always underestimates the variance of the estimates and resultsin poor CIs (Shah et al., 2014) since they do not rely on the probabilistic
model. Therefore, it is not suitable for the case applying IV method where
the results of statistical signi…cance (or hypothesis testing) are crucial (suchas economics, medical, or epidemiological research). In fact, the simulation
study conducted in this paper shows considerable biased results when missing
IV is imputed using missForest.

7

3.

Semiparametric Bayes model for instrumental variable with
nonignorable missing
In this section, we proposed semiparametric Bayesian model for instrumental variable with nonignorable missing that overcome the shortcomings ofexisting imputation methods.
3.1 Identi…cation
We are interested in recovering the joint distribution of (Y; X; Z) or possibly the conditional distribution (Y; XjZ) with missing problem of Z. Since
p (Y; X; Z) =f (Y; X; ZjR = 1) p (R = 1)
p (R = 1jY; X; Z)

and f (Y; X; ZjR = 1) and p (R = 1) can be directly estimated, identi…cationof p (Y; X; Z) depends on the identi…cation of p (R = 1jY; X; Z) : Because
the missing probability of Z depends on Z itself: p (R = 1jY; X; Z), this is
NMAR case.
Generally, the identi…cation of missing mechanism is di¢ cult under NMAR
(Little and Rubin, 2002). However, by considering following assumptions, we
can identify the model.
Theorem 1:
Assume that the following conditions hold.
(A1) the support of f (y; x; z) coincides with f (y; x; zjR = 1) ;(A2) f (y; x; z) and f (y; x; zjR = 1) are square integrable;(A3) the marginal distribution of Z, p(Z) is known.
Then the missing mechanism is identi…ed through the form ofp (Ri = 0jY = y; X = x; Z = z) = h (k0 + k1 (y; x) + k2 (z))(2)

where h is a known function which is di¤erentiable, strictly increasing withlima! 1 g(a) = 0 and lima!1 g(a) = 1, and k1 ( ) ; k2 ( ) ; k3 ( ) are unique
set of functions subject to normalization k1 (0) = k2 (0) = k3 (0) = 0:
Proof:The proof for theorem 1 is very similar to that of Hirano et al. (2001). By
regarding the set (Y; X) as time-variant completely observed variables, andZ as time-variant incompletely observed variable as in the variable de…nitionof Hirano et al. (2001), we can straightforwardly obtain the desired result.
We can also include measured confounders (exogenous variables) in thesecond stage regression model. We denote Wi as a k-dimensional vector of8measured confounders which contains Zi . We consider following system of
equations
Yi = M (Xi ; 0 ) + g (Wi ; ) + 1;i:
Xi = f (Zi ; ) + 2;i(3)where g ( ) is an unknown functions. With this speci…cation, we can identify
the missing mechanism with some assumptions.
Theorem 2
Assume that the following conditions hold.
(A10 ) the support of f (y; x; zjw) coincides with f (y; x; zjw; R = 1) ;(A20 ) f (y; x; zjw) and f (y; x; z; wjR = 1) are square integrable almost
surely with respect to w;(A30 ) the marginal distribution of Z, p(Z) is known.
Then the missing mechanism is identi…ed through the form ofp (Ri = 0jY = y; X = x; Z = z; W = w) = h(k 00 (w) +k 01 (y; x; w) +k 02 (z; w) )(4)
where h is a known function which is di¤erentiable, strictly increasing withlima! 1 g(a) = 0 and lima!1 g(a) = 1, and k1 ( ) ; k2 ( ) ; k3 ( ) are unique
set of functions subject to normalization k10 (0) = k20 (0) = k30 (0) = 0:
Proof:The proof for theorem 2 is the same as that of Hirano et al. (2001).
By regarding the set (Y; X) as time-variant completely observed variables,
Z as time-variant incompletely observed variable, and W as time-invariant
covariates as in the variable de…nition of Hirano et al. (2001), we can straightforwardly obtain the desired result.
The advantages of this missing mechanism speci…cation is: (i) we canspecify nonparametric forms on the observed variables and missing variables(IVs), and (ii) we can consider cross term e¤ect of observed variables andmissing variables. Neither of these formulations are identi…ed by the semiparametric models proposed by Kott and Chang (2010) or Kim and Yu (2011)(Because the method of Kott and Chang (2010) is developed for nonignorable nonresponse and Kim and Yu (2011) is for mean functionals, they are
not applied to IV regression).
3.2 Semiparametric formulation
We formulate the conditional distribution of the missing IV z given all
the other variables as follows.
9p (R = 0jY; X; W; Z) p (ZjY; X; W )R
p (R = 0; ZjY; X; W ) dZz
/ p (R = 0jY; X; W; Z) p (Y; XjW; Z) p (ZjW ) p (Z)p (ZjY; X; W; R = 0) =where p (R = 0jY; X; W; Z) is the missing mechanism, p (Y; XjW; Z) is thesubstantive IV model (the structural equation and the reduced-form equation), p (ZjW ) is conditional instrumental variables distribution, and p (Z)is known and true marginal distribution of Z:
As stated in the identi…cation issue, we can specify nonparametric missing mechanism for p (R = 0jY; X; W; Z) using generalized additive model inequation (4). Needless to say, parametric probit or logistic regression can beprespeci…ed to p (R = 0jY; X; W; Z).
3.2.1 Substantive model and conditional IV distribution We formulate thesubstantive model p (Y; XjZ; W ) = p (Y jX; W ) p (XjZ) in the semiparametric form. As stated, p (XjZ) represents the …rst stage: Xi = f (Zi ; ) + 2;i ,and p (Y jX; W ) represents the second stage: Yi = M (Xi ; 0 )+g (Wi ; )+ 1;iin equation (3), respectively. Therefore, we consider p (Y; XjZ; W ) to be
semiparametric IV regression model. Moreover, we specify the conditional
IV distribution p (ZjW ) to be nonparametric.
We use DPM representation to represent semiparametric IV regression
model. Fortunately, the DPM model can be estimated with a relatively
simpler MCMC algorithm by applying blocked Gibbs sampling (Ishwaran andJames, 2011). Therefore, DPM modeling has been often employed to avoid
distributional assumptions on the parameters within the Bayesian framework(e.g. Hirano, 2002; Xue et al., 2007; Conley et al., 2008; Shahbaba and Neal,2009). The theoretical properties of DPM were investigated in Shen et al.
(2013).
Dirichlet process prior is assumed for a random distribution G, denoted
as G DP ( ; G0 ) is expressed as follows:
G=

1
X

l

l

;

l

G0

l=1

where is the concentration parameter and G0 is the base distribution.
is a point mass at and
Y
(1
=
l
l
h)
h<l

10

where l Be (1; ) :By applying DPM, the resulting multivariate regression function of yi ; xi ,and zi is represented by following mixture model;p (yi ; xi ; zi jwi ) = p (yi ; xi jzi ; wi ) p (zi jwi )X1
fi el ; l N (zi j l wi ;e i jX
=
lN y
l=1

where

yi
xi

ei =
y

fi =
;X

el =

0

xi wi 00 0 zi
1;l
2;l

; el

N (0;

; el = @

l
l

l ):

l)

1

A;

4. MCMC algorithm
Blocked Gibbs sampler (Ishwaran and James, 2011) is applied to the posterior
computation of the DPM parameters. It is proved that the case with …nite
number of classes l can be used to approximate the inference that is basedon in…nite classes with satisfactory accuracy when the maximum number ofclasses L is large enough. Therefore, blocked Gibbs sampler considers thecase for truncation of the number of classes (e.g. L = 20), hence the simpler
posterior computation.
We obtain the detailed posterior computation using the MCMC estimation as follows.
1. Conditional for Ki (i = 1; :::; N )Let Ki be the indicator denoting where case i belongs, and Ki = l if
casePi belongs to class l: To assign samples to each class, generate Ki
by Ll=1 li l ( ), where li is
lN

li = XL

l=1

with

l =

l

Q

h<l (1fi el ;
yei jXlfi el ;
ei jX
lN y

h ).

11

N (zi j l wi ;ll)

N (zi j l wi ;;l)

2. Conditional for (l = 1; :::; L
Simulatel from the following normal distribution.

l

3. Update

1)

l;

N

P

i:Ki l Wil +Nl + 1;

1
Nl + 1;l

Draw l and l from the following multivariate normal and inverted
Wishart distribution.
l jrestl jrestN vec( b );W Tl W ll

IW f0 + N; G0 1 + (Z l
1

1

;

T
l W l ) (Z l
T

lW l)

;
T

T
where b = W Tl W l
W l Z l , W l = wT1 ; :::; wTN , Z l = z1T ; :::; zN,
and W l and Z l denote the subset of W l and Z l whose case i belong
to class l. f0 and G0 1 denotes the parameter of the prior distributionof l ; l IW (f0 ; G0 1 ).

4. Update the missing mechanism
If we assume parametric probit or logistic model to the missing mechanism, we can use Gibbs sampling proposed by Albert and Chib (1993)
for probit and M-H algorithm for the logistic, given the pseudo-complete
dataset. If we assume generalized additive model to the missing mechanism, we can directly employ Bayesian GAM estimation proposed by
Klein et al. (2015) given the pseudo-complete dataset.
5. Update the missing components
Draw the missing component of z i from a density proportional top (zi jyi ; xi ; wi ; r i = 0). Since it is di¢ cult to draw the missing zi , weemploy the Metropolis-Hastings algorithm and use p (zi jwi ) as a proposal density in order to draw a candidate of zi ; zi;c . We accept thecandidates with the following probability:"#p (ri = 1jyi ; xi ; wi ; zi ) p (yi ; xi jzi ; wi ) p(zi jwi )min 1;
p ri = 1jyi ; xi ; wi ; zil 1 p yi ; xi jzil 1 ; wi p(zil 1 jwi )6. Conditional forl12The standardized reduced model is used to draw. Substituting in for x
in the second stage regression, the followings are obtained
xi
wi l

yi

zi
zi

=

1 0
1
l
1
0 1
regression with unit normal errors.

where V ar ( i ) =7. Conditional for

and

l +

i

: Then we drawl using Bayesianl

The structural model parameters and can be simulated by applying
Bayesian linear regression draw with N (0; 1) error termsyi

E [ 1;i j 2;i ]=1;i j 2;ixi+

1;i j 2;iwi

+ ui ;

where=1;i j 2;iui

N (0; 1)

1;i j 2;i

2
12;! i2
11;! i11;! iandE [ 1;i j 2;i ] =12;! i2;i22;! i

When we draw , we consider linear regression model with standard
normal error
yiw =where yiw = (yiWhen we draw
normal error+ ui ;ui

N (0; 1)

1;i j 2;iE [ 1;i j 2;i ]wi l ) =1;i j 2;i:l , we consider linear regression model with standard

yix =where yix = (yi

xi

wi l

+ ui ;ui

N (0; 1)

1;i j 2;iE [ 1;i j 2;i ]xi ) =

13

1;i j 2;i:5. Simulation study
We conduct two simulation studies in order to illustrate the performanceof the proposed method when some components of instruments are missing
with NMAR and the original population information is available. Simulation 1 considers the case where the original population distribution of IV is
log-normally distributed with a linear reduced-form equation and additive
disturbance that follows the Gaussian distribution. The missing mechanism
contains the cross-term of IV and observed variables, which may not be identi…ed by existing methods that do not have population-level information.
Simulation 2 considers the case where the original population distributionof IV is log-normally distributed with the mixture of two classes’ (l = 2)
reduced-form equation. The missing mechanism also contains the cross-termof IV and observed variables. Through the simulation study, 30% of theincomplete covariates are set to be missing. We generate the missing values
based on p (ri = 0jyi ; xi ; wi ; zi ) (see Appendix for detailed missing probability
setting).
In the two simulation studies, we evaluate the …nite sample property ofthe estimators. We consider 1; 000 replications of the dataset and con…rm
if the true value of endogenous variable coe¢ cient = 1 can be recovered.
Throughout the simulations, the sample size for each dataset is set to be
N = 400. We calculate standard performance measures, such as empiricalmean, standard deviation, the coverage of nominal 95% con…dence (or credible) intervals (CI) of the estimate, and the deviation (MSE) from the true
value. Although it could be sometimes inappropriate, we take more practical perspective to compare the Bayesian credible intervals with the classicalcon…dence intervals in Conley et al. (2008).
We evaluate the proposed method by comparing it with other data imputation methods, i.e., MICE-GMM and “pseudo-proposed method” (PsProposed). As described, the MICE approach speci…es a multivariate covariate distribution by a sequence of univariate regressions for each missing
variable. In our simulation settings, where the conditional models (set of
univariate regressions) and their joint distributions are incompatible, MICE
estimates do not guarantee consistency. After creating multiple datasetsfrom MICE, the IV estimates are obtained by GMM (generalized methods of moments). We call this competitor MICE-GMM. We also employed
missForest imputation followed by IV estimation with GMM (MF-GMM).
“Pseudo-proposed method”(Ps-Proposed) is the alternative which does not
use the original population information, but where the missing instruments
are imputed from p (zi jyi ; xi ; wi ; ri = 0). Ps-Proposed might not have identi14

…cation. We also compared this method with IV estimates from other classical methods based on complete case analysis by GMM (CC-GMM), and
OLS estimates which ignore the existence of endogeneity (CC-OLS). As forBayesian methods, burn-in of 5,000 iterations followed by 10,000 iterations
were used for the posterior inference. MICE-GMM and CC-GMM assumes
an ordinary linear regression model for the reduced-form equation, that is,xi = zi + 2;i : MICE-GMM and GMM-CC used the moment condition ofE( 1;i zi ) = 0, therefore GMM based methods do not prespecify the functional
form of the reduced-form equation. We used inverse variance of the moments
as a weighting matrix, namely, the optimal weighting matrix. We calculated
the intervals based on large sample approximations for the classical estimators.
This section presents the results of simulation study 1. The detailed simulation settings and the results of simulation 2 are provided in the Appendix.
Table 1 shows the results of simulation 1, including the empirical mean,
standard deviation, the coverage of nominal 95% con…dence intervals (CI)of the estimate, and the mean squared error (MSE) from the true value of. Figure 2 shows the box plot of the estimates for obtained by 1; 000
replications.

As can be seen from the table and …gure, all estimation methods except
the proposed seem to be biased. Since Ps-Proposed do not have the information of the original distribution of IV, the missing mechanism cannot be
identi…ed and the obtained estimates are upwardly biased. The coverage of
95% credible intervals is also poor at 32.3%. As a result, the quasi-proposed15

yields an MSE over 1,000 times larger compared with the proposed. The
results indicate that population-level information is very useful to NMAR IV
regression case. Another imputation approach is MICE-GMM, which also
shows biased results. The biased results seem to be caused by incompatibility. Therefore, the standard deviation and the MSE are much larger thanthe proposed method.
The complete case analysis-based methods also show poor results. CCGMM estimates have large standard deviation and range in very wide region.
The MSE is about 900 times larger than the proposed method. The complete
case analysis results as biased are indicated in Section 2.2. The coverage of
95% con…dence intervals of CC-OLS is only 1.1%. This is because CC-OLS
ignores the existence of IV and missing values.
As can be seen from comparison with the results from the complete data,
as well as the competing methods, the proposed semiparametric approachshows very good estimates.

Fig.2 The solid horizontal line is the “true”coe¢ cient value = 1. The boxes
span the range from the 25th to the 75th percentiles, and the whiskers extend to anarea no more than 1.5 times the range from the 25th to the 75th percentiles fromthe box. The circles above and below the whiskers represent outliers. Complete,complete case analysis; Ps-Proposed, pseudo-proposed method which does not use
population-level information; MICE-GMM, IV estimates by generalized methodof moments based on multiple imputation by chained equation; MF-GMM, IV
estimates by generalized method of moments based on missForest imputation;

16

CC-GMM, IV estimates by generalized method of moments based on complete
case; CC-OLS, ordinary least squares regression based on complete case

6. Real data analysis
We applied our methods to an example from Acemoglu et al. (2001), which
consider the e¤ect of institutions on economic performance. Acemoglu et al.
(2001) surveyed the causal relationship between “institution”(such as, moresecure property rights and less distortionary policies) and “economic performance” proxied by GDP per capita. It is indicated that better institution
and GDP per capita has positive relationship, but the causal relationship had
not been proven. To solve the endogeneity problem, Acemoglu et al. (2001)
employed European settlers mortality as an instrument. They argue thatcountries with higher mortality rate are …xed to be extractive institutions.
We used the same dataset as Acemoglu et al. (2001) and complete caseanalysis results in N = 64. However, if instruments are available, 111 samples are ready for analysis given they completely observed other variables.
The outcome in this dataset is natural logarithm of GDP (logGDP ), andthe endogenous variable is Average Protection Against Ex-propriation Risk
(AP ER). Available exogenous variables are latitude and country dummy(Asia, Africa, and other countries). We use the natural logarithm of annualized European settlers’mortality per thousand mean strength (logSM ) as
an instrument, of which 36:9% are missing. Since Acemoglu et al. (2001)
constructed the IV sourced from and estimated by other research, i.e., Curtin(1989) and Curtin (1998), it tends to be missing. The F-statistic of reducedform equation is 2:25, indicating weak instruments; hence, using more samples is e¤ective in obtaining reliable results.
As population-level information, we use the information obtained from
Albouy (2012), which improved the data for logSM . Based on AIC …tted toAlbouy (2012)’s logSM , we use log-normal distribution with mean parameter
1:50 and standard deviation parameter 0:28 (AL-ln), and normal distribution
with mean 4:65 and standard deviation 1:24 (AL-n). Another candidate is thecomplete case of Acemoglu et al. (2001), which has log-normal distribution
with mean parameter 1:47 and standard deviation parameter 0:35 (AC-ln).
Figure 3 describes the density of AL-ln, AL-n, and AC-ln.

17

Fig.3 Probability density function of assumed population distribution of logSM
based on AL-ln, AL-n, and AC-ln.

Burn-in of 25; 000 iterations followed by 50; 000 iterations were used forthe posterior inference.
We also compared the results obtained from the pseudo-proposed methodwith do not consider the original population distribution (Ps-Proposed),MICE-GMM, FM-GMM, CC-GMM.
The results are described in Table 2. Three estimates based on the proposed method (AL-ln, AL-n, and AC-ln) imputed the missing values, andnone of the three have 0 for their credible intervals. However, they employed di¤erent distributional information and all the posteriors are di¤erent.
Pseudo-proposed methods, which do not have identi…cation, yield uninterpretable results. MICE-GMM and CC-GMM show signi…cant coe¢ cients,
but the means are smaller than the proposed methods (AL-ln, AL-n, andAC-ln). MF-GMM shows the smallest coe¢ cient with statistical signi…cance.
As discussed before, the estimates from missForest result in smaller standard
deviation, hence the obtained s.d. seems to abnormally small.
Given that the original population distributions are true, the causal e¤ectsof AP ER to GDP is much larger than the prior surveys.

18

7. Discussion
We have proposed semiparametric Bayes IV method which incorporates themissing instruments with NMAR using the original population information of
IV. In this paper, we assume that the missing IV is one-dimensional. Within
the …eld of economics, our method can deal with several datasets, since almost
all the data have only one IV. However, if we consider application to other…elds such as biometrics, the extension to multiple instruments might berequired. As described, for example, Mendelian randomizations use multiple
instruments, and IV used for them can be missing since they are genetic
information. We should consider another identi…cation condition for multiple
IVs with NMAR. The condition might vary by assumptions, for instance, if
joint distribution of original population is available, or if just the original
marginal distribution of each IV is available.
Acknowledgements
This work was supported by the JSPS Grant-in-Aid for Research Activity
Start-Up (18H05697) and the JSPS KAKENHI (18H03209; 17K18598).

References
Aaslund, O. and Gronquist, H. (2010). Family size and child outcomes: isthere really no trade-o¤. Labour Economics 17, 130–139.
Acemoglu, D., Johnson, S. and Robinson, J. A. (2001). The colonial origins of comparative development: An empirical investigation. American
Economic Review 91, 1369–1401.

19Albert, J. H. and Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American statistical Association 88,
669–679.
Albouy, D. Y. (2012). The colonial origins of comparative development:
an empirical investigation: comment. American Economic Review 102,
3059–3076.
Burgess, S., Seaman, S., Lawlor, D. A., Casas, J. P. and Thompson, S. G.
(2011). Missing data methods in mendelian randomization studies withmultiple instruments. American Journal of Epidemiology 174, 1069–1076.
Chaudhuri, S. and Guilkey, D. K. (2016). Gmm with multiple missing variables. Journal of Applied Econometrics 31, 678–706.
Chaudhuri, S., Handcock, M. S. and Rendall, M. S. (2008). Generalized
linear models incorporating population level information: an empiricallikelihood-based approach. Journal of the Royal Statistical Society: Series
B 70, 311–328.
Conley, T. G., Hansen, C. B., McCulloch, R. E. and Rossi, P. E. (2008). A
semi-parametric bayesian approach to the instrumental variable problem.
Journal of Econometrics 144, 276–305.
Curtin, P. D. (1989). Death by migration: Europe’s encounter with the tropical world in the nineteenth century. Cambridge University Press.
Curtin, P. D. (1998). Disease and empire: The health of European Troops in
the Conquest of Africa. Cambridge University Press.
Darolles, S., Fan, Y., Florens, J. P. and Renault, E. (2011). Nonparametric
instrumental regression. Econometrica 79, 1541–1565.
Ertefaie, A., Flory, J. H., S, H. and D., S. (2017). Instrumental variable
methods for continuous outcomes that accommodate nonignorable missing baseline values. American Journal of Epidemiology 185, 1233–1239.
Hall, P. and Horowitz, J. L. (2005). Nonparametric methods for inference inthe presence of instrumental variables. Annals of Statistics 33, 2904–2929.
Hellerstein, J. K. and Imbens, G. W. (1999). Imposing moment restrictions
from auxiliary data by weighting. Review of Economics and Statistics 81,
1–14.
Hirano, K. (2002). Semiparametric Bayesian inference in autoregressive panel
data models. Econometrica 70, 781–799.
Hirano, K., Imbens, G. W., Ridder, G. and Rubin, D. B. (2001). Combining
panel data sets with attrition and refreshment samples. Econometrica 69,
1645–1659.
Igari, R. and Hoshino, T. (2018). A bayesian data combination approachfor repeated durations under unobserved missing indicators: Application20

to interpurchase-timing in marketing. Computational Statistics and Data
Analysis 126, 150–166.
Imbens, G. W. and Lancaster, T. (1994). Combining micro and macro datain microeconometric models. Review of Economic Studies 61, 655–680.
Ishwaran, H. and James, L. F. (2011). Gibbs sampling methods for stickbreaking priors. Journal of American Statistical Association 96, 161–173.
Kato, K. (2013). Quasi-bayesian analysis of nonparametric instrumental variables models. Annals of Statistics 41, 2359–2390.
Kato, R. and Hoshino, T. (2019). Semiparametric bayesian multiple imputation for regression models with missing mixed continuous-discrete
covariates. Annals of the Institute of Statistical Mathematics, in press .
Kennedy, E. H. and Small, D. S. (2017). Paradoxes in instrumental variable studies with missing data and one-sided noncompliance. arXiv page1705.00506.
Kim, J. K. and Yu, C. L. (2011). A semiparametric estimation of meanfunctionals with nonignorable missing data. Journal of the American
Statistical Association 106, 157–165.
Klein, N., Kneib, T. and Lang, S. (2015). Bayesian generalized additivemodels for location, scale, and shape for zero-in‡ated and overdispersed
count data. Journal of the American Statistical Association 110, 405–419.
Kott, P. S. and Chang, T. (2010). Using calibration weighting to adjust
for nonignorable unit nonresponse. Journal of the American Statistical
Association 105, 1265–1275.
Liao, S. G., Lin, Y., Kang, D. D., Chandra, D., Bon, J., Kaminski, N.,
Sciurba, F. C. and Tseng, G. C. (2014). Missing value imputation in
high-dimensional phenomic data: imputable or not, and how? BMC
Bioinformatics 15, 346.
Liao, Y. and Jiang, W. (2011). Posterior consistency of nonparametric conditional moment restricted models. Annals of Statistics 39, 3003–3031.
Little, R. J. and Rubin, D. B. (2002). Bayes and multiple imputation. Statistical analysis with missing data pages 200–220.
Liu, J., Gelman, A., Hill, J., Su, Y. S. and Kropko, J. (2014). On thestationary distribution of iterative imputations. Biometrika 101, 155–
173.
Miao, W., Ding, P. and Geng, Z. (2016). Identi…ability of normal and normal
mixture models with nonignorable missing data. Journal of the American
Statistical Association 111, 1673–1683.
Mogstad, M. and Wiswall, M. (2012). Instrumental variables estimation withpartially missing instruments. Economics Letters 114, 186–189.

21

Nevo, A. (2003). Using weights to adjust for sample selection when auxiliary
information is available. Journal of Business and Economic Statistics 21,
43–52.
Newey, W. K. and Powell, J. L. (2003). Instrumental variable estimation of
nonparametric models. Econometrica 71, 1565–1578.
Palmer, T. M., Lawlor, D. A., Harbord, R. M., Sheehan, N. A., Tobias,
J. H., Timpson, N. J., Smith, G. D. and Sterne, J. A. (2012). Using multiple genetic variants as instrumental variables for modi…able risk factors.
Statistical methods in medical research 21, 223–242.
Qin, J. (2000). Combining parametric and empirical likelihoods. Biometrika
87, 484–490.
Shah, A. D., Bartlett, J. W., Carpenter, J., Nicholas, O. and Hemingway, H.
(2014). Comparison of random forest and parametric imputation models
for imputing missing data using mice: a caliber study. American Journalof Epidemiology 179, 764–774.
Shahbaba, B. and Neal, R. (2009). Nonlinear models using dirichlet process
mixtures. Journal of Machine Learning Research 10, 1829–1850.
Shen, W., Tokdar, S. T. and Ghosal, S. (2013). Adaptive Bayesian multivariate density estimation with dirichlet mixtures. Biometrika 100, 623–640.
Smith, G. D. (2006). Randomised by (your) god: robust inference froman observational study design. Epidemiology and Community Health 60,
382–388.
Stekhoven, D. J. and Buhlmann, P. (2012). Missforest–non-parametric missing value imutation for mixed-type data. Bioinformatics 28, 112–118.
van Buuren, S. (2007). Multiple imputation of discrete and continuous data
by fully conditional speci…cation. Statistical Methods in Medical Research
16, 219–242.
Waljee, A. K., Mukherjee, A., Singal, A. G., Zhang, Y., Warren, J., Balis, U.,
Marrero, J., Zhu, J. and Higgins, P. D. (2013). Comparison of imputation
methods for missing laboratory data in medicine. BMJ open 3, e002847.
Xue, Y., Liao, X., Carin, L. and Krishnapuram, B. (2007). Multi-task learning for classi…cation with dirichlet process priors. Journal of Machine
Learning Research 8, 35–63.

22

Appendix A
Detailed simulation design
A1. Setup for simulation study 1
In this study, we assume 1-dimensional endogenous regressor xi , 4-dimensional
exogenous regressor wi , whose …rst elements are set to be 1 (i.e. intercept),
and 1-dimensional instruments z i , whose …rst elements are set to be 1 andare independent of 2;i . The true relationship of the full model is as follows.
yi = xi + wi + 1;i; ixi = z i + 2;i1;iN (0; )

2;i

1
0:5
; and
0:5
1
the interested parameter= 1: As for z i , we …rst generate z 1;i from
log-normal distribution0with mean
parameter0
1 0 0 and variance
11parameter 0:5:
0:2
1 0:3 0:2
Then, w 1;i
M V N @z 1;i @ 0:2 A ; @ 0:3 1 0:3 AA ; where w 1;i0:2
0:2 0:3 1
denote the components of wi with the …rst component 1 removed, and set
z i to be z i = 1; z 1;i ; wi :After creating complete dataset, we set some components of IV to be
missing. We …rstly calculated
We seti = Ui= (1; 1; 0; 1; 1)t ;1( 1 yi += (1; 1; 1:1)t ;2 xi +

3 wi +

4 zi +=

5 yi zi +6 xi zi +7 w i zi )where Ui
unif orm(0; 1); and corresponding i of zi which has the top
30% of i are converted to be missing.
are set to be ( 1 ; 2 ; :::; 7 ) =(0:25; 0:25; 0; 0:25; 0; 0:5; 0) : Since the missing probability of IV depends
on IV, this is the NMAR case. We assume parametric probit model withcross terms.
The results are provided in the main article.
A2. Setup for simulation study 2
Simulation study 2 consider the …nite dimensional mixture of regression
model to the reduced-form equation. In this study, we assume 1-dimensional
endogenous regressor xi , 4-dimensional exogenous regressor wi , whose …rst
elements are set to be 1 (i.e. intercept), and 1-dimensional instruments z i ,whose …rst elements are set to be 1 and are independent of 2;i . The true
relationship of the model is as follows.
23

8< yi = xi + wi + 1;i2P; i
1i2l z i l + 2;i: xi =

1;i
2;i

N (0;

l)

l=1

where c is taken so that the mean of

i is zero.

For weak instrument case,= ( 1; 2; 2; 0; 0)t ; 2 = ( 1; 1; 0; 0; 0)t ;
1 0:5
2 1
= (1; 1; 1:1)t ; ( 1 ; 2 ) =;; and the inter0:5 1
1 2
ested parameter = 1: As for z i , we …rst generate z 1;i from log-normal
distribution
mean parameter0 with0
1 0 0 and variance
11parameter 0:5: Then, w 1;i0:2
1 0:3 0:2
M V N @z 1;i @ 0:2 A ; @ 0:3 1 0:3 AA ; where w 1;i denote the0:2
0:2 0:3 1
components of wi with the …rst component 1 removed, and set z i to bez i = 1; z 1;i ; wi :The data missing process and the hyperparameters are the same as simulation 1.

we set ( 1 ;2)= (0:5; 0:5);

1

A2.1 Results of Simulation study 2
Table A1 shows the results of simulation 2, including the empirical mean,
standard deviation, the coverage of nominal 95% con…dence intervals (CI)of the estimate, and the mean squared error (MSE) from the true value of. Figure 5-1 shows the box plot of the estimates for obtained by 1; 000
replications.
As can be seen from the table, this simulation study considers the casewhere MICE works relatively well. However, it shows biased results andthe coverage is poor (75.3%). The MSE is about 5.0 times larger than the
proposed. MF-GMM shows the largest MSE, and the coverage is very poor.
Ps-Proposed, which does not have the information of the original distributionof IV, contains missing mechanism that cannot be identi…ed and the obtained
estimates are slightly upwardly biased. It should be noted that standard
deviation obtained from Ps-Proposed is much lower than others, but thecoverage is poor (72.6%). As a result, the MSE of Ps-Proposed is 4.1 times
that of the proposed. The proposed obtained the smallest MSE with smallers.d. and better coverage. Comparing the results of the proposed with PsProposed and MICE-GMM, original population information combined withsemiparametric IV model speci…cation seems very useful.
The complete case analysis-based methods also shows poor results. CCGMM estimates have upwardly biased results and the MSE is about 7.1 times
24that of the proposed method. The coverage of 95% con…dence intervals of
CC-OLS is only 3:3%, since OLS-CC ignores the presence of endogeneity.

25