Quantile Regression with Multiple Proxy Variables

The natural cubic spline is a piecewise cubic polynomial function with continuous first and second derivatives at each knot, and is linear beyond the boundary knots (Green and Silverman, 1993). Suppose that $\tau_{1},\ldots,\tau_{N}$ are $N$ fixed knots covering a range of $x$ and $\bm{g}=\left(g_{1},\ldots,g_{N}\right)^{T}$ denote the values of the natural cubic spline $g_{p}(x)$ at knots $\tau_{1},\ldots\tau_{N}$, $g_{p}(\tau_{1}),\ldots g_{p}(\tau_{N})$. As a desirable property of the natural cubic spline, there is a unique natural cubic spline function $g_{p}(x)$ with knots $\tau_{1},\ldots,\tau_{N}$ satisfying $g(\tau_{i})=g_{i},\ i=1,\ldots,N$ for any given value $g_{1},\ldots,g_{N}$. Therefore, we can handle function $g(x)$ using its finite-length surrogate $\bm{g}$. In terms of Bayesian inference, we can model $g_{p}(x)$ by giving the prior to $\bm{g}$ and not $g_{p}(x)$. Following Green and Silverman (1993), the prior for $\bm{g}$ is defined by a multivariate normal distribution as follows:

$\displaystyle\pi\Big{(}\bm{g}|\lambda\Big{)}\propto\text{exp$\left(-\frac{1}{2}\lambda\bm{g}^{T}K\bm{g}\right)$},$

(5)

where $K$ is the $N\times N$ matrix with rank$(K)=N-2$, a function of the difference between the knots defined by (Eubank, 1999). $\lambda$ contributes to the smoothness of curve $g$ and has a standard conjugate gamma prior:

$$\pi\Big{(}\lambda\Big{)}=\frac{\lambda^{a_{\lambda}-1}\text{exp}(-\frac{\lambda}{b_{\lambda}})}{\Gamma(a_{\lambda})b_{\lambda}^{a_{\lambda}}},\text{ }\lambda>0.$$

The quadratic term $\bm{g}^{T}K\bm{g}$ in the exponential kernel is equivalent to the roughness penalty, $\int_{a}^{b}g_{p}^{{}^{\prime\prime}}(t)^{2}dt=\bm{g}^{T}K\bm{g}$ (Green and Silverman, 1993). This type of quadratic prior, $\exp\left(-\frac{1}{2}\lambda\int g_{p}^{{}^{\prime\prime}}(t)^{2}dt\right)$, which measures the complexity of the parameter, is a natural choice because it corresponds to the penalty term in the penalized maximum likelihood. With this prior, the posterior log density of the function $g_{p}(x_{i})$ is equal to the loss function in the regression context, with the roughness penalty added to the kernel of the log-likelihood function (Hoerl and Kennard, 1970; Green and Silverman, 1993; Tibshirani, 1996).

The final step in the Bayesian approach is defining the natural cubic spline function’s likelihood by changing the conventional polynomial part of the standard quantile regression likelihood to the natural cubic spline form (Yu and Moyeed, 2001; Thompson et al., 2010). The resulting likelihood takes the following form:

$\displaystyle L\Big{(}y\mid\bm{g},x\Big{)}=p^{n}(1-p)^{n}\text{exp}\left\{-\sum_{i=1}^{n}\rho_{p}(y_{i}-g_{p}(x_{i}))\right\}$

(6)

Notably, we explicitly specify $x$ in likelihood $L(y|\bm{g},x)$ for generalization, where $x$ is also an unknown variable.

Another general approach to spline fitting is a penalized spline or simply a P-spline. For a P-spline of degree $l$ with $N$ fixed knots, $g_{p}(x)$ is defined by $\bm{Z}(x)^{T}\bm{\beta}$ where $\bm{Z}(x)$ is $(N+l+1)$ vector composed of B-spline basis functions evaluated at observation $x$, and $\bm{\beta}$ is the coefficient of the basis functions (Eliers and Marx, 1996). A conventional basis is $\bm{Z}(x)=\left(1,x,\ldots,x^{l},(x-\tau_{1})^{l}_{+},\ldots,(x-\tau_{N})^{l}_{+}\right)^{T}$. Then, $\beta_{2+l},\ldots,\beta_{1+N+l}$ are the sizes of the jumps in the $l$th derivative of $g(x)$ at the knots.

