B Consistency of the MM estimator

On the use of aggregate survey data for estimating regional major depressive disorder prevalence

Supplemental material

Domingo Morales¹, Joscha Krause², and Jan Pablo Burgard²

1Operations Research Center, University Miguel Hern´andez de Elche, Spain

2Department of Economic and Social Statistics, Trier University, Germany

August 14, 2021

A Newton-Raphson algorithm for the MM method

This section derives the components of the Newton-Raphson algorithm for solving system of nonlinear MM equations. The updating formula is

θ^(r+1)=θ^(r)−H⁻¹(θ^(r))f(θ^(r)), whereθ₁=β₁, . . . , θ_p =β_p,θ_p+1 =φand

θ= col

1≤k≤p+1(θ_k), f(θ) = col

1≤k≤p+1(f_k(θ)), H(θ) =

∂f_a(θ)

∂θ_b

a,b=1,...,p+1

.

Recall that the moment generating function of the univariate Gaussiany ∼N(µ, σ²) is ψ(t|µ, σ²) =E

exp{ty}

= Z

R

e^tyf(y|µ, σ²)dy= exp µt+1

2σ²t² ,

and that the moment generating function of the multivariate Gaussiany∼N_p1(µ,Σ) is Ψ(t|µ,Σ) =E

exp{t⁰y}

= Z

R^p1

exp

t⁰y f(y|µ,Σ)dy= exp

t⁰µ+1 2t⁰Σt . We need to calculate the expectations appearing in f(θ) as well as its partial derivatives.

We start with the firstp MM equations that correspond to the regression parameters. The expectation (first moment) ofy_d is calculated according to

E_θ[y_d] = E_u,v

E_θ[y_d|u,v]

=E_u,v[ν_dp_d] =E_u,vh

ν_dexpn

x_dβ+β⁰u_d+φv_doi

= Z

R^p1+1

ν_dexpn

x_dβ+β⁰u_d+φv_do

f(v_d)f(u_d)dv_du_d

= ν_dexp{x_dβ}Z

R

exp{φv_d}f(v_d)dv_dZ

R^p1

exp{β⁰u_d}f(u_d)u_d

= νdexp{x_dβ}ψ(φ|0,1)Ψ(β|0,Σd) =νdexp{x_dβ}exp1

2φ² exp n1

2β⁰₁Σdβ₁ o

= ν_dexpn

x_dβ+ 1 2φ²+1

2β⁰₁Σ_dβ₁o .

Based on these developments, we can conclude that fk(θ) =

D

X

d=1

νdexp n

xdβ+ 1 2φ²+1

2β⁰₁Σdβ₁ o

xdk−

D

X

d=1

ydxdk, k= 1, . . . p.

The derivatives ofE_θ[y_d] are given by

∂E_θ[y_d]

∂βk

= ν_dexpn

x_dβ+1 2

φ²+β⁰₁Σ_dβ₁o

(x_dk+δ⁰_kΣ_dβ₁), k= 1, . . . , p₁,

∂E_θ[y_d]

∂β_k = ν_dexpn

x_dβ+1 2

φ²+β⁰₁Σ_dβ₁o

x_dk, k=p₁+ 1, . . . , p,

∂E_θ[y_d]

∂φ = ν_dexp n

x_dβ+1 2

φ²+β⁰₁Σ_dβ₁ o

φ.

Regarding the second moment, the expectation of y²_d isE_θ[y²_d] =Eu,v

E_θ[y²_d|u,v]

, where Eθ[y²_d|u,v] = varθ[yd|u,v] +E_θ²[yd|u,v] =νdpd+ν_d²p²_d.

Accordingly, it follows that E_θ[y_d²] =E_u,v

E_θ[y_d²|u,v]

= Z

R^p1+1

ν_dp_df(v_d)f(u_d)dv_ddu_d+ Z

R^p1+1

ν_d²p²_df(v_d)f(u_d)dv_ddu_d. We have

I₂ = Z

R^p1+1

p²_df(v_d)f(u_d)dv_ddu_d= Z

R^p1+1

expn

2x_dβ+ 2β⁰u_d+ 2φv_do

f(v_d)f(u_d)dv_ddu_d

