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,i=1,..,D;j=1,..,Ni$y_{ij}^{\ast },i=1,..,D;j=1,..,N_i$ satisfies y∗ij=yijδij$$y_{ij}^{\ast }=y_{ij}{\delta }_{ij}$$

• Positive part: log(yij)=β0+z′1,ijβ1+ui+eij where eiji.i.d∼N(0,σ2e).

• Binary part: Pr(δij=1)=pij, g(pij)=α0+z′2,ijα1+bi where g(⋅) 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

  • (uibi)iid∼BVN(0,Σub),Σub=(σ2uρσuσbρσuσbσ2b)

  • "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 ζi=1Ni∑Nij=1y∗ij is ˆζMMSEi=1Ni{∑j∈siy∗ij+∑j∈ˉsiˆy∗MMSEij} where ˆy∗MMSEij=Eθ(y∗ij|y∗si)

Let ˆθ be a consistent estimator of model parameter θ. ˆζEBi=1Ni{∑j∈siy∗ij+∑j∈ˉsiˆy∗EBij} where ˆy∗EBij=ˆy∗MMSEij(ˆθ)

Conditional Distribution

Result 1:

ui|bi,y∗si∼N(˜μui,˜σ2ui) where ˜μui=γiˉ˜ri+(1−γi)[ρσuσbbi],˜σ2ui=γiσ2e/˜ni and γi=(1−ρ2)σ2u[(1−ρ2)σ2u+σ2e/˜ni]−1,˜ni=∑j∈siδij, ˉ˜ri=˜n−1i∑j∈siδij{log(yij)−β0−z′1,ijβ1}

Remark: when ρ=0, γi=σ2u/(σ2u+σ2e/˜ni), ui|y∗si∼N(γiˉ˜ri,γiσ2e/˜ni).

Conditional Distribution (Cont'd)

Result 2:

f(bi|y∗si)∝πsi(bi)ϕ(bi−mi√vi) where ϕ(⋅) is the probability density function of standard normal distribution, πsi(bi)=∏j∈si[pijδij+(1−pij)(1−δij)] and pij(bi)=g−1(α0+α1zij+bi). (mi,vi)=(ρσuσbσ2u+σ2e/˜niˉ˜ri,1−ρ21−(1−γi)ρ2σ2b). Remark: when ρ=0, (mi,vi)=(0,σ2b).

Conditional Distribution (Cont'd)

Result 3:

E(y∗ij|y∗si,θ)=hij(θ)E(pij(bi)η(bi)|y∗si) where hij(θ)=exp(β0+x′ijβ1+γiˉ˜ri+˜σ2ui/2+σ2e/2) and η(bi)=exp((1−γi)ρσu/σbbi).

👉 ˆy∗EBij=E(y∗ij|y∗si,ˆθ)

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

Parameter Estimation

Given ρ≠0, how to get a consistent estimator for model parameter θ=(β′,α′,σ2u,σ2b,σ2e,ρ)′? 😣

Full likelihood function: Li(θ)=(1−γi)1/2(vi/σ2b)1/2(2πσ2e)˜ni/2exp(γiˉ˜r2i2σ2e/˜ni+m2i2vi−∑j˜r2ij2σ2e)∫πsi(bi)1√viϕ(bi−mi√vi)dbi Remark:

• A good starting value and the gradient vector ∂L/∂θ help optim find the optimizer much faster. ✌️

• Profiling out ρ takes much longer time to find the MLE.

MSE Estimator

Formally, the MSE of an EB predictor is MSE(ˆζEBi)=M1i(θ)+M2i(θ) where M1i(θ)=E(ˆζMMSEi−ζi)2=O(1), as D→∞ and M2i(θ)=E(ˆζEBi−ˆζMMSEi)2=o(1), as D→∞ Recall that ˆζMMSEi=E[ζi|y∗si,θ],ˆζEBi=E[ζi|y∗si,ˆθ]

MSE Estimator (Cont'd)

So the MSE of the optimal predictor is M1i(θ)=E(E[ζi|y∗si,θ]−ζi)2=E[V(ζi|y∗si,θ)] We propose a One-Step MSE estimator defined by V(ζi|y∗si,ˆθ)=1N2i∑j∈ˉsi∑k∈ˉsi[E{y∗ijy∗ik|y∗si,ˆθ}−E{y∗ij|y∗si,ˆθ}E{y∗ik|y∗si,ˆθ}] where E{y∗ijy∗ik|y∗si}=hijhikexp(˜σ2ui)E[pijpikη(2bi)|y∗si],j≠k and E{y∗ijy∗ik|y∗si}=h2ijexp(˜σ2ui+σ2e)E[pijη(2bi)|y∗si],j=k

Alternate MSE Estimators

• Denote original parameter estimate ˆθ, bootstrap population Y∗(b), bootstrap sample y∗(b), parameter estimate ˆθ(b)

• For b=1,…,B, we obtain Y∗(b)⟶ˉy∗(b)Ni y∗(b),ˆθ⟶ˆˉy∗(b)MMSENi y∗(b),ˆθ(b)⟶ˆˉy∗(b)EBNi

1. ^MSEBooti=B−1∑Bb=1(ˆˉy∗(b)EBNi−ˉy∗(b)Ni)2

2. ^MSESemi-Booti=ˆM1i+ˆMBoot2i where a bootstrap estimator of M2i is defined by ˆMBoot2i=B−1∑Bb=1(ˆˉy∗(b)EBNi−ˆˉy∗(b)MMSENi)2

Simulation Setting

• Number of areas: D=60

• Sample rate: (Ni,ni)=(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: zij∼N(4.45,0.055)

• β=(−13,2)′,α=(−20,5)′, (σ2u,σ2e,σ2b)=(0.22,1.23,0.52)

• ρ∈{−0.9,−0.6,−0.3,0,0.3,0.6,0.9}

• Compare the proposed EB predictor with the EB(0) predictor ˆζEBi∣ρ=0

Remark: the EB(0) predictor (Lyu 2018) is based on ˆθ0 obtained from fitting the two parts separately

Simulation Results

RDMSEi=MSE(ˆζEB(0)i)−MSE(ˆζEBi)MSE(ˆζEBi)

Simulation Results

RBMSEi=^MSEi−MSE(ˆζi)MSE(ˆζi),CIi=ˆζi±1.96√^MSEi

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