Eliers and Marx (1996) suggested a roughness penalty based on differences of adjacent spline coefficient to guarantee sufficient smoothness. In Bayesian analysis, the prior of $\bm{\beta}$ replaces the roughness penalty term of the penalized likelihood as their stochastic analogs (Lang and Brezger, 2004). Assuming a first-order random walk for $\bm{\beta}$, that is,

$\displaystyle\pi\Big{(}\beta_{k}|\beta_{k-1},\lambda\Big{)}\propto N\Big{(}\beta_{k-1},\lambda\Big{)}$

the joint conditional distribution of $\bm{\beta}$ is

$\displaystyle\Big{(}\bm{\beta}|\lambda\Big{)}\propto\exp\left(-\frac{1}{2}\lambda\bm{\beta}^{T}K\bm{\beta}\right)$

where $K=\lambda R^{T}R$ is $(N+l+1)\times(N+l+1)$ penalty matrix with rank$(K)=N+l$ and $R$ is a first-order difference matrix (Lang and Brezger, 2004; Waldmann et al., 2013). This prior on $\bm{\beta}$ induces a prior on $g_{p}$ owing to the deterministic relationship between $g_{p}$ and $\bm{\beta}$, $g_{p}(x)=\bm{Z}(x)\bm{\beta}$. The precision parameter $\lambda$, again, contributes to the smoothness of curve $g$ and has a standard gamma prior:

$$\pi\Big{(}\lambda\Big{)}=\frac{\lambda^{a_{\lambda}-1}\text{exp}(-\frac{\lambda}{b_{\lambda}})}{\Gamma(a_{\lambda})b_{\lambda}^{a_{\lambda}}},\text{ }\lambda>0.$$

The P-spline function’s likelihood can be defined by changing the conventional polynomial part of the standard quantile regression likelihood. The same result is obtained in (6). However, unlike the natural cubic spline, Waldmann et al. (2013) suggested exploiting the stochastic representation of the likelihood for more efficient sampling (Kozumi and Kobayashi, 2011).

Suppose we have two proxies: one with an additive error and another with a quadratic relationship.

	$\displaystyle y_{i}$	$\displaystyle=g_{p}(x_{i})+e_{i},$
	$\displaystyle w_{1i}$	$\displaystyle=x_{i}+u_{1i},$
	$\displaystyle w_{2i}$	$\displaystyle=\alpha_{0}+\alpha_{1}x_{i}+\alpha_{2}x_{i}^{2}+u_{2i},$		(7)

where $u_{ki}\sim N(0,\sigma_{k}^{2}),\ k=1,2$. Here, parameter $\bm{\theta_{2}}$ for $h_{2}(x)$ is $\bm{\alpha}=(\alpha_{0},\alpha_{1},\alpha_{2})$. We use normal prior for $\pi(\alpha)\sim N(\mu_{\alpha},V_{\alpha})$. Following Thompson et al. (2010), we model a quantile function of covariate $x$ using the natural cubic spline, with $N$ evenly spaced fixed knots covering a range of $x$.

A Gibbs sampling algorithm for the quantile regression model is constructed by sampling each component of $\Theta$ from the full conditional distributions. Following Thompson et al. (2010)’s initialization, $\bm{g}^{(0)}$ is set as the posterior mean value of the quantile regression curve (Yu and Moyeed, 2001) at $\tau_{1},\dots\tau_{N}$, and $\lambda^{(0)}$ is obtained by applying generalized cross-validation of the usual smoothing spline (Green and Silverman, 1993). Additionally, $x^{(0)}$ is set as a multiple proxies $w_{1}$ because it is a more reliable proxy in the initialization step with no information regarding $\bm{\alpha}$.

One iteration of the Gibbs sampling algorithm at iteration $t$ is as follows:

