Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
103,359 characters · 8 sections · 68 citation commands
Seemingly Unrelated Regression with Measurement Error: Estimation via Markov chain Monte Carlo and Mean Field Variational Bayes Approximation
The seemingly unrelated regression (SUR) model consists of a system of linear multiple regression equations such that each equation has a different continuous dependent variable with a potentially different set of exogenous explanatory variables (covariates) and the errors are correlated across equations Zellner-1962. When the conditions of the SUR model apply, estimators obtained from SUR are more efficient relative to ordinary least squares estimators. The optimality feature and other theoretical properties of the SUR estimator within the frequentist framework are well studied in Srivastava-Dwivedi-1979, Srivastava-Giles-1987 and Fiebig-2001. The Bayesian approach to estimating SUR model was introduced in Zellner-1971, where the author analytically derived the conditional posterior densities of the parameters. Given the conditional posteriors, the model can then be estimated using a Markov chain Monte Carlo (MCMC) technique, known as Gibbs sampling Geman-Geman-1984, Casella-George-1992. Since the introduction in Zellner-1971, the literature on Bayesian analysis of SUR has grown considerably in various directions, including estimation via MCMC Percy-1992,Griffiths-Chotikapanich-1997,Griffiths-Valenzuela-2006 and direct Monte Carlo approach Zellner-Ando-2010,Ando-Zellner-2010, prediction in SUR model Percy-1992 and several model extensions that include restricted SUR Steel-1992, SUR with serially correlated errors and time varying parameters Chib-Greenberg-1995 and semiparametric inference in SUR model Koop-etal-2005.
The existing literature on SUR models including the quoted articles have worked based on the assumption that the covariates are measured correctly. Nonetheless, in practice there can emerge situations where one or more of the covariates are recorded with error, thus giving rise to SUR with measurement error (hereafter SURME). Modeling measurement error within a SUR structure or more generally in a multi-equation system has largely gone unnoticed in the literature (both frequentist and Bayesian), the only exception is Carroll-etal-Book-2006,Carroll-etal-2006 explained in the next paragraph. In contrast, there has been considerable work on single equation models with measurement error. Within a linear regression framework, it is well known that measurement error in the data leads to bias and inconsistency in ordinary least squares (OLS) estimator (see for instance Cheng-VanNess-1999, Fuller-1987, Wansbeek-Meijer-2000, Rao-etal-2008 and Hu-Wansbeek-2017). To achieve consistency of OLS estimator, side assumptions are required such as known measurement error variance or known reliability ratio.\footnote{ If $w$ and $z$ are two random variables such that $w=z+u$ and the error $u$ is independent of $z$, then the reliability ratio $R_{z}$ is defined as the true variance divided by the total variance, i.e., $ R_{z}=Var(z)/(Var(z)+Var(u))$. By definition $0\leq R_{z}\leq 1$.} However, consistent estimator of regression parameters without the side assumptions can be constructed when measurement errors have replicated observations Shalabh-2003. Measurement error in nonlinear models is discussed in Carroll-etal-2006 along with the Bayesian analysis of linear and non-linear measurement error models.
Within the multi-equation framework, Carroll-etal-2006 consider a combination of linear mixed measurement error model and SUR model to understand the properties of measurement error in food frequency questionnaire data for protein and energy. They adopt the frequentist estimation approach and use a nearby adaptive method based on weighted Akaike information criterion (AIC) to select the best fitting model, a form of model averaging which is popular in the Bayesian literature. Carroll-etal-2006 find that a fully parameterized model in which measurement errors in the two nutrients are modeled jointly, offers no gain in efficiency compared to fitting each model separately. However, when some parameters are set to zero resulting in a reduced model, considerable gains in efficiency is attained. We may adopt the frequentist approach to estimating SURME model with structural measurement error, but it is fraught with difficulty because the number of parameters become larger than the number of normal equations derived from the likelihood function. In such cases, side assumptions can be used to identify the model as done in linear regression, but even then deriving the maximum likelihood estimators for SURME model is a challenging task. Besides, ignoring measurement error in the data can lead to a poor model fit.
In this paper, we introduce two novel methods---a pure Bayesian algorithm and a mean field variational Bayes (MFVB) technique---to estimate the SURME model where each equation can potentially have a different covariate that is measured with error. Both the approach employs a classical structural form of measurement error and the link between the covariate measured with error and the other covariates (with no measurement error) is modeled through an exposure equation. Identification of parameters is achieved by placing a prior distribution on the measurement error variance. The pure Bayesian approach is analytically simpler and produces tractable conditional distributions which enables the use of Gibbs sampling. However, the MCMC draws of the parameters corresponding to the covariate measured with error tend to be highly correlated. To reduce autocorrelation in MCMC draws, one may consider thinning i.e., use every $l$-th draw in estimating the parameter. Thinning is debatable and while some authors such as Owen-2017 recommend thinning, others such as Link-Eaton-2012 advise against the use of thinning. So, we explore other methods and come up with a more elegant solution to the problem of high autocorrelation i.e., the MFVB approach to estimating SURME model.
We illustrate both the techniques in multiple simulation studies and compare the results to a standard SUR model, where we ignore or do not model the measurement error. In the first set of simulation studies, data are generated from a SURME model using different values of the variance of the true unobserved variable, while holding the reliability ratio fixed. In the second set of simulations, data are generated using different values of reliability ratio, holding the variance of the true unobserved variable at a fixed value. The results suggest that both the proposed methods perform well and correctly highlight the importance of modeling measurement error within the SUR structure when variables are measured with error. In addition, the SURME model is implemented in an application drawn from the health literature and estimated using the two proposed methods. Specifically, weight and high density lipoprotein (HDL) are jointly modeled as a function of several covariates and blood pressure, which is common to both equations and considered to have measurement error. Blood pressure is modeled as a function of the covariates in the exposure equation. Model selection exemplify the practical utility of the SURME model compared to the standard SUR model.
The remainder of the paper is organized as follows. Section 2 presents the SURME model, derives the joint posterior density and proposes a Gibbs sampling algorithm to estimate the model. Section 3 develops the MFVB approximation of the MCMC algorithm. Section 4 demonstrates the two algorithms in several Monte Carlo simulation exercises and Section 5 presents an application drawn from health literature. Section 6 concludes.
The seemingly unrelated regression with measurement error (SURME) model incorporates measurement error for covariates in the SUR model and can be expressed in terms of the following equations,
where the response $y_{mi}$ is a scalar, $x_{mi}^{\prime }$ is $\left( 1\times k_{m}\right) $ vector of covariates, $z_{mi}$ is a true unobserved scalar covariate that is prone to measurement error, and the subscripts $m$ and $i$ denote the equation number and individual/observation, respectively. Stacking the equations for each $i$, we can write model ((ref)) as follows,
where $y_{i}=\left( y_{1i},...,y_{Mi}\right) ^{\prime }$ and $\gamma =\left( \gamma _{1},...,\gamma _{M}\right) ^{\prime }$ are vectors of dimension $\left( M\times 1\right)$, and $\beta =\left( \beta _{1},...,\beta _{M}\right) ^{\prime }$ is of dimension $(K \times 1)$, where $K = k_{1}+\cdots + k_{M}$. The matrices,
are of dimension $(M \times K)$ and $(M \times M)$, respectively. In addition, the error $\varepsilon_{i}$ is assumed to be independently and identically distributed (i.i.d.) as a normal distribution i.e., $\varepsilon_{i} \sim N(0, \Sigma_{\varepsilon})$ for $i=1,\cdots,N$, where the covariance,
is a symmetric matrix that permits nonzero correlation across equations (or first subscript) for any given individual (or second subscript) and ties each independent regression into a system of equations, hence the phrase seemingly unrelated regression. Measurement error in reference to model ((ref)) arises because $Z_{i}$ is not observed, instead we observe $W_{i}$ which is a sum of the true unobserved quantity $Z_{i}$ and a measurement error term $u_{i}$. This definition implies a classical measurement error Fuller-1987. Additionally, we assume that the true unobserved quantity $Z_{i}$ follows a distribution, so that the measurement error model is of the structural form. This can be represented as follows,
where for algebraic simplification, we use the notations $\widetilde{W}_{i}=\left( w_{1i},...,w_{Mi}\right) ^{\prime }$, $\widetilde{Z}_{i}=\left( z_{1i},...,z_{Mi}\right) ^{\prime }$, $\widetilde{u}_{i}=\left( u_{1i},...,u_{Mi}\right) ^{\prime }$, then $W_{i}=diag(\widetilde{W}_{i})$, $ Z_{i}=diag(\widetilde{Z}_{i})$, $u_{i}=diag(\widetilde{u}_{i})$ are $\left( M\times M\right) $ diagonal matrices and $I_{M}$ is a ($M \times M$) identity matrix.
An interesting addition to equation (ref) is to relate the primary explanatory variable of interest (here $Z_{i}$) to other covariates ($X_{i}$), giving rise to the exposure model. The term “exposure model” comes from epidemiology, where the primary explanatory variable is affected by exposure to “toxicants” or “risk factors”. Therefore, the potential links between the latent variable $Z$ and the other covariates $X$ can be expressed as follows,
The three equations (ref), (ref) and (ref) together define our SURME model and the resulting likelihood is derived as follows,
where $\Delta \equiv (\beta, \gamma, \Sigma_{\varepsilon}, \omega, \sigma_{Z}^{2},\sigma_{u}^{2})$ and as mentioned earlier, $\widetilde{W}_{i}$ and $\widetilde{Z}_{i}$ are column vectors that contain the diagonal elements of the matrices $W_{i}$ and $Z_{i}$, respectively.
Before proceeding with estimation, we add a few words on identification issues that typically arise with measurement error models. In linear regression with measurement error, identification of parameters require additional assumptions. Such assumptions can be constant measurement error variance, known reliability ratio or some other conditions as presented in Cheng-VanNess-1999. The same identification conditions are also applicable to the proposed SURME model under the existing distributional assumptions. Nonetheless, we follow a purely Bayesian approach and employ prior distributions to identify the parameters of the model Zellner-1971.
The Bayesian estimation method combines the likelihood of the model with suitable prior distributions to obtain the joint posterior distribution. We utilize the following prior distributions:
where $W_{M}$ denotes a Wishart distribution of dimension $M$ and $IG$ denotes an inverse gamma distribution. Here we note that if one is not interested in the exposure equation, it can be dropped from the model. In such a case, $\tilde{Z}_{i} \sim N(\mu, \sigma_{Z}^{2} I_{M})$ and $\mu$ can be given a normal prior as $\mu \sim N(\mu_{0}, \sigma_{\mu}^{2} I_{M})$. Coming back to the SURME model, the joint posterior distribution can be obtained by combining the likelihood (ref) with the prior distributions (ref) as follows,
Typical with the Bayesian approach, the joint posterior density ((ref)) is not tractable and the parameters are sampled using MCMC techniques. To this purpose, conditional posterior densities of the parameters are derived (see Appendix A in the supplementary material) and Gibbs sampling is employed to estimate the model as exhibited in Algorithm (ref). Note that some of the conditional posteriors are conditioned on a subset of parameters, but these are full conditionals that just do not depend on the full set of parameters. Conditional posteriors that depend on a subset of parameters have also been referred to as reduced conditional posteriors and Gibbs sampling as partially collapsed Gibbs sampling Liu-1994,vanDyk-Park-2008.
The sampling algorithm, presented in Algorithm (ref), shows that $\beta$ and $\gamma$ are sampled from an updated multivariate normal distribution. Standard result is obtained for the precision matrix $\Sigma_{\varepsilon}^{-1}$, which is sampled from an updated Wishart distribution. All the three parameters ($\beta,\gamma,\Sigma_{\varepsilon}$) follow their respective distributions marginally of $\omega$, $\sigma_{u}^{2}$ and $\sigma _{Z}^{2}$. The true unobserved quantity $Z$ is drawn from an updated multivariate normal distribution conditional on all the remaining model parameters. Similarly, $\omega $ is sampled from an updated multivariate normal distribution conditional on $\left( Z,\sigma _{Z}^{2}\right) $. The two variance parameters are drawn from updated inverse gamma distributions with $\sigma _{Z}^{2}$ conditioned on $\left(Z,\omega \right) $ and $\sigma _{u}^{2}$ conditioned on $\left(W, Z \right)$. Note that if we drop the exposure equation from the SURME model, Algorithm (ref) only requires a slight modification. In this context, $\mu$ replaces $X_{i}\omega$ and is sampled from an updated normal distribution as $\mu|Z,\sigma_{Z}^{2} \sim N_{M}(\bar{l}, \bar{\Lambda})$, where $\bar{l} = \bar{\Lambda} \Big( \sum_{i=1}^{N} \widetilde{Z}_{i} / \sigma_{Z}^{2} + \frac{\mu_{0}}{\sigma_{\mu}^{2}} \Big)$ and $\bar{\Lambda}^{-1} = \Big( \frac{N}{\sigma_{Z}^{2}} + \frac{1}{\sigma_{\mu}^{2}}\Big) I_{M}$ are the posterior mean and posterior precision, respectively.
We note that the model presented in this paper utilizes the structural measurement error model which assumes that $Z$ follows a distribution. Hence, the distribution of $Z$ was introduced as a part of the model. However, in the measurement error literature, there is another form of measurement error known as functional form. The functional measurement error model assumes that the true unobserved quantity $Z$ is fixed. In our modeling and estimation framework, we can easily incorporate the functional form of measurement error by modeling the distribution of $Z$ as a part of the subjective prior information Zellner-1971. This implies that the joint posterior distribution (ref) will be unchanged and derivations of the conditional posterior distributions will proceed in exactly the same way as described in Appendix A of the supplementary file. To reiterate, the fundamental difference in analyzing SUR model with structural and functional forms of measurement error lies in the interpretation given to the distribution of $Z$, the true unobserved quantity.
In the MCMC estimation of SURME model, one consideration that arise is that $Z$ and $\gamma $ are both unknown, and drawing them conditional on each other lead to high autocorrelation in MCMC draws. This occurrence is a general problem and happens when two or more unknown variables/parameters that appear in product form are drawn conditional on each other. To reduce the autocorrelation in MCMC draws (and consequently reduce the inefficiency factors) some authors\footnote{See for instance Jeliazkov-2013 for the case of latent variables in a non parametric VAR specification.} propose to improve mixing by sampling $\gamma $ from the marginal distribution and then sampling $Z|\gamma $ or vice versa. However, deriving the marginal posterior distribution of $\gamma$ (or of $Z$) is not straightforward and the marginalization trick do not improve the results in our modeling context. \footnote{Many thanks to Ivan Jeliazkov and the participants of the UCI seminar for the suggestion to sample $\gamma$ marginally of $Z$ and then sampling $Z|\gamma$. See appendix E in the supplementary material. However, the several tests we conducted did not improve our initial results with standard Gibbs sampling.} As a solution to reduce autocorrelation, many researchers have employed thinning to improve the mixing of the draws. The thinning of MCMC draws has been criticized by some authors including MacEachern-Berliner-1994, Link-Eaton-2012), but others such as Geyer-1991 acknowledges that thinning can increase statistical efficiency. In a recent paper, Owen-2017 shows that the usual advice against thinning can be misleading. We employ thinning to improve the mixing properties of the MCMC draws of $\gamma$ in our simulation studies and application. However, given the controversy around thinning, we explore other methods and come up with the MFVB approximation to estimate SURME model.
Variational Bayes is an alternative to MCMC methods that provides a locally-optimal, exact analytical solution to an approximation of the posterior distribution. The parameters of the approximate distribution are selected to minimize the Kullback-Leibler divergence (a distance measure) between the approximation and the posterior. The MFVB approximation is a deterministic optimization approach and so is particularly useful for big data sets and/or models with large sparse covariance matrices. Besides, it is similar to Gibbs sampling for conjugate models. Some recent articles on MFVB approach include Bishop-2006, Ormerod-Wand-2010, Pham-etal-2013, Lee-Wand-2016 and Blei-etal-2017.
Suppose, $y$ denotes an observed data vector and $\theta$ is a parameter vector defined over the parameter space $\Theta$. Following the Bayes theorem, the posterior distribution can be written as:
where $p\left( y\right) =\int_{\Theta }p\left( \theta ,y\right) d\theta$ is the marginal likelihood. Let $q$ be an arbitrary density function over $\Theta $. Then, the logarithm of the marginal likelihood is,
where $\log \underline{p}\left(y,q\right) = E_{q(\theta)} \left[ \log \left( \frac{ p\left( \theta, y\right)}{ q\left( \theta \right)} \right) \right]$ denotes the lower bound on the marginal log-likelihood and $KL(q,p)= E_{q\left( \theta \right) }\left[ \log q\left( \theta \right) \right] -E_{q\left( \theta \right) }\left[ \log p\left( \theta |y\right) \right]$ is the Kullback-Leibler divergence $q\left(\theta \right) $ and $p\left( \theta|y\right)$. Since $\log p(y)$ is a constant, the minimization of $KL(q,p)$ is equivalent to maximizing the scalar quantity $\log \underline{p}\left(y,q\right)$, typically known as evidence lower bound (ELBO) or variational lower bound. In practice, the maximization of the ELBO is often preferred to minimization of the KL divergence since it does not require knowledge of the posterior.
The MFVB approximates the posterior distribution $p\left( \theta|y\right) $ by the product of the $q$-densities,
Each optimal $q$-density minimizes the Kullback-Leibler divergence and is given by,
where $E_{q\left( -\theta _{j}\right)}$ denotes the expectation with respect to $\prod_{k\neq j}q_{k}\left( \theta _{k}\right)$, $ \Omega\equiv \left\{ y,\theta _{1},...,\theta _{j-1},\theta _{j+1},...,\theta _{P}\right\} $ is the set containing all random vectors in the model except $\theta _{j}$, and $p\left(\theta _{j}|\Omega\right) $ are the full conditional distributions of the parameters.
For the SURME model, we now consider a MFVB approximation based on the following factorization:
These optimal $q$-densities can be derived, as presented in Appendix B of the supplementary file, to have the following form,
where $f$ denotes the density function of the distribution given in the subscript. The parameters of the optimal densities are updated according to Algorithm (ref). When exposure is dropped, $q(\omega)$ is replaced with $q(\mu)$ and the optimal density is an updated normal distribution. Convergence of Algorithm (ref) is assessed using the evidence lower bound $\ell$ on the marginal log-likelihood (see Appendix C in the supplementary material) that is guaranteed to reach a local optima based on the convexity property. This algorithm belongs to the family of coordinate ascent variational inference (CAVI) and iteratively optimizes each factor of the mean field variational density, while holding the remaining fixed Bishop-2006, Blei-etal-2017.
The MFVB technique provides computational advantages compared to MCMC because it is deterministic and does not require a large number of iterations Pham-etal-2013,Lee-Wand-2016. Besides, existing works including Bishop-2006, Ormerod-Wand-2010, Faes-etal-2011, Pham-etal-2013, and Lee-Wand-2016 suggest that the accuracy scores of the MFVB approximation, relative to MCMC, generally exceed $95-97\%$ and rarely goes below $90\%$. Given these advantages, the MFVB approach can be gainfully utilized for large data models. However, some authors have reported that covariance matrices from variational approximation may be typically “too small” relative to the sampling distribution of the maximum likelihood estimator. In this regard, Blei-etal-2017 opine that underestimation of the variance should be judged in relation to the task at hand. However, evidence from empirical research indicate that variational inference typically do not suffer in accuracy.
This section examines the performance of the two proposed methods in multiple simulation studies. The first set of simulations (Case I) employ different values of $\sigma _{Z}^{2}$ to generate the simulated data. The second set of simulations (Case II) use different values of reliability ratio defined as $R_{Z}=\sigma _{Z}^{2}/(\sigma _{Z}^{2}+\sigma _{u}^{2})$. In both sets of simulations, we use a two equation structure represented as follows,
where the first, second and third subscripts in $x_{mij}$ denote the equation number ($m=1,2$), observation ($i=1,\cdots,N$) and variable number ($j=2,3$), respectively. The first variable is common to both the equations (i.e., $x_{1i2}=x_{2i2}$ for all $i=1,...,N$) and the remaining covariates are exclusive to the respective equations. Moreover, we assume the error prone covariate $Z_{i}=diag(z_{1i},z_{2i})$ for all $i$ is unobserved, but is defined by an exposure model as follows,
The unobserved $Z_{i}$ is related to the observed $W_{i}=diag(w_{1i},w_{2i})$ by the equations below,
Note that the estimation of the SURME model solely relies on $W$ and the role of $Z$ is limited to generating values for $(W,y)$.
To proceed with data generation, we assign specific values to the parameters $\beta $, $\gamma $, $\omega$, $\Sigma _{\varepsilon }$, $\sigma _{Z}^{2}$, $\sigma _{u}^{2}$ and generate $N=300$ observations in each simulation study for all the variables in the model. Let $\beta _{11}=3$, $\beta _{12}=5 $, $\beta _{13}=4$, $\beta _{21}=4$, $\beta _{22}=3.8$, $\beta _{23}=3$, $ \gamma _{1}=4$, $\gamma _{2}=4$, $\omega_{11}=1.5$, $\omega_{12}=0.75$, $\omega_{13}=0.3$, $\omega_{21}=1.5$, $\omega_{22}=1.05$, and $\omega_{23}=0.45$. For all values of $i$, the error vector $ \varepsilon _{i}=\left( \varepsilon _{1i},\varepsilon _{2i}\right) ^{\prime } $ is generated from a bivariate normal distribution $N(0_{M},\Sigma _{\varepsilon })$, where $\Sigma _{\varepsilon }=[1$ $0.5$; $0.5$ $1]$. Values for the common covariate ($x_{1i2}=x_{2i2}$) are generated from $U(0,2)$ and values for the exclusive covariates $ x_{1i3}$ and $x_{2i3}$ are generated from $U(0,4)$, where $U$ denotes an uniform distribution. Values for $\widetilde{Z}_{i}=\left( z_{1i},z_{2i}\right) ^{\prime }$ are generated as $\widetilde{Z}_{i}\sim N$($X_i\omega ,\sigma_{Z}^{2}I_{M}$), and the $\widetilde{W}_{i}$'s are generated as $\widetilde{W} _{i}=\widetilde{Z}_{i}+\widetilde{u}_{i}$ where $\widetilde{u}_{i}\sim N$($0_{M},\sigma _{u}^{2}I_{M} $). The above setting remains the same in the following subsections, with change occurring only in values of $\sigma _{u}^{2}$ (through $R_{Z}$) or $ \sigma _{Z}^{2}$.
In Case I, we investigate the performance of the proposed algorithms in two simulation studies where the reliability ratio $R_{Z}$ is fixed ($R_{Z}=0.8$) and $\sigma _{Z}^{2}$ is gradually decreased. Specifically, two values are considered $\sigma_{Z}^{2}=\{1, 0.0625 \}$. The definition of $R_{Z}$ is used to generate the corresponding values for $\sigma_{u}^{2}=\sigma_{Z}^{2}(1-R_{Z})/R_{Z}$, which leads to a noise-to-true variance ratio $(1-R_{Z})/R_{Z}$ of 25%. In Case II, we again examine the performance of the proposed algorithms in two simulation studies by keeping $\sigma_{Z}^{2}$ fixed ($\sigma_{Z}^{2}=0.0625$) and using two values of reliability ratio $R_{Z}= \{0.8, \; 0.5714 \} $. The chosen values are similar to those used in Pham-etal-2013 and leads to noise-to-true variance ratios of 25% and 75%, respectively. We could define a noise-to-true variance of $100\%$, $150\%$ or more, but in those cases we will be dealing more with outliers than with measurement errors.
Bayesian procedures require prior distribution on the parameters of the model. For the SURME model, we stipulate the following priors: $\beta \sim N_{K}\left( \beta _{0},B_{0}\right) $ with $ \beta _{0}=\iota _{K}$, $B_{0}=I_{K}$, and $\iota _{M}$ is a $ \left( M\times 1\right) $ vector of ones; $\gamma \sim N_{M}\left( \gamma_{0},G_{0}\right) $ with $\gamma_{0}=\iota_{M}$, $G_{0}=I_{M}$; $\omega \sim N_{K}\left( \omega_{0}, O_{0}\right) $ with $\omega_{0}=\iota_{K}$, $O_{0}=I_{K}$; $\Sigma_{\varepsilon}^{-1} \sim W_{M}\left( \nu_{0}, S_{0}\right)$ with $\nu_{0}=50$ and $S_{0}= \nu_{0}[1$ $0.5$; $0.5$ $1]$; $\sigma_{Z}^{2}\sim IG\left( \delta_{1},\delta_{2}\right) $ and $\sigma_{u}^{2}\sim IG\left( \delta_{3},\delta _{4}\right)$ with $\delta_{1}=\delta_{2}=\delta_{3}=\delta_{4}=1/100$. All these priors are proper yet specify vague information about the parameters mainly for the measurement error $u_{i}$ and the error prone covariate $Z_{i}$. In addition, the same priors are used in all the simulations to highlight the effect of changing $R_{Z}$ or $\sigma _{Z}^{2}$ in estimation of the parameters and consequently on the performance of the algorithms.
The MCMC results are obtained from $50,000$ draws, after a burn-in of $1,000$ draws. We replicate these simulations $100$ times and report the means over these $100$ replications of the posterior means\footnote{ To save time, we only run $51,000$ draws per replication. Higher number of MCMC draws, such as $100,000$ or $200,000$, exponentially increase the computing time without any increase in precision. As an example for $\sigma_{Z}^{2}=0.0625$ and $R_{Z}=0.909$, the MFVB takes only $11.52$ seconds per replication. If $51,000$ (resp. $101,000$ and $201,000$) draws are used, the computing time per replication for the Gibbs sampling of the BSURME model is about $65.64$ (resp. $142.07$ and $352.29$) seconds using a MacBook Pro, 2.8 GHz core i7 with 16Go 1600 MHz DDR3 RAM.}. We also compare the results with the usual frequentist SUR estimation and the standard Bayesian estimation of SUR model. The Gibbs sampling algorithm for the latter is presented in Appendix D of the supplementary material.
Amongst the first set of simulation studies labeled Case I, Table (ref) presents the results from the frequentist estimation of SUR\footnote{Without any prior information on the measurement error, the SUR model for $M$ equations is the following: $y_{i}=X_{i}\beta +W_{i}\gamma +\varepsilon _{i}$ , $\varepsilon _{i}\sim N\left( 0,\Sigma _{\varepsilon }\right) $ , $i=1,..,N $ where $W_{i}$ is the covariate with measurement error.} model for the case $R_{Z}=0.8$ and $\sigma _{Z}^{2}=1$. Results show that estimates are strongly biased mainly for the intercepts $\beta _{11}$ and $\beta _{21}$, and for $\gamma _{1}$ and $\gamma _{2}.$ The relative biases ($\hat{\beta} / \beta - 1$) of the intercepts (resp. the $\gamma$'s) are $38.1\%$ and $28.8\%$, (resp. $-19.8\%$ and $-19.6\%$). The $\gamma$'s are strongly under-estimated. On the contrary, slope coefficients $\beta _{12}$, $\beta _{13}$ and $\beta _{23}$ are less contaminated by the measurement error and have a lower dispersion of the estimated coefficients than the intercepts. The relative bias of $\beta _{22}$ ($22.3\%$) is close (in absolute value) to that of $\gamma$'s. Elements of the variance-covariance matrix $\Sigma _{\varepsilon }$ are strongly over-estimated with a relative error of $317\%$ and $314\%$ for the variances and $4.6\%$ for the covariance $\sigma_{12}$. It leads to a strong under-evaluation of the coefficient of correlation $\rho_{\varepsilon_1 \varepsilon_2} =0.126$ far from the true value $\left(0.5\right) $. The Bayesian estimation of SUR model, presented in Table G1 of the supplementary material, give similar results. The posterior means of the coefficients (resp. posterior standard errors) are similar to the frequentist coefficient estimates (resp. standard errors) of the SUR model. The $95\% $ highest posterior density intervals (HPDI) are also close to the 95% confidence interval of the frequentist estimation. The estimated correlation coefficient $\rho_{\varepsilon_1 \varepsilon_2} =0.125$ is similar to the frequentist estimate. We also report Geweke's convergence diagnostic (CD), which tests for the equality of means of the first and last part of a Markov chain on the basis of samples drawn from the stationary distribution of the chain. In more than $98\%$ of cases, Geweke's CD (under the null hypothesis, $CD \sim N(0,1)$) accepts the null hypothesis at $5\%$ level, which suggests that a sufficiently large number of draws has been taken. Moreover, inefficiency factors (reported in Table G1) are also close to 1, which confirms that the chain is mixing well.
In the upper panel of Table G2 (of the supplementary material), we see how the Bayesian estimation of SURME model improves the results. The intercepts $\beta_{11}$ and $\beta_{21}$ are now less biased as compared to that of SUR model in Table G1. Their relative errors are $-3.5\%$ and $-7.1\%$, respectively. This is also true for the other slope coefficients $\beta$. Moreover, the model neatly corrects the measurement errors and results in better estimates of $\gamma_{1}$ and $\gamma_{2}$, their relative errors being $2.1\%$ and $2.6\%$, respectively. We also note that the variance-covariance matrix is precisely estimated leading to a correlation coefficient $\rho_{\varepsilon_1 \varepsilon_2} =0.494$. The parameters $\sigma_{Z}^{2}$ and $\sigma_{u}^{2}$ are well estimated with small posterior standard deviations and small relative errors ($-2.6\%$ and $0.6\%$, respectively). But inefficiency factors are large indicating strong autocorrelation in MCMC draws, particularly for $\gamma_1$ and $\gamma_2$ whose inefficiency factors are $8.62$ and $10.52$, respectively. In more than $90\%$ of cases, the Geweke's CD confirms that a sufficiently large number of draws has been taken. The improvement obtained with a SURME (as compared to the SUR) is interesting and emphasizes the need to model measurement error. The lower panel of Table G2 (in the supplementary material) presents the results of the exposure equation from the SURME model. They show that the biases are negligible, the posterior standard deviations are small and so are the inefficiency factors.
Overall, the SURME model is well estimated, but the high autocorrelation in MCMC draws of $\gamma$ needs additional consideration. According to Owen-2017, the problem of high autocorrelation can be dealt with thinning, which itself can be optimized according to the cost of computing the quantities of interest (after advancing the Markov chain) and the speed at which autocorrelations decay. As shown in Table G3 (see the supplementary material), autocorrelations between the successive draws of $\gamma_1$ and $\gamma_2$, denoted $\rho_{\tau} (\gamma_1)$ and $\rho_{\tau} (\gamma_2)$, are close to one and the rate of decay is very slow. For example, $\rho_{1} (\gamma_1)= 0.98$, $\rho_{10} (\gamma_1)= 0.82$ and $\rho_{1} (\gamma_2)= 0.98$, $\rho_{10} (\gamma_2)= 0.87$. The autocorrelations of some latent variables $Z_i$ are slightly higher ($0.995$) than those of the $\gamma$'s, but not reported for the sake of brevity. The cost of computing of $\tilde{Z}_i$ is on average $2.71$ and an autocorrelation of $0.995$ leads to an optimal thinning of factor $k=86$ (see appendix F and Table F1 in the supplementary material). Henceforth, we use a thinning of factor $k=100$ for all simulations.
We re-estimate the Bayesian SUR and SURME models with a thinning of 100, but only report the results for SURME. The results, presented in the upper panel of Table (ref), show that the posterior means and standard deviations are close (or identical) to those of Table G2 (in the supplementary material). Values of the inefficiency factors are small and are all between $(1.004, 1.23)$. Specifically, the reduction in inefficiency factor is tremendous for $\gamma$'s (e.g., $1.11$ versus $8.62$ for $\gamma_1$ and $1.23$ versus $10.52$ for $\gamma_2$). The lower panel of Table (ref) presents the results for the exposure equation in the SURME model. Once again, the results show that the biases are negligible, the posterior standard deviations are small and the inefficiency factors and Geweke's CD suggest good mixing of the MCMC draws. Specifically, the autocorrelations of $\gamma_1$ and $\gamma_2$ are now small ($\rho_{1} (\gamma_1)= 0.15$, $\rho_{1} (\gamma_2)= 0.27$) and quickly converge towards zero ($\rho_{10} (\gamma_1)= -0.007$, $\rho_{10} (\gamma_2)= -0.003$) confirming a good mixing of Markov chains (see Table G4 in the supplementary material).
We next discuss the results from MFVB approach, which on average takes about $145$ cycles to get the maximum of the evidence lower bound $l$ and the algorithm is terminated when the relative increase in the evidence lower bound $l$ is less than $10^{-7}$. The results from the MFVB estimation of SURME model are presented in Table (ref), which shows that the MFVB approach gives better results compared to Gibbs sampling. All parameters have similar or lower biases, mainly for the intercepts $\beta_{11}$, $\beta_{21}$ and for $\gamma$. However, the relative biases of the $\gamma$'s are now reduced ($0.7\%$ and $-0.3\%$) as compared to Bayesian estimation of SURME model. The estimates for $\sigma_{Z}^{2}$ and $\sigma_{u}^{2}$ show that the model accurately estimates the variances and their relative biases are small ($0.4\%$ and $-3.5\%$, respectively). The standard deviation of all the parameters are smaller compared to those from MCMC estimation leading to slightly narrower $95\%$ credible intervals (as compared to the $95\%$ HPDI)\footnote{When calibrating this Monte Carlo study, we found a significant underestimation of the variances of the coefficients $\gamma _{1}$ and $\gamma_{2}$, echoing the previous discussion around the work of Blei et al. (2017) (Section 3). After several trials (and to avoid embarking on more complex approaches such as linear response variational Bayes Giordano-etal-2018 or $\alpha $-variational inference Yang-etal-2018, we decided to use the following simple trick to correct this undervaluation: $\sigma_{\gamma_{j}}$ is replaced by $\sigma_{\gamma_{j}} \times \sqrt{MK/ E_{q\left( \sigma _{Z}^{2} \right) }}$, for $j=1,..,M$ (see Section B2 of the supplementary material).}. Additionally, estimates of $\sigma _{mm^{\prime }}$ are closer to their theoretical values and the estimated correlation coefficient $\rho_{\varepsilon_1 \varepsilon_2} =0.488$ is close to 0.5, the actual value. The lower panel of Table (ref) presents the results from the exposure equation which emphasizes the accurate estimation of the $\omega$ parameters. The MFVB approximation of both the classical structural form and the exposure model shows that there are definite advantages in adopting the MFVB approach to estimate measurement error models as compared to the pure Bayesian method.
We next decrease the variance $\sigma_{Z}^{2}$ from $\sigma_{Z}^{2}=1$ to $\sigma_{Z}^{2}=0.0625$ leading to $\sigma_{u}^{2}= 0.0156$. The results are presented in Tables G6 to G9 of the supplementary material. Results from the frequentist and Bayesian estimation of SUR model always reveal strong over-estimation of the intercepts, $\beta_{22}$ and strong under-estimation of the slopes $\gamma $ of the error prone covariate $Z_{i}$. However, over-estimation of the variances $\sigma_{11}$ and $\sigma_{22}$ ($ \simeq 19\%$ for both) are largely reduced as compared to those of $\sigma_{Z}^{2}=1$, but leads to a slightly under-estimated correlation coefficient $\rho =0.42$. When we incorporate the measurement error in the model, i.e., SURME model with a very small variance of $\sigma_{Z}^{2}$, the Bayesian estimates show a less accurate estimation of the intercepts (increasing the negative relative biases $-29.4\%$ and $-25.6\%$), of the $\gamma$'s ($15\%$ and $16.8\%$) and of all the $\beta$'s. Moreover, inefficiency factors rise to about $2$ indicating a relative loss of efficiency due to slightly correlated samples. To neutralize this effect, we can increase the thinning appropriately\footnote{We relaunched the simulations for this case with a thinning of $120$ and we find inefficiency factors close to $1$. The results are available upon request for the sake of brevity.}. Results for the exposure equation in the SURME model do not seem to be affected by the strong decrease of the variance $\sigma_{Z}^{2}$. The use of the MFVB approximation significantly attenuates the biases observed with the Bayesian estimation of SURME. The relative biases for the intercepts are now $-17\%$ and $-11\%$, and those for the $\gamma$'s are $8.6\%$ and $7\%$, respectively. The relative errors for the variances $\sigma _{mm},\, (m =1,2)$ reduce to $-2.5\%$ approximately and we get an estimated correlation coefficient $\rho =0.52$. The MFVB approximation accurately estimates parameters of the exposure equation. Once again, the MFVB method reveals its advantages in estimating a SUR model with measurement error although this advantage tends to be attenuated when a very small variance $\sigma_{Z}^{2}$ occurs.
In summary, for a fixed measurement error of $25\%$, increasing the variance $\sigma _{Z}^{2}$ of the error prone covariate $Z_{i}$ strongly biases the estimated variances $\sigma _{mm} \, (m=1,2)$ as well as the whole set of coefficients (intercepts and slopes) in the SUR model irrespective of the method of estimation. But, taking into account the measurement errors through SURME model neutralizes the negative effects of the increasing uncertainty on the error prone covariate $Z_{i}$ and thus strongly reduces, or even eliminates the biases to obtain satisfactory estimates. This conclusion is further reinforced with the use of the MFVB approximation.
We now investigate the performance of the proposed algorithms where $\sigma_{Z}^{2}=0.0625$ and the reliability ratio $R_{Z}$ is gradually decreased. Specifically, we consider $R_{Z}=\left\{\, 0.8, \, 0.5714\right\}$, which leads to $\sigma_{u}^{2}=\sigma _{Z}^{2}(1-R_{Z})/R_{Z} = \left\{0.0156, \, 0.0469\right\} $ and noise-to-true variance ratio $(1-R_{Z})/R_{Z}$ of $\left\{25\%, \, 75\%\right\}$. In the previous subsection, we have already studied the case where $\sigma_{Z}^{2}=0.0625$ and $R_{Z}=0.8$, therefore the focus is only on the case where the reliability is reduced to 0.5714.
The results presented in Tables G14-G17 of the supplementary material are poorer than those of $R_{Z}=0.8$ for both the frequentist and Bayesian estimates of SUR model, with stronger over-estimation of the intercepts ($84.3\%$ and $25.6\%$) and stronger under-estimation of the slopes $\gamma$ ($-42.5\%$). The relative biases of the intercepts are larger than in the previous cases and the same is true for the $\gamma$'s and even more obvious for the $\sigma_{mm} \, (m=1,2)$ (approximately $42\%$). The estimated correlation coefficient $\rho =0.38$ is far from the true value $0.5$. The Bayesian estimates of SURME model show a significant improvement, reducing the biases for the intercepts ($33\%$ and $-10.8\%$) and the $\gamma$'s ($16.2\%$ and $18.4\%$), but with slightly larger posterior standard deviations. The variance-covariance matrix $\Sigma_{\varepsilon}$ is well estimated, with small relative biases of $-7\%$ and $-9\%$ for $\sigma_{mm}\, (m=1,2)$. The estimated correlation coefficient turns out to be $\rho =0.53$. Both $\sigma_{Z}^{2}$ and $\sigma_{u}^{2}$ are also close to the true values (their relative biases are $-15.2\%$ and $15.1\%$, respectively). Once again, the improvement with the MFVB approximation is more noticeable as we get better results for the parameters with slightly smaller standard deviations. The relative biases for the intercepts are $-16.6\%$ and $-5.7\%$, and those of the $\gamma$'s are $8.5\%$ and $6.5\%$. For the slope coefficients, the relative biases range between $-11\%$ and $-2.5\%$. Both $\sigma_{Z}^{2}$ and $\sigma_{u}^{2}$ are also better estimated (their relative biases are respectively $-5.6\%$ and $4.4\%$). This is also true for the $\sigma_{mm}\, (m=1,2)$ (with small relative biases of $-1.9\%$ and $-1.6\%$) leading to an estimated correlation coefficient $\rho =0.51$.
To summarize, a change in the reliability ratio $R_{Z}$ --- for example, increasing the measurement error from $25\%$ to $75\%$ --- strongly biases the whole set of coefficients (intercepts and slopes), including the estimated variances $\sigma_{mm} \, (m=1,2)$ in the SUR model. This is true both for the frequentist and Bayesian approach. On the other hand, accounting for measurement error through SURME model largely eliminates the negative effects of this alteration and strongly reduces the biases in SUR estimation. Moreover, the use of MFVB approximation improves the results beyond those obtained with the Bayesian estimation of SURME model.\footnote{To get results with the Bayesian estimation of the SURME equivalent to those obtained with the MFVB approximation, it should be necessary to greatly increase the number of MCMC draws resulting in a very important cost in terms of computing time. Going from $51,000$ draws to $201,000$ draws leads to a relative increase in computing time per replication from $16$ to $86$ times that of MFVB (see note \footref{footnote_label_1}). There is therefore an obvious trade-off against the Bayesian estimation of the SURME and in favor of the MFVB approximation.}
Finally, we present in Table (ref) the relative errors of parameters for all the cases\footnote{For $\sigma_{Z}^{2}=1$ and $R_{Z}=0.5714$, results are given in Tables G13-G18. Last, Table G25 gives a summary of DICs and $p_D$s.} of $\sigma_{Z}^{2}=\left\{1, \, 0.0625\right\} $ and $R_{Z}=\left\{0.8, \, 0.5714\right\}$. At a glance, this table allows us to compare and contrast all of the previously discussed results and another case provided in the supplementary material. To reiterate, the results show that for a fixed measurement error, increasing the variance of the error prone covariate $Z_{i}$ strongly biases the whole set of coefficients (intercepts and slopes) as well as the estimated variances $\sigma_{mm} \, (m=1,2)$ in the SUR model, irrespective of the method of estimation. This is also the case when, for a fixed variance $\sigma_{Z}^{2}$, the reliability ratio $R_{Z}$ is reduced. Fortunately, taking into account the measurement error through SURME model attenuates or even neutralises the undesirable effects of the increasing uncertainty on the error prone covariate $Z_{i}$ or reducing the reliability ratio. This conclusion is further reinforced with the use of the MFVB approximation.
The current study utilizes data from the National Health and Nutrition Examination Survey (NHANES) for 2007-2008, a widely used survey designed to assess health and nutritional status of civilians, non-institutionalized adults and children in the United States. NHANES collects data by interviewing individuals at home, who then report to mobile examination centers (MECs) to complete the health examination component of the survey. The MEC's provide a standardized environment for the collection of high quality data, thus favoring dependable statistical estimation and interpretations. The survey is unique in the sense that it combines interviews and physical examinations of the respondents.
The dependent variables in the model are log of weight and high density lipoprotein ($HDL$), which is also known as `good cholesterol'. The covariates that are common to both equations and assumed to be measured without error are as follows: age, gender, smoking status, hours of sedentary activities, sleep disorder and low density lipoprotein plus 20 percent of Triglyceride ($LDL20T$). The variable `height' is only expected to effect weight, not $HDL$, and therefore only included in the log weight equation. Observed $SBP$ is assumed to have measurement error and transformed as $\ln (SBP-50)$ to avoid scaling problems, as done in Carroll-etal-2006. The third reading on $SBP$ is used as data and the first two readings are utilized to form priors on relevant parameters. Focussing on adults and removing missing observations on all variables of interest leaves us with a total of $N=1,001$ observations. Table (ref) presents the definition and descriptive statistics of all the variables used in the study.
To estimate the different SUR models with and without measurement error, we utilize the following relatively vague priors on the parameters: $\beta \sim N_{15}\left( \beta _{0},B_{0}\right) $ with $\beta _{0}=0_{15}$, $B_{0}=10 I_{15}$, $\gamma \sim N_{2}\left( \gamma _{0},G_{0}\right) $ with $\gamma _{0}=0_{2}$, $G_{0}=10 I_{2}$, $\Sigma _{\varepsilon }\sim IW_{2}\left( \nu _{0},S_{0}\right) $ with $\nu _{0}=10$ and $S_{0}=10 I_{2}$, $\omega \sim N\left( \omega _{0},O_{0} \right) $ with $\omega_{0}=0_{15}$, $O_{0}=I_{15}$, $\sigma _{Z}^{2}\sim IG\left( 50,10\right) $ and $\sigma _{u}^{2}\sim IG\left( 50,5\right)$. The prior distribution for $\sigma _{Z}^{2}$ is specified such that the prior mean ($0.2$) is close to the mean difference between first and second readings on transformed SBP ($0.024$). Similarly, the prior distribution for measurement error variance $\sigma _{u}^{2}$ is stipulated such that prior mean ($0.1$) is near the mean difference in variance from first and second readings on transformed SBP ($0.002$). Note that some of the parameters only appear in the measurement error model and the priors are used accordingly.
We first look at the results for the Bayesian estimation of SUR model presented in Table (ref) from 400,000 draws after a burn-in of 50,000 draws with a thinning factor of 100 (optimized following the approach in Owen-2017). The posterior estimates show that $\ln (age)$ is not statistically different from zero in both the $\ln (weight)$ and $HDL$ equations. Male indicator variable positively affects $\ln (weight)$, but negatively affects HDL. Height has a strong positive effect on $\ln (weight)$ and this is typically anticipated for all adults. Smoking daily or some days is negatively associated with $\ln (weight)$. This outcome is not surprising since smoking is well known to reduce appetite. On the contrary, smoking seems to have no significant effect on $HDL$. Number of hours of sedentary activities is positively associated with $\ln (weight)$ and negatively associated with HDL. The result confirms the generally held belief that being inactive increases weight and is negatively associated with good cholesterol. Sleep disorder is also known to be positively associated with weight gain and this is confirmed in our findings but it has no effect on HDL. LDL20T has a positive (negative) effect on $\ln (weight)$ (\emph{HDL}), which is expected since \emph{LDL} is commonly referred as `bad cholesterol' and is associated with weight gain. Transformed \emph{SBP} has a positive effect on $\ln (weight)$, but a negative effect on \emph{HDL} which is statistically different from zero at 90% probability level. As the first equation is a log-log specification, the coefficient of $\ln(SBP-50)$ is an elasticity. Then, a $10\%$ increase in the transformed \emph{SBP} leads to a $0.967\%$ growth of weight. The second equation is a semi-log specification, so the elasticity of \emph{HDL} relative to the transformed \emph{SBP} at the mean of the sample is: $0.1012/1.32 = 0.076$. A $10\%$ increase in the transformed $SBP$ leads to a $0.76\%$ growth of $HDL$. The estimated correlation coefficients of the residuals between the two equations is $\rho_{\varepsilon_1 \varepsilon_2} =-0.304$. Inefficiency factors are close to 1 suggesting a good mix of the draws and the Geweke's CD ($CD \sim N(0,1)$) confirms that a sufficiently large number of draws has been taken.
The results from the Bayesian estimation of SURME model (which accounts for the measurement error in the covariate SBP) is presented in the upper panel of Table (ref). A quick glance shows that the results for the covariates measured without error are similar to those in Table (ref), except for the intercept and the male indicator in the $\ln (weight)$ equation. The $\ln (height)$ coefficient in the $\ln (weight)$ equation is now slightly higher but its $ 95\%$ credible interval overlaps with that of the SUR model. Posterior estimates corresponding to transformed $SBP$ increase in both equations (from $0.097$ to $0.141$ in the $\ln(weight)$ equation and from $0.101$ to $0.152$ in the $HDL$ equation) leading to the following elasticities at the mean of the sample: $0.141$ and $0.115 (= 0.152/1.32)$, for weight and $HDL$, respectively. However, as the posterior standard errors become larger (from $0.027$ to $0.044$ in the $\ln(weight)$ equation and from $0.056$ to $0.092$ in the $HDL$ equation)---as in the Monte Carlo study---the $95\%$ HPDI of the posterior means of $\ln(SBP-50)$ overlap even if the distribution moves to the right when we go from SUR to SURME (see Figure (ref)). In particular, the $95\%$ HPDI of the posterior means of $\ln(SBP-50)$ are $\left[ 0.051;0.143\right] $ and $\left[ 0.009;0.194\right] $ in the $\ln (weight)$ and the $HDL$ equations, respectively, in the SUR model, and $\left[ 0.067;0.213\right]$ and $\left[ 0.005;0.304\right] $ in the $\ln (weight)$ and the $HDL$ equations in the SURME model. In the $HDL$ equation, the coefficient of SBP is positive but statistically equivalent to zero. Posterior estimate of measurement error variance is $0.029$, which leads to an estimated reliability ratio of about $ 59.78\%$ and a noise-to-true variance ratio of $67.27\%$. The posterior variances of the disturbances from $\ln (weight)$ and $HDL$ are close to those of SUR model and lead to an error correlation $\rho_{\varepsilon_1 \varepsilon_2} =-0.25$. Inefficiency factors are close to 1 and Geweke's CD confirms, for most parameters, that a sufficiently large number of draws has been taken.
The lower panel of Table (ref) presents the results for the exposure equation from the Bayesian estimation of SURME model. In the first equation, only three variables have a positive effect on $\ln(SBP-50)$: $\ln(age)$, LDL20T and $\ln (height)$ (only at the $10\%$ probability level). In the second equation, four variables have a positive effect on $\ln(SBP-50)$: $\ln(age)$, LDL20T, male and smokers (the last two variables are different from zero only at the $10\%$ probability level).
We next estimate the SURME model using the MFVB approach, which takes $1565$ cycles to get the maximum of the evidence lower bound $l$. The results, presented in the upper panel of Table (ref), shows that for the transformed SBP both the coefficient ($0.159$) and the probability interval ($\left[ 0.076;0.242\right]$) are similar to those obtained from the Bayesian estimation of SURME model. In the $HDL$ equation, the marginal effect at $0.186$ is higher compared to the Bayesian estimate, with a $95\%$ probability interval of $\left[0.03;0.34\right]$. Taking into account measurement errors using MFVB allows to get significantly larger and more accurate elasticities of weight ($0.159$) and $HDL$ ($0.1414= (0.186/1.32)$) for the transformed $SBP$ compared to the other method.
Posterior estimate of measurement error variance from the MFVB approach is $0.029$, which is similar to the Bayesian estimate, and leads to an estimated reliability ratio of about $60.2\%$ and to a noise-to-true variance ratio of $66.1\%$. The posterior variances of the disturbances from $\ln (weight)$ and $HDL$ are close to the Bayesian estimates and lead to the same correlation between the errors of the two equations $\rho_{\varepsilon_1 \varepsilon_2} =-0.25$. The lower panel of Table (ref) presents results for the exposure model estimated using with MFVB approximation method. In the first equation, four variables have a positive effect on $\ln(SBP-50)$: $\ln(age)$, smokers, $LDL20T$ and $\ln (height)$. Similarly, in the second equation four variables have a positive effect on $\ln(SBP-50)$: $\ln(age)$, $male$, $smokers$ and $LDL20T$. The exposure equation in the SURME model allows to define the implicit links between the true systolic blood pressure and the “risk factors” such as age, gender, smokers and “bad cholesterol” ($LDL20T$).
Figure (ref) gives the posterior densities of the parameter corresponding to $\ln(SBP-50)$ in the $\ln(weight)$ equation from the Bayesian estimation of SUR model, SURME model and the MFVB estimation of SURME model. We note a shift of the marginal effect of $\ln(SBP-50)$ on $\ln(weight)$ to the right of the distribution from a mode established around $0.097$ for SUR model to a mode of $0.141$ for SURME model but with a wider dispersion. The estimated probability density function (pdf) with MFVB is slightly to the right and centered around the mode ($0.159$) but with a surface under the curve globally equivalent to that of Bayesian estimation of SURME model. In Figure (ref), we observe similar shifts in the posterior density of the parameter corresponding to $\ln(SBP-50)$ in the $HDL$ equation, when we move from SUR ($0.101$) to SURME ($0.152$) or to MFVB ($0.187$) estimation of SURME model.
We note that similar to the Monte Carlo study, the MFVB approach to estimating the SURME model improves the results compared to those from Gibbs sampling. The results actually lend credibility to the proposed MFVB algorithm since coefficient estimates for variables which do not have measurement error are almost unaltered. However, when measurement error in $SBP$ is ignored as in the SUR model, the posterior estimates are underestimated relative to the MFVB estimates. So, accounting for measurement error potentially corrects or reduces the bias in parameter estimates Carroll-etal-2006.
The paper considers a SURME model (seemingly unrelated regression where some covariates have classical measurement error of the structural form) and introduces two novel estimation methods: a pure Bayesian algorithm based on MCMC and a second algorithm based on mean field variational Bayes approximation. The proposed algorithms use a prior distribution on measurement error variance to resolve identification issues in the model. In the MCMC estimation, Gibbs sampling is employed to sample the parameters from the conditional posterior distributions. While most of the conditional posterior densities have the standard form and are easily derived, the conditional posterior density for the true unobserved quantity associated with covariates having measurement error requires extensive attention to arrive at a manageable form. We also note that the proposed SURME model as explained is based on the structural form of measurement error, but the functional form of measurement error can be easily incorporated by introducing the distribution of the true unobserved quantity as a part of the subjective prior information. The expression for the joint and conditional posteriors will remain unchanged. However, estimating the SURME model using MCMC leads to high autocorrelation in the draws corresponding to the covariate measured with error. While this is easily dealt using thinning, the paper also proposes the MFVB approach as an alternative to get around the problem of high autocorrelation.
The proposed estimation algorithms are illustrated in multiple Monte Carlo simulation studies. While, the first set of 2 simulations (labeled Case I) investigate the effect on estimates by varying the variance of the true unobserved variable (for a fixed reliability ratio), the second set of 2 simulations examine the effect for a changing reliability ratio (for a fixed variance of the true unobserved variable). The results from all the simulations show that the Bayesian and MFVB estimation of SURME model reduce the biases to obtain satisfactory estimates as compared to estimates from SUR model. Moreover, the MFVB approach turns out as an excellent alternative to the MCMC and its poor mixing properties in the presence of latent variables. Besides, the MFVB approach has slightly better estimation accuracy and can be advantageous with large data sets.
The proposed models and techniques are also implemented in a health study where the two dependent variables, log of weight and high density lipoprotein (HDL), are regressed on a set of covariates measured without error and on systolic blood pressure (SBP) known to have measurement error. The model is estimated using the two algorithms and the results obtained reveal that the sign of the estimated coefficients are mostly consistent with what is typically found in the literature. Specifically, SBP has a positive effect on both $\ln(weight)$ and HDL, measurement error variance is small with an estimated reliability ratio of about $60\%$ and a noise-to-true variance ratio of $66\%$. To offer a baseline comparison, a SUR model that ignores measurement error in SBP is also estimated using Gibbs sampling. Comparing the results across models, we see that posterior estimates for covariates without measurement error are almost identical, but that of SBP is lower and hence underestimated both in the weight and the $HDL$ equations.
The combination of SUR and measurement error models is attractive and the proposed model can be generalized in several directions. One straightforward extension is the introduction of multiple covariates with measurement error in each SUR equation. However, the challenge here is to keep track of measurement errors arising from different covariates. The proposed SURME model can also be modified by introducing classical measurement error in the response variable or nonclassical measurement error models, where the errors may be correlated with the latent true values. Beyond the SUR models, these Bayesian approaches may be useful for measurement error in simultaneous equation models. We leave these possibilities for future research.
\pdfbookmark[1]{References}{unnumbered}