Heteroscedastic replicated measurement error models under asymmetric heavy-tailed distributions

Abstract

We propose a heteroscedastic replicated measurement error model based on the class of scale mixtures of skew-normal distributions, which allows the variances of measurement errors to vary across subjects. We develop EM algorithms to calculate maximum likelihood estimates for the model with or without equation error. An empirical Bayes approach is applied to estimate the true covariate and predict the response. Simulation studies show that the proposed models can provide reliable results and the inference is not unduly affected by outliers and distribution misspecification. The method has also been used to analyze a real data of plant root decomposition.

Appendix

1.1 Appendix A: Some heavy-tailed SMSN distributions

The pdf and the conditional moments of some heavy-tailed SMSN distributions are listed here:

1. The multivariate skew-t distribution $\text {ST}_m(\varvec{\mu }, \varvec{\Sigma }, \varvec{\lambda };\nu )$:

$\kappa (u)=1/u$, $U\sim \text {Gamma}(\nu /2,\nu /2)$ with $\nu > 0$. The pdf of $\varvec{Y}$ is given by

$$\begin{aligned} p(\varvec{y})=2t_m(\varvec{y}|\varvec{\mu },\varvec{\Sigma };\nu )T\bigg (\sqrt{\frac{m+\nu }{d+\nu }}A;\nu +m\bigg ), \end{aligned}$$

where $d=(\varvec{y}-\varvec{\mu })^{\top }\varvec{\Sigma }^{-1}(\varvec{y}-\varvec{\mu })$, $t_m(\cdot |\varvec{\mu }, \varvec{\Sigma };\nu )$ and $T(\cdot ;\nu )$ denotes the pdf of m-variate Student-t distribution and the cdf of standard univariate t distribution, respectively. When $\nu \rightarrow +\infty$, one obtains the SN distribution.

The conditional moments are given by

$$\begin{aligned}&u_r=\frac{p_0(\varvec{y})}{p(\varvec{y})}\frac{2^{r+1}\Gamma ((\nu +m+2r)/2)(\nu +d)^{-r}}{\Gamma ((\nu +m)/2)}T\bigg (\sqrt{\frac{m+\nu +2r}{d+\nu }}A;\nu +m+2r\bigg ),\\&\eta _r=\frac{p_0(\varvec{y})}{p(\varvec{y})}\frac{2^{(r+1)/2}\Gamma ((\nu +m+r)/2)}{\pi ^{1/2}\Gamma ((\nu +m)/2)}\frac{(\nu +d)^{(\nu +m)/2}}{(\nu +d+A^2)^{(\nu +m+r)/2}}. \end{aligned}$$

where , i.e. the pdf of the class of SMN distribution when $\varvec{\lambda }=\varvec{0}$.

2. The multivariate skew-slash distribution $\text {SS}_m(\varvec{\mu }, \varvec{\Sigma }, \varvec{\lambda };\nu )$:

$\kappa (u)=1/u$, $U\sim \text {Beta}(\nu , 1)$ with $0<u<1$ and $\nu > 0$. The pdf of $\varvec{Y}$ is given by

$$\begin{aligned} p(\varvec{y})=2\nu \int \nolimits _0^{1}u^{\nu -1}\varphi _m(\varvec{y}|\varvec{\mu },u^{-1}\varvec{\Sigma })\Phi (u^{1/2}A)\text {d}u. \end{aligned}$$

When $\nu \rightarrow +\infty$, it reduces to the SN one. The conditional moments are

$$\begin{aligned} u_r&=\frac{p_0(\varvec{y})}{p(\varvec{y})}\frac{2\Gamma ((2\nu +m+2r)/2)}{\Gamma ((2\nu +m)/2)}\Big (\frac{2}{d}\Big )^{r}\frac{P_1((2\nu +m+2r)/2,d/2)}{P_1((2\nu +m)/2,d/2)}\text {E}[\Phi (S^{1/2}A)],\\ \eta _r&=\frac{p_0(\varvec{y})}{p(\varvec{y})}\frac{2^{{(r+1)}/2}\Gamma ((2\nu +m+r)/2)}{\pi ^{1/2}\Gamma ((2\nu +m)/2)}\frac{d^{(2\nu +m)/2}}{(d+A^2)^{(2\nu +m+r)/2}}\\&\quad \ \times \frac{P_1((2\nu +m+r)/2,(d+A^2)/2)}{P_1((2\nu +m)/2,d/2)}, \end{aligned}$$