= exp

2x_dβ ψ(2φ|0,1)Ψ(2β|0,Σ_d) = exp

2x_dβ exp

2φ² exp n

2β⁰₁Σ_dβ₁ o

= expn

2x_dβ+ 2φ²+ 2β⁰₁Σ_dβ₁o . and

E_θ[y²_d] =ν_dexp

x_dβ+1

2β⁰₁Σ_dβ₁+1

2φ² +ν_d²exp

2x_dβ+ 2β⁰₁Σ_dβ₁+ 2φ² . We can conclude that

fp+1(θ) =

D

X

d=1

νdexp

n

xdβ+1 2φ²+ 1

2β⁰₁Σdβ₁ o

+ν_d²exp n

2xdβ+ 2φ²+ 2β⁰₁Σdβ₁ o

−

D

X

d=1

y²_d.

The derivatives ofE_θ[y_d²] are

∂E_θ[y²_d]

∂β_k = νdexp n

xdβ+1 2φ²+1

2β⁰₁Σdβ₁ o

(xdk+δ⁰_kΣdβ₁) + 2ν_d²exp

n

2xdβ+ 2φ²+ 2β⁰₁Σdβ₁ o

(xdk+δ⁰_kΣdβ₁), k= 1, . . . , p1,

∂E_θ[y²_d]

∂βk

= ν_dexpn

x_dβ+1 2φ²+1

2β⁰₁Σ_dβ₁o x_dk + 2ν_d²exp

n

2xdβ+ 2φ²+ 2β⁰₁Σdβ₁ o

xdk, k=p1+ 1, . . . , p,

∂E_θ[y²_d]

∂φ = ν_dexpn

x_dβ+1 2φ²+1

2β⁰₁Σ_dβ₁o

φ+ 4ν_d²expn

2x_dβ+ 2φ²+ 2β⁰₁Σ_dβ₁o φ.

The elements of the Jacobian matrix are H_rk = ∂fr(θ)

∂θ_k =

D

X

d=1

∂E_θ[y_d]

∂θ_k x_dk, r= 1, . . . , p, k= 1, . . . , p+ 1, H_p+1k = ∂fp+1(θ)

∂θ_k =

D

X

d=1

∂Eθ[y²_d]

∂θ_k , k= 1, . . . , p+ 1.

B Consistency of the MM estimator

This section deals with the consistency of the MM estimator as D → ∞. It presents the adaptations of Lemma 1 and Theorem 1 of Jiang (1998) to the new MEPM model. Jiang (1998) gives a proof for GLMMs where the moments that appear in the MM equations must be calculated by the Monte Carlo method. In the case of the Poisson distribution, these moments can be calculated explicitly and this simplifies the proofs because it is not necessary to impose regularity conditions for the Monte Carlo approach.

Note that for models with intercept, the MM equation 1 is

D

X

d=1

Eθ[yd] =X

d=1

yd.

As

D

X

d=1

y_d 2

=

D

X

d=1

y²_d+

D

X

d6=`

y_dy_`,

we can follow Jiang (1998) and substitute the MM equationp+ 1 by

D

X

d6=`

E_θ[y_dy_`] =

D

X

d6=`

y_dy_`.

LetX ^L

2

−→0 denote L2-convergence, i.e. E[X²]−→0. Denote η0d =ηd(θ0), d= 1, . . . , D, whereθ0 is the true parameter. Define the subset ofR⁴

Q=

(1,0,0,0),(1,1,0,0),(2,1,0,0),(2,2,0,0),(1,1,1,1) .

We first give a lemma that states the convergence of equations (3) and (4) of the main text to zero. In Lemma 1, condition (1) states the requirements for the expectations of the products of derivatives of b(η) =e^η. The conditions (2) and (3) specify the orders of the normalizing constantsa_Dk andb_D. The conclusions appear in (4) and (5).

Lemma 1. Suppose that K_D = max

1≤d,`,d<sup>0</sup>,`⁰≤D max

(a,b,c,r)∈QE_θ0

e^aη0de^{bη0`</sup>e<sup>cη0d0</sup>e<sup>rη0`0}

<∞, lim sup

D→∞

