Basic Setting

• Predictions for small domains are our primary interest.

  • Standard survey estimators are often unreliable given small sample sizes. (Rao, Molina, 2015)

  • Use model to incorporate auxiliary information (known for the whole population).

• Data motivation (CEAP RUSLE2)

  • skewed data with heavy zeros

  • potential correlation between the positive part and the binary part

Data Model

Suppose response variable $y_{ij}^{\ast },i=1,..,D;j=1,..,N_i$ satisfies $$y_{ij}^{\ast }=y_{ij}{\delta }_{ij}$$

• Positive part: $$\log \left(y_{ij}\right)={\beta }_0+z_{1,ij}^{\prime }{\beta }_1+u_i+e_{ij}$$ where $e_{ij}\stackrel{i.i.d}{\sim }N(0,{\sigma }_e^2)$.

• Binary part: $\text{Pr}({\delta }_{ij}=1)=p_{ij}$, $$g\left(p_{ij}\right)={\alpha }_0+z_{2,ij}^{\prime }{\alpha }_1+b_i$$ where $g(\cdot )$ is the parametric link function to be specified.

Data Model

• Lognormal extension

  • Mathematically tractable and computationally simple.

  • SAE lognormal model - Berg and Chandra (2013).

  • SAE zero-inflated lognormal - Lyu (2018).

• Correlated random effects

  • $$\left(\begin{array}{c}{u}_i\\ b_i\end{array}\right)\stackrel{iid}{\sim }\text{BVN}(0,{\Sigma }_{ub}),{\Sigma }_{ub}=\left(\begin{array}{cc}{\sigma }_u^2& \rho {\sigma }_u{\sigma }_b\\ \rho {\sigma }_u{\sigma }_b& {\sigma }_b^2\end{array}\right)$$

  • "plug-in" approach of Chandra, Chambers (2016) => implicitly assume independence

  • Bayesian approach of Pfeffermann et al (2018) => requires specifying prior distributions

Empirical Bayes general

The minimum mean squared error (MMSE) predictor of ${\zeta }_i=\frac{1}{N_i}{\sum }_{j=1}^{N_i}{y}_{ij}^{\ast }$ is $${\hat{\zeta }}_i^{MMSE}=\frac{1}{N_i}\{\sum _{j\in s_i}{y}_{ij}^{\ast }+\sum _{j\in {\overline{s}}_i}{\hat{y}}_{ij}^{\ast MMSE}\}$$ where ${\hat{y}}_{ij}^{\ast MMSE}=E_{\theta }\left(y_{ij}^{\ast }|y_{s_i}^{\ast }\right)$

Let $\hat{\theta }$ be a consistent estimator of model parameter $\theta$. $${\hat{\zeta }}_i^{EB}=\frac{1}{N_i}\{\sum _{j\in s_i}{y}_{ij}^{\ast }+\sum _{j\in {\overline{s}}_i}{\hat{y}}_{ij}^{\ast EB}\}$$ where ${\hat{y}}_{ij}^{\ast EB}={\hat{y}}_{ij}^{\ast MMSE}\left(\hat{\theta }\right)$

Conditional Distribution

Result 1:

$$u_i|b_i,y_{s_i}^{\ast }\sim N({\tilde{\mu }}_{u_i},{\tilde{\sigma }}_{u_i}^2)$$ where $${\tilde{\mu }}_{u_i}={\gamma }_i{\overline{\tilde{r}}}_i+(1-{\gamma }_i)[\rho \frac{{\sigma }_u}{{\sigma }_b}{b}_i],{\tilde{\sigma }}_{u_i}^2={\gamma }_i{\sigma }_e^2/{\tilde{n}}_i$$ and ${\gamma }_i=(1-{\rho }^2){\sigma }_u^2\left[\right(1-{\rho }^2){\sigma }_u^2+{\sigma }_e^2/{\tilde{n}}_i{]}^{-1},{\tilde{n}}_i={\sum }_{j\in s_i}{\delta }_{ij}$, $${\overline{\tilde{r}}}_i={\tilde{n}}_i^{-1}\sum _{j\in s_i}{\delta }_{ij}\left\{\log \right(y_{ij})-{\beta }_0-z_{1,ij}^{\prime }{\beta }_1\}$$

Remark: when $\rho =0$, ${\gamma }_i={\sigma }_u^2/\left({\sigma }_u^2+{\sigma }_e^2/{\tilde{n}}_i\right)$, $u_i|y_{s_i}^{\ast }\sim N({\gamma }_i{\overline{\tilde{r}}}_i,{\gamma }_i{\sigma }_e^2/{\tilde{n}}_i)$.