where $S\sim \text {Gamma}((2\nu +m+2r)/2,d/2)\text {I}_{(0,1)}$ and $P_x(a,b)$ denotes the cdf of the $\text {Gamma}(a,b)$ distribution evaluated at x.

3. The multivariate skew-contaminated normal distribution $\text {SCN}_m(\varvec{\mu }, \varvec{\Sigma }, \varvec{\lambda };\nu ,\gamma )$:

When $\kappa (u)=1/u$ and U follows a discrete random probability function $h(u;\nu ,\gamma )=\nu \text {I}_{(u=\gamma )}+(1-\nu ) \text {I}_{(u=1)}$ with parameter $0<\nu<1, 0<\gamma \leqslant 1$, one get the multivariate skew-contaminated normal distribution with the pdf as

$$\begin{aligned} p(\varvec{y})=2\{\nu \varphi _m(\varvec{y}|\varvec{\mu },\gamma ^{-1}\varvec{\Sigma })\Phi (\gamma ^{1/2}A)+(1-\nu )\varphi _m(\varvec{y}|\varvec{\mu },\varvec{\Sigma })\Phi (A)\}. \end{aligned}$$

It reduces to the SN distribution when $\gamma =1$. In this case, we have

$$\begin{aligned}&u_r=\frac{2}{p(\varvec{y})}\{\nu \gamma ^{r}\varphi _m(\varvec{y}|\varvec{\mu },\gamma ^{-1}\varvec{\Sigma })\Phi (\gamma ^{1/2}A)+(1-\nu )\varphi _m(\varvec{y}|\varvec{\mu },\varvec{\Sigma })\Phi (A)\},\\&\eta _r=\frac{2}{p(\varvec{y})}\{\nu \gamma ^{r/2}\varphi _m(\varvec{y}|\varvec{\mu },\gamma ^{-1}\varvec{\Sigma })\phi (\gamma ^{1/2}A)+(1-\nu )\varphi _m(\varvec{y}|\varvec{\mu },\varvec{\Sigma })\phi (A)\}. \end{aligned}$$

1.2 Appendix B: EM algorithm for SMSN-HRME model without equation error

Denoting the complete data set of model (4) without equation error by $\varvec{Z}_c=\{\varvec{Z}_{ct}=(\varvec{Z}_t^{\top }, x_t, u_t, v_t)^{\top }|t=1,\ldots ,n\}$, an equivalent form of structure (5) is given by

$$\begin{aligned} \begin{array}{c} \varvec{Z}_t|x_t,U_t=u_t\mathop {\sim }\limits ^{ind}\text {N}_{m_t}(\varvec{a}_t+x_t\varvec{b}_t, \kappa (u_t)\mathbf D (\varvec{\phi }_{t})),\\ x_t|U_t=u_t,V_t=v_t\mathop {\sim }\limits ^{ind}\text {N}(\mu _x+\tau _x v_t, \kappa (u_t)\gamma _x),\\ V_t|U_t=u_t\mathop {\sim }\limits ^{ind}\text {HN}(0,\kappa (u_t)),\\ U_t\mathop {\sim }\limits ^{ind}\text {H}(u_t; \varvec{\nu }),\ t=1,\ldots , n, \end{array} \end{aligned}$$

(11)

The complete log-likelihood function based on $\varvec{Z}_c$ omitting items unrelated with $\varvec{\theta }$ can be written as

$$\begin{aligned} l(\varvec{\theta }\big |\varvec{Z}_c)=\sum _{t=1}^n {\big (l_{x_t|v_t,u_t}+l_{\varvec{Z}_t|x_t,u_t}\big )}, \end{aligned}$$

where,

$$\begin{aligned}&l_{x_t|v_t,u_t}=-\frac{1}{2}\log (\gamma _x)-\frac{1}{2}\gamma _x^{-1}\kappa ^{-1}(u_t)\big (x_t-\mu _x-\tau _x v_t\big )^2,\\&l_{\varvec{Z}_t|x_t,u_t} =-\frac{1}{2}\kappa ^{-1}(u_t){(\varvec{Z}_t-\varvec{a}_t-\varvec{b}_tx_t)}^{\top }\mathbf D ^{-1}(\varvec{\phi }_{t}){(\varvec{Z}_t-\varvec{a}_t-\varvec{b}_tx_t)}. \end{aligned}$$