1. 1.

   Generate candidate $\bm{g}^{*}$ from the multivariate normal distributions,

   $$\bm{g}^{*}|\bm{g}^{(t-1)}\sim\text{MVN}(\bm{g}^{(t-1)},\Sigma),$$

   and accept $\bm{g}^{*}$ with probability,

   $$r=\text{min}\left\{1,\frac{L(y|\bm{g}^{*},x^{(t-1)})\pi(\bm{g}^{*}|\lambda)q(\bm{g}^{(t-1)}|\bm{g}^{*})}{L(y|\bm{g}^{(t-1)},x^{(t-1)})\pi(\bm{g}^{(t-1)}|\lambda)q(\bm{g}^{*}|\bm{g}^{(t-1)})}\right\},$$

   where $q$ represents the proposal density function.

2. 2.

   Generate candidate $\lambda^{*}$ from the log-normal distribution,

   $$\eta^{*}\sim\text{N}(\text{log}(\lambda^{(t-1)}),\sigma_{\lambda}^{2}),$$

   where $\lambda^{*}=\exp(\eta^{*})$, and accept $\lambda^{*}$ with probability,

   $$r=\min\left\{1,\frac{\pi(\bm{g}^{(t)}|\lambda^{*})\pi(\lambda^{*})q(\lambda^{(t-1)}|\lambda^{*})}{\pi(\bm{g}^{(t)}|\lambda^{(t-1)})\pi(\lambda^{(t-1)})q(\lambda^{*}|\lambda^{(t-1)})}\right\}.$$

3. 3.

   Generate $x^{*}$ from the multivariate normal distribution,

   $$x^{*}|x^{(t-1)}\sim\text{MVN}(x^{(t-1)},\Sigma_{xx}),$$

   and accept $x^{*}$ with probability,

   $$r=\min\left\{1,\frac{L(y|\bm{g}^{(t)},x^{*})\pi(w_{1}|w_{2},x^{*})\pi(w_{2}|x^{*})\pi(x^{*}|\mu_{x}^{(t-1)},\left(\sigma_{x}^{2}\right)^{(t-1)})q(x^{(t-1)}|x^{*})}{L(y|\bm{g}^{(t)},x^{(t-1)})\pi(w_{1}|w_{2},x^{(t-1)})\pi(w_{2}|x^{(t-1)})\pi(x^{(t-1)}|\mu_{x}^{(t-1)},\left(\sigma_{x}^{2}\right)^{(t-1)})q(x^{*}|x^{(t-1)})}\right\}.$$

   From the independent assumption, this step can be separated as generating $x_{i}^{*}$ from the normal $x_{i}^{*}|x_{i}^{(t-1)}\sim N(x_{i}^{(t-1)},\sigma_{xx}^{2})$ with acceptance probability reduced to the term of the $i$-th data.

4. 4.

   Sample $\left(\sigma_{1}^{2}\right)^{(t)}\sim\text{Inv-Gamma}\left(\frac{n}{2}+a_{1},b_{1}+\frac{1}{2}\sum_{i=1}^{n}\left(w_{1i}-x_{i}^{(t)}\right)^{2}\right)$, where $\text{Inv-Gamma}(a,b)$ indicates the inverse gamma with shape parameter $a$ and scale parameter $b$, and $a_{1}$ and $b_{1}$ are the corresponding parameters for the prior of $\sigma_{1}^{2}$.

5. 5.

   Sample $\left(\sigma_{2}^{2}\right)^{(t)}\sim\text{Inv-Gamma}\left(\frac{n}{2}+a_{2},b_{2}+\frac{1}{2}\sum_{i=1}^{n}\left(w_{2i}-\left(\alpha_{0}^{(t-1)}+\alpha_{1}^{(t-1)}x_{i}^{(t)}\right)\right)^{2}\right)$ with $a_{2}$ and $b_{2}$ be the parameters for the prior of $\sigma_{2}^{2}$.

6. 6.

   Sample $\left(\sigma_{x}^{2}\right)^{(t)}\sim\text{Inv-Gamma}\left(\frac{n}{2}+a_{x},b_{x}+\frac{1}{2}\sum_{i=1}^{n}\left(x_{i}^{(t)}-\mu_{x}^{(t-1)}\right)^{2}\right)$ with $a_{x}$ and $b_{x}$ be parameters for the prior of $\sigma_{x}^{2}$.