Conditional Distribution (Cont'd)

Result 2:

$$f\left(b_i|y_{s_i}^{\ast }\right)\propto {\pi }_{s_i}\left(b_i\right)\varphi \left(\frac{b_i-m_i}{\sqrt{v_i}}\right)$$ where $\varphi (\cdot )$ is the probability density function of standard normal distribution, $${\pi }_{s_i}\left(b_i\right)=\prod _{j\in s_i}[p_{ij}{\delta }_{ij}+(1-p_{ij})(1-{\delta }_{ij})]$$ and $p_{ij}\left(b_i\right)=g^{-1}({\alpha }_0+{\alpha }_1{z}_{ij}+b_i)$. $$(m_i,v_i)=(\frac{\rho {\sigma }_u{\sigma }_b}{{\sigma }_u^2+{\sigma }_e^2/{\tilde{n}}_i}{\overline{\tilde{r}}}_i,\frac{1-{\rho }^2}{1-(1-{\gamma }_i){\rho }^2}{\sigma }_b^2).$$ Remark: when $\rho =0$, $(m_i,v_i)=(0,{\sigma }_b^2)$.

Conditional Distribution (Cont'd)

Result 3:

$$E(y_{ij}^{\ast }|y_{s_i}^{\ast },\theta )=h_{ij}\left(\theta \right)E\left(p_{ij}\right(b_i\left)\eta \left(b_i\right)|y_{s_i}^{\ast }\right)$$ where $h_{ij}\left(\theta \right)=\exp ({\beta }_0+x_{ij}^{\prime }{\beta }_1+{\gamma }_i{\overline{\tilde{r}}}_i+{\tilde{\sigma }}_{u_i}^2/2+{\sigma }_e^2/2)$ and $\eta \left(b_i\right)=\exp \left(\right(1-{\gamma }_i\left)\rho {\sigma }_u/{\sigma }_b{b}_i\right).$

👉 ${\hat{y}}_{ij}^{\ast EB}=E(y_{ij}^{\ast }|y_{s_i}^{\ast },\hat{\theta })$

Remark: The proposed predictor can accommodate any parametrized link (transformation) function for the positive part and the binary part.

Parameter Estimation

Given $\rho \ne 0$, how to get a consistent estimator for model parameter $\theta =({\beta }^{\prime },{\alpha }^{\prime },{\sigma }_u^2,{\sigma }_b^2,{\sigma }_e^2,\rho {)}^{\prime }$? 😣

Full likelihood function: $$\begin{array}{ccc}{L}_i\left(\theta \right)& =& \frac{(1-{\gamma }_i{)}^{1/2}(v_i/{\sigma }_b^2{)}^{1/2}}{(2\pi {\sigma }_e^2{)}^{{\tilde{n}}_i/2}}\exp (\frac{{\gamma }_i{\overline{\tilde{r}}}_i^2}{2{\sigma }_e^2/{\tilde{n}}_i}+\frac{m_i^2}{2v_i}-\frac{{\sum }_j{\tilde{r}}_{ij}^2}{2{\sigma }_e^2})\\ & & \int {\pi }_{s_i}\left(b_i\right)\frac{1}{\sqrt{v_i}}\varphi \left(\frac{b_i-m_i}{\sqrt{v_i}}\right)d{b}_i\end{array}$$ Remark:

• A good starting value and the gradient vector $\partial L/\partial \theta$ help optim find the optimizer much faster. ✌️

• Profiling out $\rho$ takes much longer time to find the MLE.

MSE Estimator

Formally, the MSE of an EB predictor is $$MSE\left({\hat{\zeta }}_i^{EB}\right)=M_{1i}\left(\theta \right)+M_{2i}\left(\theta \right)$$ where $$M_{1i}\left(\theta \right)=E({\hat{\zeta }}_i^{MMSE}-{\zeta }_i{)}^2=O(1),\text{as}D\to \infty$$ and $$M_{2i}\left(\theta \right)=E({\hat{\zeta }}_i^{EB}-{\hat{\zeta }}_i^{MMSE}{)}^2=o(1),\text{as}D\to \infty$$ Recall that $${\hat{\zeta }}_i^{MMSE}=E[{\zeta }_i|y_{s_i}^{\ast },\theta ],{\hat{\zeta }}_i^{EB}=E[{\zeta }_i|y_{s_i}^{\ast },\hat{\theta }]$$

MSE Estimator (Cont'd)

So the MSE of the optimal predictor is $$M_{1i}\left(\theta \right)=E\left(E\right[{\zeta }_i|y_{s_i}^{\ast },\theta ]-{\zeta }_i{)}^2=E[V({\zeta }_i|y_{s_i}^{\ast },\theta )]$$ We propose a One-Step MSE estimator defined by $$V({\zeta }_i|y_{s_i}^{\ast },\hat{\theta })=\frac{1}{N_i^2}\sum _{j\in {\overline{s}}_i}\sum _{k\in {\overline{s}}_i}[E\{y_{ij}^{\ast }{y}_{ik}^{\ast }|y_{s_i}^{\ast },\hat{\theta }\}-E\{y_{ij}^{\ast }|y_{s_i}^{\ast },\hat{\theta }\}E\{y_{ik}^{\ast }|y_{s_i}^{\ast },\hat{\theta }\}]$$ where $$E\left\{y_{ij}^{\ast }{y}_{ik}^{\ast }|y_{s_i}^{\ast }\right\}=h_{ij}{h}_{ik}\exp \left({\tilde{\sigma }}_{u_i}^2\right)E\left[p_{ij}{p}_{ik}\eta \left(2b_i\right)|y_{s_i}^{\ast }\right],j\ne k$$ and $$E\left\{y_{ij}^{\ast }{y}_{ik}^{\ast }|y_{s_i}^{\ast }\right\}=h_{ij}^2\exp ({\tilde{\sigma }}_{u_i}^2+{\sigma }_e^2)E\left[p_{ij}\eta \left(2b_i\right)|y_{s_i}^{\ast }\right],j=k$$

Alternate MSE Estimators

• Denote original parameter estimate $\hat{\theta }$, bootstrap population $Y^{\ast \left(b\right)}$, bootstrap sample $y^{\ast \left(b\right)}$, parameter estimate ${\hat{\theta }}^{\left(b\right)}$

• For $b=1,\dots ,B$, we obtain $$Y^{\ast \left(b\right)}⟶{\overline{y}}_{N_i}^{\ast \left(b\right)}$$ $$y^{\ast \left(b\right)},\hat{\theta }⟶{\hat{\overline{y}}}_{N_i}^{\ast \left(b\right)MMSE}$$ $$y^{\ast \left(b\right)},{\hat{\theta }}^{\left(b\right)}⟶{\hat{\overline{y}}}_{N_i}^{\ast \left(b\right)EB}$$

1. ${\widehat{\text{MSE}}}_i^{\text{Boot}}=B^{-1}{\sum }_{b=1}^B({\hat{\overline{y}}}_{N_i}^{\ast \left(b\right)\text{EB}}-{\overline{y}}_{N_i}^{\ast \left(b\right)}{)}^2$

2. ${\widehat{\text{MSE}}}_i^{\text{Semi-Boot}}={\hat{M}}_{1i}+{\hat{M}}_{2i}^{\text{Boot}}$ where a bootstrap estimator of $M_{2i}$ is defined by ${\hat{M}}_{2i}^{\text{Boot}}=B^{-1}{\sum }_{b=1}^B({\hat{\overline{y}}}_{N_i}^{\ast \left(b\right)\text{EB}}-{\hat{\overline{y}}}_{N_i}^{\ast \left(b\right)\text{MMSE}}{)}^2$

Simulation Setting

• Number of areas: $D=60$

• Sample rate: $(N_i,n_i)=(71,5),(143,10),(286,20)$ for every 20 areas so that $(N,n)=(10000,700)$

• logit link for the binary part

• One-dimensional covariates: $z_{ij}\sim N(4.45,0.055)$

• $\beta =(-13,2{)}^{\prime },\alpha =(-20,5{)}^{\prime }$, $({\sigma }_u^2,{\sigma }_e^2,{\sigma }_b^2)=(0.22,1.23,0.52)$

• $\rho \in \{-0.9,-0.6,-0.3,0,0.3,0.6,0.9\}$

• Compare the proposed EB predictor with the EB(0) predictor ${\hat{\zeta }}_i^{EB}{\mid }_{\rho =0}$

Remark: the EB(0) predictor (Lyu 2018) is based on ${\hat{\theta }}_0$ obtained from fitting the two parts separately

Simulation Results

$${RDMSE}_i=\frac{MSE\left({\hat{\zeta }}_i^{EB\left(0\right)}\right)-MSE\left({\hat{\zeta }}_i^{EB}\right)}{MSE\left({\hat{\zeta }}_i^{EB}\right)}$$

Simulation Results

$${RBMSE}_i=\frac{{\widehat{MSE}}_i-MSE\left({\hat{\zeta }}_i\right)}{MSE\left({\hat{\zeta }}_i\right)},{CI}_i={\hat{\zeta }}_i\pm 1.96\sqrt{{\widehat{MSE}}_i}$$

Application to CEAP data

Application to CEAP data

Application to CEAP data

Conclusion

• A frequentist approach has been proposed to fit zero-inflated lognormal model with correlated random effects in small area prediction.

• When the true correlation deviates far from 0, the EB(0) predictor (Lyu 2018) has moderately larger MSE than the EB predictor (without model mis-specification).

• The data analysis of the cropland CEAP RUSLE2 measurements collected in South Dakota shows the model assumptions are reasonable.

Appendix B: Simulation results