The EM algorithm is listed as follows.

E-step: Let $\varvec{\theta }^{(k)}$ be the estimates of $\varvec{\theta }$ at the k-th iteration. By calculating the conditional expectation $\text {E}\big [l(\varvec{\theta }|\varvec{Z}_c)\big |\widehat{\varvec{\theta }}^{(k)}, \varvec{Z}\big ]$, we get the Q-function

$$\begin{aligned} Q(\varvec{\theta }|\widehat{\varvec{\theta }}^{(k)})=\sum _{t=1}^n\big [Q_{1t}(\varvec{\theta }|\widehat{\varvec{\theta }}^{(k)})+Q_{2t}(\varvec{\theta }|\widehat{\varvec{\theta }}^{(k)})\big ], \end{aligned}$$

with

$$\begin{aligned} \begin{aligned} Q_{1t}(\varvec{\theta }|\widehat{\varvec{\theta }}^{(k)})&=-\frac{1}{2}\log (\gamma _x)-\frac{1}{2}\gamma _x^{-1}\big (\widehat{ux^2}_t^{(k)}+\mu _x^2\widehat{u}_t^{(k)}+\tau _x^2\widehat{uv^2}_t^{(k)}\\&\quad -\,2\mu _x\widehat{ux}_t^{(k)}-2\tau _x\widehat{uxv}_t^{(k)}+2\mu _x\tau _x\widehat{uv}_t^{(k)}\big ), \end{aligned} \end{aligned}$$

and

$$\begin{aligned} \begin{aligned} Q_{2t}(\varvec{\theta }|\widehat{\varvec{\theta }}^{(k)})&=-\,\frac{1}{2}\big [\widehat{u}_t^{(k)}(\varvec{Z}_t-\varvec{a}_t)^{\top }\mathbf D ^{-1}(\varvec{\phi }_{t})(\varvec{Z}_t-\varvec{a}_t)\\&\quad -\,2\widehat{ux}_t^{(k)}(\varvec{Z}_t-\varvec{a}_t)^{\top }\mathbf D ^{-1}(\varvec{\phi }_{t})\varvec{b}_t +\widehat{ux^2}_t^{(k)}\varvec{b}_t^{\top }\mathbf D ^{-1}(\varvec{\phi }_{t})\varvec{b}_t\big ], \end{aligned} \end{aligned}$$

where $\widehat{ux}_t^{(k)}=\text {E}[\kappa ^{-1}(U_t)x_t|\widehat{\varvec{\theta }}^{(k)},\varvec{Z}_t]$, $\widehat{ux^2}_t^{(k)}=\text {E}[\kappa ^{-1}(U_t)x_t^2|\widehat{\varvec{\theta }}^{(k)},\varvec{Z}_t]$, $\widehat{uxv}_t^{(k)}=\text {E}[\kappa ^{-1}(U_t)x_t V_t|\widehat{\varvec{\theta }}^{(k)},\varvec{Z}_t]$, and their computational expressions are as follows

$$\begin{aligned} \widehat{ux}_t^{(k)}&=\widehat{u}_t^{(k)}\widehat{r}_t^{(k)}+\widehat{uv}_t^{(k)}\widehat{s}_t^{(k)},\\ \widehat{ux^2}_t^{(k)}&=\widehat{u}_t^{(k)}{\widehat{r}_t}^{2(k)}+2\widehat{uv}_t^{(k)}\widehat{r}_t^{(k)}\widehat{s}_t^{(k)} +\widehat{uv^2}_t^{(k)}{\widehat{s}}_t^{2(k)}+\widehat{\gamma }_x^{(k)}/\widehat{c}_{1t}^{(k)},\\ \widehat{uxv}_t^{(k)}&=\widehat{uv}_t^{(k)}\widehat{r}_t^{(k)}+\widehat{uv^2}_t^{(k)}\widehat{s}_t^{(k)}, \end{aligned}$$

with $\widehat{r}_t=\widehat{\mu }_x+\widehat{\gamma }_x\widehat{a}_t/\widehat{c}_{1t}$, $\widehat{s}_t=\widehat{\tau }_x/\widehat{c}_{1t}$. Here $a_t=(\varvec{Z}_t-\varvec{\mu }_t)^{\top }\mathbf D ^{-1}(\varvec{\phi }_{t})\varvec{b}_t$, $c_t=1+\phi _x\varvec{b}_t^{\top }\mathbf D ^{-1}(\varvec{\phi }_{t})\varvec{b}_t$, $c_{1t}=1+\gamma _x\varvec{b}_t^{\top }\mathbf D ^{-1}(\varvec{\phi }_{t})\varvec{b}_t$.