7. 7.

   Sample $\mu_{x}^{(t)}\sim\text{N}(M_{*},V_{*})$, where

   $$V_{*}=\left(\frac{n}{\left(\sigma_{x}^{2}\right)^{(t)}}+\frac{1}{\sigma_{\mu}^{2}}\right)^{-1}\text{and }M_{*}=\left(\frac{\sum_{i=1}^{n}x_{i}^{(t)}+\frac{M_{\mu}}{\sigma_{\mu}^{2}}}{\left(\sigma_{x}^{2}\right)^{(t)}}\right)/\left(\frac{n}{\left(\sigma_{x}^{2}\right)^{(t)}}+\frac{1}{\sigma_{\mu}^{2}}\right)$$

   with $M_{\mu}$ and $\sigma_{\mu}^{2}$ the prior mean and variance for $\mu_{x}$.

8. 8.

   Sample $\alpha^{(t)}\sim\text{MVN}(\mu_{*},V_{*})$, where

   $$V_{*}=\left(\frac{X^{T}X}{\left(\sigma_{2}^{2}\right)^{(t)}}+V_{\alpha}^{-1}\right)^{-1}\text{and }\mu_{*}=\left(\frac{X^{T}X}{\left(\sigma_{2}^{2}\right)^{(t)}}+V_{\alpha}^{-1}\right)^{-1}\left(\frac{X^{T}W_{2}}{\left(\sigma_{2}^{2}\right)^{(t)}}+V_{\alpha}^{-1}\mu_{\alpha}\right)$$

   with $M_{\alpha}$ and $V_{\alpha}$ as the prior mean vector and covariance matrix for $\alpha$ and $X$ as the vector of $x_{i}^{(t)}$s, $i=1,\dots,n$.

Steps 1–3 require the Metropolis-Hastings algorithm, and the other steps can be sampled from the conjugate distribution. The inference regarding the unobserved regressor, quantile spline function, or regression coefficient is based on these posterior samples.

For most observed data in applications, prespecifying the relationship $h_{k}$ between the observed proxy and unobserved covariates in advance is challenging . Usually, a polynomial parametric relationship is applied with a certain purpose, such as interpretability. However, one might want a more flexible model with few prior assumptions regarding the relationship between the proxy and unobserved covariates. Provided that a benchmark variable is available for fixing the scale of unobserved $x$, the proposed method can adopt a nonlinear relationship for the other $h_{k}$.

Suppose we have two proxies—one with an additive error and another with a nonlinear relationship.

	$\displaystyle y_{i}$	$\displaystyle=g_{p}(x_{i})+e_{i},$
	$\displaystyle w_{1i}$	$\displaystyle=x_{i}+u_{1i},$
	$\displaystyle w_{2i}$	$\displaystyle=h_{2}(x_{i})+u_{2i},$		(8)

where $u_{ki}\sim N(0,\sigma_{k}^{2}),\ k=1,2$.

Following Waldmann et al. (2013), we used a stochastic representation of the likelihood $L\Big{(}y\mid\bm{g},x\Big{)}$, that is, $y\mid\bm{g},x\propto\mbox{N}(g_{p}(x)+As,B\frac{s}{\delta^{2}})$, where $A=\frac{1-2p}{p(1-p)},\ B=\frac{2}{p(1-p)}$ and $s\sim\text{Exp}(\delta^{2})$ with the conjugate gamma prior on $\delta^{2}\sim\mbox{GA}(a_{\delta},b_{\delta})$. This hierarchical representation of the likelihood in (6) enables efficient Gibbs sampling. For additional details, refer to Kozumi and Kobayashi (2011); Waldmann et al. (2013).