K_D <∞, (1) and that the sequences{a_Dk},k= 1, . . . , p, and {b_D} satisfy

ε_D,k = max

1≤d≤D

1 a²_Dk max

D

X

d=1

x²_dk,

D

X

d=1 D

X

`=1

|x_dkx_d`|

D→∞−→ 0, k= 1, . . . , p. (2)

and

εD,p+1 = max

1≤d≤D

D⁴ b²_D −→

D→∞0. (3)

Then

1 a_Dk

D

X

d=1

x_dk y_d−E_θ0[y_d] L²

D→∞−→ 0, k= 1, . . . p, (4) and

1 bD

D

X

d6=`

ydy`−Eθ0[ydy`] L²

D→∞−→ 0. (5)

Proof. We use the notation E_θ[·] and E[·] for expectations with respect to distributions depending and not depending on θ, respectively. More concretely, we use E_θ[·] for the distributions ofyd andyd|(v_d,ud), andE[·] for the distribution of (v_d,u⁰_d).

For the equationk,k= 1, . . . , p, it holds that Ek=Eθ0

D

X

d=1

xdk yd−E_θ0[yd] 2

=Eθ0

D

X

d=1

xdk yd−b(η˙ 0d) +

D

X

d=1

xdk b(η˙ 0d)−E_θ0[yd] 2

.

We recall that

E_θ0[y_d] =E

E_θ0[y_d|v_d,u_d]

=E e^η0d

=Eb(η˙ _0d) and that

(a+b)² =a²+b²+ 2ab= 2a²+ 2b²−(a²+b²−2ab) = 2a²+ 2b²−(a−b)²≤2a²+ 2b². Take

a=

D

X

d=1

xdk yd−b(η˙ 0d) , b=

D

X

d=1

xdk b(η˙ 0d)−Eθ0[yd] . Then, we have

Ek≤2

Eθ0

D

X

d=1

xdk yd−b(η˙ 0d) 2

+E D

X

d=1

xdk b(η˙ 0d)−Eθ0[yd] 2

= 2(S1+S2).

The first summand is S1=

D

X

d=1

x²_dkE_θ0 h

y_d−b(η˙ _0d)2i +

D

X

d6=`

x_dkx_d`E_θ0 h

y_d−b(η˙ _0d)

y_{`</sub>−b(η˙ <sub>0`})i .

The first expectations are E_θ0h

y_d−b(η˙ _0d)2i

=Eh E_θ0h

y_d−E_θ0[y_d|v_d,u_d]2

v_d,u_dii

=E¨b(η_0d) . Because of the independence, conditioned tov and u, the second expectations are

E_θ0 h

y_d−b(η˙ _0d)

y_{`</sub>−b(η˙ <sub>0`})i

=E

E_θ0 h

y_d−E_θ0[y_d|v_d,u_d]

y_{`</sub>−E<sub>θ0</sub>[y<sub>`}|v_{`</sub>,u<sub>`}] v,u

i

=E

E_θ0h

y_d−E_θ0[y_d|v_d,u_d] v_d,u_di

E_θ0h

y_{`</sub>−E<sub>θ0</sub>[y<sub>`}|v_{`</sub>,u<sub>`}]

v_{`</sub>,u<sub>`}i

= 0.

Therefore, the first summand takes the form

S1 =

D

X

d=1

x²_dkE¨b(η0d) .

The second summand is S₂=

D

X

d=1

x²_dkEh

b(η˙ _0d)−Eb(η˙ _0d)2i +

D

X

d6=`

x_dkx_d`E_θ0h

b(η˙ _0d)−Eb(η˙ _0d) b(η˙ _{0`</sub>)−Eb(η˙ <sub>0`})i .

The expectations are

Eh

b(η˙ _0d)−Eb(η˙ _0d)2i

= var ˙b(η_0d) , E

h b(η˙ _0d)−Eb(η˙ _0d) b(η˙ _{0`</sub>)−Eb(η˙ <sub>0`})i

= cov ˙b(η_0d),b(η˙ _0`) . Therefore, the second summand takes the form

S2 =

D

X

d=1 D

X

`=1

x_dkx_{d`</sub>cov ˙b(η<sub>0d</sub>),b(η˙ <sub>0`}) .

Going back toEk, we have E_k ≤ 2

D

X

d=1

x²_dkE¨b(η_0d) +

D

X

d=1 D

X

`=1

x_dkx_{d`</sub>cov ˙b(η<sub>0d</sub>),b(η˙ <sub>0`})

≤ 2KD

D

X

d=1

x²_dk+

D

X

d=1 D

X

`=1

|x_dkxd`|

.

Thus (4) follows by (1) and (2).

Concerning the equationp+ 1, it holds that Ep+1 =

D

X

d6=`

y_dy_{`</sub>−E<sub>θ0</sub>[y<sub>d</sub>y<sub>`}]

=

D

X

d6=`

(y_d−b(η˙ _0d))(y_{`</sub>−b(η˙ <sub>0`})) +

D

X

d6=`

b(η˙ _0d)(y_{`</sub>−b(η˙ <sub>0`}))

+

D

X

d6=`

(y_d−b(η˙ _0d))˙b(η_0`) +

D

X

d6=`

b(η˙ _0d)˙b(η_{0`</sub>)−E<sub>θ0</sub>[y<sub>d</sub>y<sub>`}] =I₁+I₂+I₃+I₄.

In what follow, we apply the following property. IfX₁, . . . , X_nare independent withE[X_i] = 0 andE[X_i⁴]<∞, and if A= aij

1≤i,j≤n is a symmetric matrix, then E

n

X

i=1 n

X

j=1

a_ijX_iX_j−

n

X

i=1

a_iiE[X_i²] 2

=

n

X

i=1

a²_iivar(X_i²) + 2

n

X

i6=j

a²_ijE[X_i²]E[X_j²]

≤ 2

n

X

i6=j

a²_ijE[X_i²]E[X_j²].

Using this fact and recalling that

E_θ0[y_d|v_d,u_d] = ˙b(η_0d) =e^η0d, var_θ0[y_d|v_d,u_d] = ¨b(η_0d) =e^η0d, leads to

E_θ0[I₁²] = E

E_θ0 D

X

d=1 D

X

`=1

δ_{d`</sub>(y<sub>d</sub>−b(η˙ <sub>0d</sub>))(y<sub>`}−b(η˙ _0`)) 2

v,u

≤ 2E D

X

d6=`

E_θ0 h

(y_d−b(η˙ _0d))² v_d,u_d

i E_θ0

h

(y_{`</sub>−b(η˙ <sub>0`}))² v_d,u_d

i

= 2

D

X

d6=`

Eh

¨b(η_0d)¨b(η_0`)i

= 2

D

X

d6=`

Eh

e^η0de^η0`i

≤2D(D−1)K_D ≤2D²K_D,

whereδ_d` is the Kronecker’s delta.

E_θ0[I₂²] = E

E_θ0 D

X

d=1 D

X

`=1

δ_{d`</sub>b(η˙ <sub>0d</sub>)(y<sub>`}−b(η˙ _0`)) 2

v,u

≤ 2E D

X

d6=`

b˙²(η_0d)E_θ0h

(y_{`</sub>−b(η˙ <sub>0`}))²

v_{`</sub>,u<sub>`}i

= 2

D

X

d6=`

Eh

b˙²(η_0d)¨b(η_0`)i

= 2

D

X

d6=`

Eh

e^2η0de^η0`i

≤2D(D−1)K_D ≤2D²K_D.

SimilarlyE_θ0[I₃²]≤2D²K_D.

E_θ0[I₄²] = E D

X

d6=`

b(η˙ _0d)˙b(η_{0`</sub>)−E<sub>θ0</sub>[y<sub>d</sub>y<sub>`}] 2

=

D

X

d6=`

D

X

d⁰6=`⁰

E n

b(η˙ 0d)˙b(η0`)−Eb(η˙ 0d)˙b(η0`)on

b(η˙ _0d⁰)˙b(η_{0`</sub><sup>0</sup>)−Eb(η˙ <sub>0d</sub><sup>0</sup>)˙b(η<sub>0`}⁰)o

=

D

X

d6=`

D

X

d⁰6=`⁰

cov

b(η˙ _0d)˙b(η_{0`</sub>),b(η˙ <sub>0d</sub><sup>0</sup>)˙b(η<sub>0`}⁰)

≤

D

X

d6=`

D

X

d⁰6=`⁰

Eh

b(η˙ _0d)˙b(η_{0`</sub>)˙b(η<sub>0d</sub><sup>0</sup>)˙b(η<sub>0`}⁰)i

=

D

X

d6=`

D

X

d⁰6=`⁰

Eh

e^η0de^{η0`</sup>e<sup>η0d0</sup>e<sup>η0`0}i

≤D⁴K_D.

Thus (5) follows by (1) and (3).

Let us define

M_D,k(θ) = 1 a_Dk

D

X

d=1

x_dkE_θ[y_d], M_D,p+1(θ) = 1 b_D

D

X

d=1

E_θ[y_d²],

Mˆ_D,k = 1 a_Dk

D

X

d=1

x_dky_d, Mˆ_D,p+1= 1 b_D

D

X

d=1

y²_d,

M˜_D,k(θ) = 1 aDk

D

X

d=1

x_dkb(η˙ _d), M˜_D,p+1(θ) = 1 bD

D

X

d6=`

b(η˙ _d)˙b(η_`), and

MD(θ) = MD,k(θ)

1≤k≤p+1,Mˆ D = MˆD,k

1≤k≤p+1,M˜ D(θ) = M˜D,k(θ)

1≤k≤p+1. Lemma 1 gives sufficient conditions for

kMˆ _D−M_D(θ₀)k −→^L2

D→∞0

Let {A_Dk}, 1 ≤ k ≤ p, and {B_D} any sequences such that A_Dk → ∞ and B_D → ∞ as D→ ∞. Let ˆθ = ( ˆβ⁰,φ) be anyˆ θ = (β⁰, φ)∈Θ_D =

θ : |β_k| ≤ A_Dk,1≤k≤p; 0< φ≤ BD satisfying the inequality

M˜ D(θ)−Mˆ D

≤δD, (6) whereδD → ∞ asD→ ∞.

The following theorem states the consistency of the MM estimators of the parameters of the MEPM model. The proof follows from Theorem 1 of Jiang (1998) by doing the particularization to the current MEPM model and by noting that:

(1) Jiang (1998) gives a proof for GLMMs with q random effects. Here we consider the Poisson mixed model with one random effect, but we addp random measurement errors.

The random errors have a known multivariate normal distribution and are independent of the random effect.

(2) The moments M_k(θ) are calculated analytically and not approximated by the Monte Carlo method, and

(3) The expectations E[·] are taken with respect to the joint distribution of vd and ud, d= 1, . . . , D, as it was shown in Lemma 1.

Because of (1)-(3), there are not remarkable difficulties in adapting and particularizing the proof of Jiang. We present the sketch of the proof. For more details, see Jiang (1998).

Theorem 1. Suppose that the conditions of Lemma 1 are satisfied. Let θ₀ be the true parameter and letεD = max1≤k≤p+1{ε_D,k}.

(a) If ε_D/δ_D² → ∞, then θˆ exists with probability tending to one as D→ ∞.

(b) If, furthermore, all the first derivatives of the expectations Eθ[yd] = E[˙b(ηd(θ))] and E_θ[y_dy_{`</sub>] = E[˙b(η<sub>d</sub>(θ))˙b(η<sub>`}(θ))], with respect to the components of θ, can be taken under the expectation sign; the quantities

sup

kθk≤B

Eh

b(η˙ _d(θ))4i , E

sup

kθk≤B

¨b(η_d(θ))

, E

sup

kθk≤B

¨b(η_d(θ))

b(η˙ _`(θ))

, d, `= 1, . . . , D,

are finite for anyB >0; and lim inf

D→∞ inf

kθ⁻θ0k>ε

k|M_D(θ)−MD(θ0)k>0 for any ε >0, then ˆθ converges to θ0 with probability tending to one.

Proof.

(a) By the proof of Lemma 1, we have E_θ0

Mˆ _D−M_D(θ₀)

2

=

p+1

X

k=1

E_k≤(2p+ 4)K_Dε_D −→

D→∞0.

On the other hand Eθ0

( ˜MD,k(θ0)−M_D,k(θ0))²

=E 1

a_Dk

D

X

d=1

xdk b(η˙ 0d)−E[˙b(η0d)]

2

≤ KD

a²_Dk

D

X

d6=`

|x_dkxd`| ≤KDεD.

and similarly

E_θ0

( ˜M_D,p+1(θ₀)−M_D,p+1(θ₀))²

≤ K_D

b²_D D⁴ ≤K_Dε_D. Thus, we have

P

M˜ D(θ0)−Mˆ D

≤δD

≤ P

M˜ D(θ0)−MD(θ0) ≤ δD

2

+P

MD(θ0)−Mˆ D

≤ δD

2

≤ [(p+ 1) + (2p+ 4)]K_Dε_D δ_D −→

D→∞0.

Therefore, (6) holds and ˆθ exists with probability tending to one at θ=θ₀.

(b) Because of the continuity ofM with respect toθ, it enough to prove that for anyδ >0, it holds

P kM_D(ˆθ)−MD(θ0)k ≥δ

D→∞−→ 0.

We have that

kM_D(ˆθ)−M_D(θ0)k ≤ kM_D(ˆθ)−M˜ _D(ˆθ)k+kM˜ _D(ˆθ)−Mˆ _Dk+kMˆ _D−M_D(θ0)k

≤ sup θ^∈ΘD

kM_D(θ)−M˜ D(θ)k+δD +kMˆ D−MD(θ0)k.

Because of Lemma 1, the third summand becomes bounded and close to zero as D→ ∞.

By expanding M_D(ˆθ) and ˜M_D(ˆθ) in Taylor’s series around some θ∗ in a neighborhood of θ₀, we mimic the same steps of the proof of Theorem 1(b) of Jiang (1988) to obtain a bound for the first summand as D → ∞. This fact allows applying the the Chebyshev’s inequality for bounded random variables and we get

P kM_D(ˆθ)−M_D(θ₀)k ≥δ

≤ 1 δE

kM_D(ˆθ)−M_D(θ₀)k

−→

D→∞0.

This completes the proof.

C Additional simulation experiments

C.1 Set up

This section presents additional simulation experiments, in which we compare the explicit measurement error modelling implemented by the MEPM model against the simulative approach used by SIMEX. For the latter, we use the R-packagesimex provided by Lederer et al. (2019). Please note that this package does not yet allow to fit generalized linear mixed models with measurement errors. However, it allows for the fitting of generalized linear models with measurement errors (without random effects). Thus, in order to avoid giving SIMEX an unfair disadvantage, we alter the data generating process presented in Section 5 of the main text by setting the random effect standard deviation to zero. Accordingly, we consider the measurement error Poisson (MEP), the standard Poisson (P) model, and the SIMEX approach that uses the standard Poisson as naive model. Since the objective of the paper is the prediction of regional prevalence figures, we focus on mean parameter prediction in the simulation hereafter.

A Monte Carlo simulation withI = 750, i= 1, . . . , I, iterations is conducted. We generate a population ofD domains, whereDvaries over scenarios. Ford= 1, ..., D, we define

µ_d=ν_dp_d, y_d∼Poiss(µ_d), p_d= exp{β₀+x_1,dβ₁+x_2,dβ₂+u⁰_1,dβ₁+u⁰_2,dβ₂}, where ν_d = 300, β₀ = −4 β₁ = (0.5,0.5)⁰, and β₂ = −β₁. Accordingly, we have an intercept and four covariates that are measured with error. Note that the random effects are generated in every iteration individually. The unbiased covariate predictors are drawn from uniform distributions according tox_dj ∼U(1.0,1.4), j = 1, ...,4, and held fixed over all Monte Carlo iterations. The covariate measurement errorsud= (u⁰_1,d,u⁰_2,d)⁰ are drawn in each iteration individually according to

u_d∼N4(0,Σ_d), Σ_d=











σ²_1,d σ_12,d σ_13,d σ_14,d σ_21,d σ²_2,d σ_23,d σ_24,d σ_31,d σ_32,d σ_3,d² σ_34,d σ_41,d σ_42,d σ_34,d σ_4,d²











, d= 1, . . . , D,

whereσ²_j,d ∼U(0.05,0.15), σ_jk,d = ρ_jkσ²_j,dσ_k,d² , ρ_jk = 0.5 forj = 1 andk = 2,3,4, as well asρ_jk =−0.3 forj = 2,3 and k= 3,4, j 6=k. Just like in Section 5 of the main text, we consider four simulation scenarios arising from the four different values forD. The scenarios are defined as in Table 1 of the main text.

C.2 Results

We consider relative mean squared error (RMSE), relative root mean squared error (RRMSE), absolute bias (ABIAS), as well as relative absolute bias (RABIAS) as performance measures.

They are given as follows:

RM SE_d= 1 I

I

X

i=1

(µ⁽ⁱ⁾_d −µˆ⁽ⁱ⁾_d )²

!1/2

, RRM SE_d= RM SE_d

¯ µd

, µ¯_d= 1 I

I

X

i=1

µ⁽ⁱ⁾_d ,

ABIAS_d= 1 I

I

X

i=1

|µ⁽ⁱ⁾_d −µˆ⁽ⁱ⁾_d |, RABIAS_d= ABIAS_d

¯ µd

.

We further define the subsequent aggregated measures RM SE = 1

D

D

X

d=1

RM SE_d, ABIAS = 1 D

D

X

d=1

ABIAS_d,

RRM SE= 1 D

R

X

d=1

RRM SEd·100, RABIAS = 1 D

D

X

d=1

RABIASd·100%, to allow for compact presentation. The results are summarized in Table 1 and Figure 1.

Method Scenario RM SE RRM SE ABIAS RABIAS

P 1 0.7086 19.2049 0.5631 15.4376

SIMEX 1 0.6963 18.8702 0.5548 15.2096

MEP 1 0.5770 15.6044 0.4460 12.2295

P 2 0.6376 17.1027 0.5059 13.5730

SIMEX 2 0.6226 16.7053 0.4958 13.4092

MEP 2 0.5357 14.3590 0.4163 11.1577

P 3 0.5747 15.4733 0.4560 12.2784

SIMEX 3 0.5615 15.1224 0.4471 12.1412

MEP 3 0.5477 14.7210 0.4236 11.3904

P 4 0.5449 14.7829 0.4318 11.7164

SIMEX 4 0.5294 14.3652 0.4212 11.5733

MEP 4 0.4788 12.9803 0.3746 10.1575

Table 1: Results of mean parameter prediction

141822

Scenario 1

RRMSE 121620

Scenario 2

RRMSE

101418

Scenario 3

RRMSE 1216

Scenario 4

RRMSE

P SIMEX MEP

P SIMEX MEP

P SIMEX MEP

P SIMEX MEP Figure 1: Domain-level mean parameter prediction performance

In Table 1, we see that the SIMEX approach outperforms the standard Poisson model in all scenarios and for all performance measures. This could be expected from theory, since SIMEX has been introduced as an effective and flexible technique to deal with measure- ment errors. However, the SIMEX approach does not reach the efficiency level of the MEP approach. The latter dominates in all scenarios and all performance measures by a consid- erable margin. This is likely due to the explicit consideration of the measurement errors in the model equation. All inferential steps are strictly derived under this premise. That is to say, for this particular setting, the proposed MEP approach is the best considered option.

Figure 1 displays the measureRRM SE_d for all considered domains and simulation scenar- ios. The predictions obtained under the standard Poisson model are marked in blue. The predictions from the SIMEX approach are plotted in light red. And finally, the predictions obtained under the proposed MEP approach are displayed in red. In accordance with Table 1, we see that the MEP approach outperforms the other methods significantly. Its RRMSE on domain-level is always smaller. However, we further see that the range of RRMSE values is also much narrower for the MEP approach. This implies that the method overall obtains more stable results relative to standard Poisson and SIMEX.

References

Jiang, J. (1998). Consistent estimators in generalized linear models. Journal of the Amer- ican Statistical Association,93, 720-729.

Lederer, W., Seibold, H., K¨uchenhoff, H., Lawrence, C., Brondum, R. F. (2019). Package

‘simex’. URL: https://cran.r-project.org/web/packages/simex/simex.pdf.