M-step: Maximizing $Q(\varvec{\theta }|\widehat{\varvec{\theta }}^{(k)})$ with respect to $\varvec{\theta }$, we achieve the updated estimates $\widehat{\varvec{\theta }}^{(k+1)}$ by the following iterative equations:

$$\begin{aligned} \widehat{\mu }_x^{(k+1)}&=\frac{\sum _{t=1}^n \widehat{ux}_t^{(k)}\sum _{t=1}^n \widehat{uv^2}_t^{(k)}-\sum _{t=1}^n \widehat{uv}_t^{(k)}\sum _{t=1}^n \widehat{uxv}_t^{(k)}}{\sum _{t=1}^n \widehat{u}_t^{(k)}\sum _{t=1}^n \widehat{uv^2}_t^{(k)}-\big (\sum _{t=1}^n \widehat{uv}_t^{(k)}\big )^2},\\ \widehat{\tau }_x^{(k+1)}&=\frac{\sum _{t=1}^n \widehat{u}_t^{(k)}\sum _{t=1}^n \widehat{uxv}_t^{(k)}-\sum _{t=1}^n \widehat{uv}_t^{(k)}\sum _{t=1}^n \widehat{ux}_t^{(k)}}{\sum _{t=1}^n \widehat{u}_t^{(k)}\sum _{t=1}^n \widehat{uv^2}_t^{(k)}-\big (\sum _{t=1}^n \widehat{uv}_t^{(k)}\big )^2},\\ \widehat{\gamma }_x^{(k+1)}&=\frac{1}{n}\sum _{t=1}^n \bigg (\widehat{ux^2}_t^{(k)}+(\widehat{\mu }_x^{(k+1)})^2\widehat{u}_t^{(k)}+(\widehat{\tau }_x^{(k+1)})^2\widehat{uv^2}_t^{(k)}\\&-\,2\widehat{\mu }_x^{(k+1)}\widehat{ux}_t^{(k)}-2\widehat{\tau }_x^{(k+1)}\widehat{uxv}_t^{(k)} +2\widehat{\mu }_x^{(k+1)}\widehat{\tau }_x^{(k+1)}\widehat{uv}_t^{(k)}\bigg ), \end{aligned}$$

$$\begin{aligned} \widehat{\alpha }^{(k+1)}&=\frac{\sum _{t=1}^n \big (\widehat{u}_t^{(k)}\bar{Y}_tq_t/\phi _{\varepsilon _t}\big )\sum _{t=1}^n \big (\widehat{ux^2}_t^{(k)}q_t/\phi _{\varepsilon _t}\big )-\sum _{t=1}^n \big (\widehat{ux}_t^{(k)}q_t/\phi _{\varepsilon _t}\big )\sum _{t=1}^n \big (\widehat{ux}_t^{(k)}\bar{Y}_tq_t/\phi _{\varepsilon _t}\big )}{\sum _{t=1}^n \big (\widehat{u}_t^{(k)}q_t/\phi _{\varepsilon _t}\big )\sum _{t=1}^n \big (\widehat{ux^2}_t^{(k)}q_t/\phi _{\varepsilon _t}\big )-\big (\sum _{t=1}^n \widehat{ux}_t^{(k)}q_t/\phi _{\varepsilon _t}\big )^2},\\ \widehat{\beta }^{(k+1)}&=\frac{\sum _{t=1}^n \big (\widehat{ux}_t^{(k)}\bar{Y}_tq_t/\phi _{\varepsilon _t}\big )\sum _{t=1}^n (\widehat{u}_t^{(k)}q_t/\phi _{\varepsilon _t})-\sum _{t=1}^n \bigg (\widehat{ux}_t^{(k)}q_t/\phi _{\varepsilon _t}\bigg )\sum _{t=1}^n \bigg (\widehat{u}_t^{(k)}\bar{Y}_tq_t/\phi _{\varepsilon _t}\bigg )}{\sum _{t=1}^n \bigg (\widehat{u}_t^{(k)}q_t/\phi _{\varepsilon _t}\bigg )\sum _{t=1}^n \bigg (\widehat{ux^2}_t^{(k)}q_t/\phi _{\varepsilon _t}\bigg )-\bigg (\sum _{t=1}^n \widehat{ux}_t^{(k)}q_t/\phi _{\varepsilon _t}\bigg )^2}, \end{aligned}$$