As we have two nonlinear functions $g_{p}$ and $h_{2}$ to be fitted, we use two P-splines. To discriminate the coefficients of the basis function for $g_{p}$ and $h_{2}$, we use subscripts $\bm{\beta}_{g}$ and $\bm{\beta}_{h}$. For a P-spline with $N$ fixed knots, $g_{p}$ is defined by $\bm{\beta}_{g}$, and for any realization of $\bm{\beta}_{g}$, there exists a corresponding realization $g_{p}(x)=\bm{Z}(x)\bm{\beta}_{g}$. As both splines $g_{p}$ and $h_{2}$ share covariate $x$, they share the same knots and penalty matrix $K$, which reduces the computation. Thereafter, the prior is as specified in Section 3.1.2; $\pi(\bm{\beta}_{g})\propto\exp\Big{(}-\frac{1}{2}\lambda_{g}\bm{\beta}_{g}^{T}K\bm{\beta}_{g}\Big{)}$, $\pi(\bm{\beta}_{h})\propto\exp\Big{(}-\frac{1}{2}\lambda_{h}\bm{\beta}_{h}^{T}K\bm{\beta}_{h}\Big{)}$, $\pi(\lambda_{g})\sim\mbox{GA}(a_{\lambda_{g}},b_{\lambda_{g}})$ and $\pi(\lambda_{h})\sim\mbox{GA}(a_{\lambda_{h}},b_{\lambda_{h}})$.

One iteration of the Gibbs sampling algorithm at iteration $t$ is as follows:

1. 1.

   Generate $x_{i}^{*}$ from the normal distribution, $x_{i}^{*}|x_{i}^{(t-1)}\sim\text{N}(x^{(t-1)},\sigma_{xx}^{2}),$ and accept $x_{i}^{*}$ with probability,

   	$\displaystyle r$	$\displaystyle=\min\Big{\{}1,\frac{N(y_{i}|Z(x_{i}^{*})^{T}\beta_{g}^{(t-1)}+As_{i}^{(t-1)},B\frac{s_{i}^{(t-1)}}{\delta^{2^{(t-1)}}})}{N(y_{i}|Z(x_{i}^{(t-1)})^{T}\beta_{g}^{(t-1)}+As_{i}^{(t-1)},B\frac{s_{i}^{(t-1)}}{\delta^{2^{(t-1)}}})}$
   		$\displaystyle\times\frac{\pi(w_{1i}|x_{i}^{*})\pi(w_{2i}|\beta_{h}^{(t-1)},x_{i}^{*})\pi(x_{i}^{*}|\mu_{x}^{(t-1)},\left(\sigma_{x}^{2}\right)^{(t-1)})q(x_{i}^{(t-1)}|x_{i}^{*})}{\pi(w_{1i}|x_{i}^{(t-1)})\pi(w_{2i}|\beta_{h}^{(t-1)},x_{i}^{(t-1)})\pi(x_{i}^{(t-1)}|\mu_{x}^{(t-1)},\left(\sigma_{x}^{2}\right)^{(t-1)})q(x_{i}^{*}|x_{i}^{(t-1)})}\Big{\}}$

2. 2.

   Sample $\bm{\beta_{g}}^{(t)}\sim\text{N}(M_{*},V_{*})$, where

   $$V_{*}=\left(\lambda_{g}^{(t-1)}K+\frac{\delta^{2^{(t-1)}}}{B}Z^{T}D^{-1}Z\right)^{-1}\text{and }M_{*}=V_{*}^{-1}\left(\frac{\delta^{2^{(t-1)}}}{B}Z^{T}D^{-1}\left({y}-A{s}^{(t-1)}\right)\right)$$

   with $D=\text{diag}(s_{1}^{(t-1)},\ldots s_{n}^{(t-1)})$ and $Z$ represents the design matrix with $Z(x_{i}^{(t)}),\ i=1,\ldots,n$.

3. 3.

   Sample $\lambda_{g}^{(t)}\sim\text{GA}(a_{\lambda}+0.5\text{rank}(\lambda_{g}^{(t-1)}K),b_{\lambda}+0.5\bm{\beta_{g}}^{(t)T}\lambda_{g}^{(t-1)}K\bm{\beta_{g}}^{(t)})$, where $\text{GA}(a,b)$ is gamma distribution with shape parameter $a$ and rate parameter $b$, and $a_{\lambda_{g}}$ and $b_{\lambda_{g}}$ are the corresponding parameters for the prior of $\lambda_{g}$.

4. 4.

   Sample $s_{i}^{{-1}^{(t)}}\sim\text{Inv-Gauss}\left(\frac{\sqrt{A^{2}+2B}}{y_{i}-Z(x_{i}^{(t)})^{T}\bm{\beta_{g}^{(t)}}},\frac{\delta^{2^{(t-1)}}(A^{2}+2B)}{B}\right)$, where $\text{Inv-Gauss}(a,b)$ is an inverse gaussian distribution with mean parameter $a$ and shape parameter $b$.

5. 5.

   Sample $\delta^{2^{(t)}}\sim GA(a_{\delta}+\frac{3n}{2},b_{\delta}+\frac{1}{2B}\sum_{i=1}^{n}s_{i}^{(t)^{-1}}\left(y_{i}-Z(x_{i}^{(t)})^{T}\beta_{g}^{(t)}-As_{i}^{(t)}\right)^{2}+\sum_{i=1}^{n}s_{i}^{(t)})$, where $a_{\delta}$ and $b_{\delta}$ are the corresponding parameters for the prior of $\delta^{2}$.

6. 6.

   Sample $\mu_{x}^{(t)}\sim\text{N}(M_{*},V_{*})$, where

   $$V_{*}=\left(\frac{n}{\left(\sigma_{x}^{2}\right)^{(t-1)}}+\frac{1}{\sigma_{\mu}^{2}}\right)^{-1}\text{and }M_{*}=\left(\frac{\sum_{i=1}^{n}x_{i}^{(t)}+\frac{M_{\mu}}{\sigma_{\mu}^{2}}}{\left(\sigma_{x}^{2}\right)^{(t-1)}}\right)/\left(\frac{n}{\left(\sigma_{x}^{2}\right)^{(t-1)}}+\frac{1}{\sigma_{\mu}^{2}}\right)$$

   with $M_{\mu}$ and $\sigma_{\mu}^{2}$ the prior mean and variance for $\mu_{x}$.

7. 7.

   Sample $\left(\sigma_{x}^{2}\right)^{(t)}\sim\text{Inv-Gamma}\left(\frac{n}{2}+a_{x},b_{x}+\frac{1}{2}\sum_{i=1}^{n}\left(x_{i}^{(t)}-\mu_{x}^{(t)}\right)^{2}\right)$ with $a_{x}$ and $b_{x}$ as the parameters for the prior of $\sigma_{x}^{2}$

8. 8.

   Sample $\left(\sigma_{1}^{2}\right)^{(t)}\sim\text{Inv-Gamma}\left(\frac{n}{2}+a_{1},b_{1}+\frac{1}{2}\sum_{i=1}^{n}\left(w_{1i}-x_{i}^{(t)}\right)^{2}\right)$ with $a_{1}$ and $b_{1}$ as the parameters for the prior of $\sigma_{1}^{2}$.

9. 9.

   Sample $\bm{\beta_{h}}^{(t)}\sim\text{N}(M_{*},V_{*})$, where

   $$V_{*}=\left(\lambda_{h}^{(t-1)}K+\frac{1}{\sigma_{2}^{2}}Z^{T}Z\right)^{-1}\text{and }M_{*}=V_{*}^{-1}\left(\frac{1}{\sigma_{2}^{2}}Z^{T}{w_{1}}\right)$$

10. 10.

    Sample $\lambda_{h}^{(t)}\sim\text{GA}(a_{\lambda_{h}}+0.5\text{rank}(\lambda_{h}^{(t-1)}K),b_{\lambda_{h}}+0.5\bm{\beta_{h}}^{(t)T}\lambda_{h}^{(t-1)}K\bm{\beta_{h}}^{(t)})$, where $a_{\lambda_{h}}$ and $b_{\lambda_{h}}$ are the corresponding parameters for the prior of $\lambda_{h}$.

11. 11.

    Sample $\left(\sigma_{2}^{2}\right)^{(t)}\sim\text{Inv-Gamma}\left(\frac{n}{2}+a_{2},b_{2}+\frac{1}{2}\sum_{i=1}^{n}\left(w_{2i}-Z(x_{i}^{(t)})^{T}\beta_{h}^{(t)}\right)^{2}\right)$ with $a_{2}$ and $b_{2}$ as the parameters for the prior of $\sigma_{2}^{2}$.

The inference regarding the unobserved regressor, quantile spline function, or regression coefficient is based on these posterior samples.