where $\bar{X}_t=\frac{1}{p_t}\sum _{i=1}^{p_t} X_t^{(i)}$, $\bar{Y}_t=\frac{1}{q_t}\sum _{j=1}^{q_t} Y_t^{(j)}$. Note that the estimators $\widehat{\lambda }_x$ and $\widehat{\phi }_x$ can be inferred from the one-to-one transformation $\widehat{\lambda }_x=\widehat{\tau }_x/\sqrt{\widehat{\gamma }_x}$ and $\widehat{\phi }_x=\widehat{\gamma }_x+\widehat{\tau }_x^2$.

1.3 Appendix C: The related derivatives

By direct calculations, we have the first derivatives of $d_t$ and $A_t$ respect to $\varvec{\theta }$ as follows:

$$\begin{aligned} \frac{\partial d_t}{\partial \theta _i}&=-2{(\varvec{Z}_t-\varvec{\mu }_t)}^{\top }\varvec{\Sigma }_t^{-1}\frac{\partial \varvec{\mu }_t}{\partial \theta _i} -{(\varvec{Z}_t-\varvec{\mu }_t)}^\top \varvec{\Sigma }_t^{-1}\frac{\partial \varvec{\Sigma }_t}{\partial \theta _i}\varvec{\Sigma }_t^{-1}(\varvec{Z}_t-\varvec{\mu }_t),\\ \frac{\partial A_t}{\partial \theta _i}&=\bigg (\frac{\partial \psi _t}{\partial \theta _i}\varvec{b}_t^{\top } +\psi _t\frac{\partial \varvec{b}_t^{\top }}{\partial \theta _i}-\psi _t\varvec{b}_t^{\top }\varvec{\Sigma }_t^{-1} \frac{\partial \varvec{\Sigma }_t}{\partial \theta _i}\bigg )\varvec{\Sigma }_t^{-1}(\varvec{Z}_t-\varvec{\mu }_t) -\psi _t\varvec{b}_t^{\top }\varvec{\Sigma }_t^{-1}\frac{\partial \varvec{\mu }_t}{\partial \theta _i}, \end{aligned}$$

where $\psi _t=\frac{\lambda _x\phi _x}{\sqrt{\phi _x+\lambda _x^2\Lambda _{xt}}}$.

In addition, we also need to calculate the following derivatives:

$$\begin{aligned} \frac{\partial \varvec{\mu }_t}{\partial \mu _x}&=\varvec{b}_t,\quad \frac{\partial \varvec{\mu }_t}{\partial \alpha }=\varvec{c}_t,\quad \frac{\partial \varvec{\mu }_t}{\partial \beta }=\mu _x\varvec{c}_t,\quad \frac{\partial \varvec{b}_t}{\partial \beta }=\varvec{c}_t,\\ \frac{\partial \varvec{\Sigma }_t}{\partial \beta }&=\phi _x(\varvec{c}_t\varvec{b}_t^{\top }+\varvec{b}_t\varvec{c}_t^{\top }),\quad \frac{\partial \varvec{\Sigma }_t}{\partial \phi _x}=\varvec{b}_t\varvec{b}_t^{\top },\quad \frac{\partial \varvec{\Sigma }_t}{\partial \phi _e}=\varvec{c}_t\varvec{c}_t^{\top },\\ \frac{\partial \psi _t}{\partial \beta }&=q_t\beta c_t^{-2}\psi _t^3/(\phi _{\varepsilon _t}+q_t\phi _e),\quad \frac{\partial \psi _t}{\partial \lambda _x}=\phi _x^2(\phi _x+\lambda _x^2\Lambda _{xt})^{-3/2},\\ \frac{\partial \psi _t}{\partial \phi _x}&=\frac{1}{2}\psi _t\phi _x^{-1}+\frac{1}{2}\psi _t^3\phi _x^{-1}c_t^{-2}\varvec{b}_t^{\top }\varvec{\Sigma }_{1t}^{-1}\varvec{b}_t,\\ \frac{\partial \psi _t}{\partial \phi _e}&=-\frac{1}{2}q_t^2\psi _t^3c_t^{-2}\beta ^2{(\phi _{\varepsilon _t}+q_t\phi _e)}^{-2}. \end{aligned}$$