Information Principles in Random Effects Models

Information Principles in Random Effects Models

Herwig Friedl G¨ oran Kauermann August 14, 2013

1 Information Matrices

In this note, data (y, z) are considered, where y denotes the observable part and z refers to that part which is unobservable. Later we will concentrate on a response variabley which is modelled in terms of a fixed predictor variable and an additional random effectz.

1.1 Complete Likelihood

Letθdenote all unknown parameters in the model and λ_c(y, z;θ) = logf(y, z;θ)

be the complete log-likelihood function corresponding to the joint distribution of the response and the random effect. We consider the model

f(y, z;θ) =f(y|z;θ)×f(z;θ)

wheref(y|z;θ) is the conditional model, given the random effectz, andf(z;θ) denotes the random effect density. Both,f(y|z;θ) andf(z;θ) possibly depend on unknown parametersθ.

The respective complete score vector is

S_c(y, z;θ) = ∂

∂θλ_c(y, z;θ)

= f^′(y, z;θ)

f(y, z;θ) (1)

with complete negative second derivative Ic(y, z;θ) = − ∂²

∂θ∂θ^tλc(y, z;θ)

= −f^′′(y, z;θ)

f(y, z;θ) +f^′(y, z;θ) f(y, z;θ)

f^′t(y, z;θ) f(y, z;θ)

= −f^′′(y, z;θ)

f(y, z;θ) +S_c(y, z;θ)S_c^t(y, z;θ). (2)

1.2 Observed Likelihood

The MLE ˆθis constructed by maximizing the observed log-likelihood λ(y;θ) = log

∫

f(y, z;θ)dz

which depends on the observationsyonly. In what follows we need the exchangeability of integra- tion and differentiation of the complete density function. From the respective observed scores

S(y;θ) = ∂

∂θλ(y;θ)

=

∫ f^′(y, z;θ)dz

∫f(y, z;θ)dz

=

∫ f^′(y, z;θ) f(y, z;θ)

f(y, z;θ)

∫ f(y, z;θ)dzdz

= E_θ(S_c(y, z;θ)|y) (3)

we get the matrix of observed negative second derivatives I(y;θ) = − ∂²

∂θ∂θ^tλ(y;θ)

= −

∫ f^′′(y, z;θ)dz

∫ f(y, z;θ)dz +

∫ f^′(y, z;θ)dz

∫ f(y, z;θ)dz

∫f^′t(y, z;θ)dz

∫f(y, z;θ)dz

= −

∫ f^′′(y, z;θ) f(y, z;θ)

f(y, z;θ)

∫ f(y, z;θ)dzdz+S(y;θ)S^t(y;θ)

= Eθ

(

−f^′′(y, z;θ) f(y, z;θ)

y )

+S(y;θ)S^t(y;θ).

Because−f^′′(y, z;θ)/f(y, z;θ) =Ic(y, z;θ)−Sc(y, z;θ)S_c^t(y, z;θ) from (2), this can be rewritten as I(y;θ) = Eθ(Ic(y, z;θ)|y)−Eθ

(Sc(y, z;θ)S_c^t(y, z;θ)|y)

+Eθ(Sc(y, z;θ)|y)Eθ(S^t_c(y, z;θ)|y)

= E_θ(I_c(y, z;θ)|y)−V ar_θ(S_c(y, z;θ)|y) (4) often called the missing information principle. This result is due to Louis (1982) and used to extract the observed information matrix when the EM algorithm is applied to find MLEs in incomplete data problems. Moreover, it provides a means of estimating the information which is associated with the MLEs and requires only the computation of a complete-data gradient vector and a second derivative matrix but not those associated with the incomplete-data likelihood.

1.3 Missing Information Principle

Denote the observed and complete Fisher information matrices by I(θ) = Eθ(I(y;θ)) =Eθ(S(y;θ)S^t(y;θ))

Ic(θ) = E_θ(I_c(y, z;θ)) =E_θ(S_c(y, z;θ)S_c^t(y, z;θ)),

where expectations are taken over the marginal and the joint density, respectively. From λ(y;θ) =λc(y, z;θ)−logf(z|y;θ)

we get the identity

− ∂²

∂θ∂θ^tλ(y;θ) = − ∂²

∂θ∂θ^tλ_c(y, z;θ) + ∂²

∂θ∂θ^tlogf(z|y;θ) I(y;θ) = Ic(y, z;θ) + ∂²

∂θ∂θ^tlogf(z|y;θ).

By taking the expectation with respect to the conditional density ofz, given the datay, we have I(y;θ) =E_θ(I_c(y, z;θ)|y)−E_θ

(

− ∂²

∂θ∂θ^tlogf(z|y;θ)y )

.

McLachlan and Krishnan (1997) define

Ic(y;θ) := Eθ(Ic(y, z;θ)|y) Im(y;θ) := Eθ

(

− ∂²

∂θ∂θ^tlogf(z|y;θ)y )

to rewrite the above identity as a difference of conditionally expected informations, namely I(y;θ) =Ic(y;θ)− Im(y;θ), (5) whereIm(y;θ) is considered as the missing information as a consequence of observingy only and not alsoz. As stated in Louis (1982) the comparison of (5) with (4) results in an easy to interpret expression for the missing information, i.e.

Im(y;θ) =V arθ(Sc(y, z;θ)|y). (6) Note that with respect to the marginal density we have

E_θ(Ic(y;θ)) =

∫ (∫

I_c(y, z;θ)f(z|y;θ)dz )

f(y;θ)dy=

∫ ∫

I_c(y, z;θ)f(y, z;θ)dy dz=Ic(θ).

By taking marginal expectations in (5) and by using the above result, we get

I(θ) =Ic(θ)−E_θ(Im(y;θ)) (7)

for the a priori expected information.

1.4 EM Estimates

Dempster, Laird and Rubin (1977) iteratively maximizes Q(θ|θ^(t)) =

∫

λc(y, z;θ)f(z|y;θ^(t))dz=E_θ(t)(λc(y, z;θ)|y),

the conditional expected joint log-likelihood given the data. They also showed that ifθ^(t)converges to a point ˆθ then the observed score function has a zero at θ^(∞) = ˆθ, the marginal MLE. Each iteration of the EM algorithm solves

∂

∂θQ(θ|θ^(t)) θ=θ^(t+1)

= 0. (8)

Louis (1982) and later on Meilijson (1989) consideredQand its derivatives when evaluated in θ=θ^(t), e.g.

Q(θ|θ^(t))

θ=θ^(t)

= E_θ(t)(λ_c(y, z;θ^(t))|y)

∂

∂θQ(θ|θ^(t)) θ=θ^(t)

= E_θ(t)(S_c(y, z;θ^(t))|y) =S(y;θ^(t))

− ∂²

∂θ∂θ^tQ(θ|θ^(t)) θ=θ^(t)

= E_θ(t)(I_c(y, z;θ^(t))|y) =Ic(y;θ^(t)).

Evaluating the missing information (6) at ˆθ gives us the simplification

Im(y; ˆθ) =E_θ(S_c(y, z;θ)S_c^t(y, z;θ)|y)|_{θ= ˆθ} (9) becauseEθ(Sc(y, z;θ)|y)

θ= ˆθ=S(y; ˆθ) = 0 at convergence, a consequence of (8).

A software which provides EM estimates by directly applying successive E and M steps will automatically provide the inverse of the matrixIc(y; ˆθ), an estimate of the complete information at convergence. But this matrix should not be used for estimating the standard errors of the MLE θ, because it does not account for the missing information. Instead of that, Efron and Hinkleyˆ (1978) suggest to use the inverse of the observed informationI(y; ˆθ) to serve as an estimate of the covariance matrix of the MLE.

If we take the observed score vector S(y;θ) and calculate its negative derivative then we get the desired observed information

−∂

∂θS(y;θ) =I(y;θ).

The Newton-Raphson procedure

θ^(t+1)=θ^(t)+I(y;θ^(t))⁻¹S(y;θ^(t))

will giveI(y; ˆθ)⁻¹at convergence. One remaining open problem is concerned with the question on how we can compute or approximate all the expected values that we need.

2 A Model to Handle Overdispersion

Now we consider a Generalized Linear Model for the conditional meanµ=µ(z) =E(y|z;β) of an observationy given an unobservable random effect zof the form

g(µ) =x^tβ+z

where xincludes all explanatory variables, y is a response variable and z is the random effect.

β is the vector of unknown parameters associated to the fixed effects and we like to construct estimates for the standard error of the maximum likelihood estimate (MLE) ˆβ. We also assume that distributional assumptions can be made on the response variable conditional on the random effect.

2.1 Gaussian Random Effects

First let us assume that datay= (y1, . . . , yn)^tare available that conditionally follow a Generalized Linear Model where the random effectsz = (z1, . . . , zn)^t are independently drawn from a normal distribution with zero mean and varianceσ_z², e.g.

g(µ_i) =x^t_iβ+σ_zz_i,

whereg is a (canonical) link function andz_i^iid∼ N(0,1). Here,θ= (β^t, σ_z)^tand the density of the random effectsf(z) =φ(z) does not depend on any unknown parameter. Let ˜X = (X|z) denote the design matrixX extended by the vector of random effectszand write

g(µ) = ˜Xθ.

For f(y|z;θ) from the exponential family and a canonical link model we have the well known results

S_c(y, z;θ) = f^′(y|z;θ)×φ(z)

f(y|z;θ)×φ(z) =f^′(y|z;θ) f(y|z;θ) = ∂

∂θlogf(y|z;θ)

= X˜^t(y−µ) Ic(y, z;θ) = − ∂²

∂θ∂θ^tlogf(y|z;θ)

= X˜^tVX.˜

Here,V =V(µ(z)) denotes the conditional variance function, i.e. the variance function tof(y|z;θ) as in the usual GLM setting fory|z. Therefore, the information matrices in (5) are

Ic(y;θ) = E_θ (

I_c(y, z;θ)y )

=E_θ

(X˜^tVX˜y ) Im(y;θ) = Eθ

(

Sc(y, z;θ)S_c^t(y, z;θ)y )−Eθ

(

Sc(y, z;θ)y )

Eθ

(

S_c^t(y, z;θ)y )

= Eθ

(X˜^t(y−µ)(y−µ)^tX˜y )−Eθ

(X˜^t(y−µ)y )

Eθ

(

(y−µ)^tX˜y )

.

We approximate the above conditional expectations through Gauss-Hermite quadrature. That means that the unobservable effects zi are replaced by some known masspoints ζk with known massesπk,k= 1, . . . , K. Further, let ˜xik= (x^t_i, ζk)^tand

wik= f(yi|ζk;θ)πk

∑K

l=1f(yi|ζl;θ)πl

the respective approximation tof(zi|yi;θ). This gives an approximation to thecomplete Fisher information

Ic(y;θ) = ∑

i

Eθ(˜xix˜^t_iVi|yi)

≈ ∑

i

∑

k

˜

x_ikx˜^t_ikV_ikw_ik,

which is automatically provided at convergence of the EM algorithm. For the first term inIm(y;θ) we get

Eθ



∑

i

˜

xix˜^t_i(yi−µi)²+∑

i̸=j

˜

xix˜^t_j(yi−µi)(yj−µj)y





≈∑

i

∑

k

˜

xikx˜^t_ik(yi−µik)²wik+∑

i̸=j

∑

k

∑

l

˜

xik˜x^t_jl(yi−µik)(yj−µjl)wikwjl. The second term ofIm(y;θ) is zero atθ= ˆθ. Generally, it can be approximated by E_θ

(∑

i

˜

x_i(y_i−µ_i)y )

E_θ



∑

j

˜

x^t_j(y_j−µ_j)y



≈∑

i

∑

k

˜

x_ik(y_i−µ_ik)w_ik∑

j

∑

l

˜

x^t_jl(y_j−µ_jl)w_jl. Subtracting the last from the previous result gives the approximation to themissing information

Im(y;θ) = ∑

i

Eθ

(x˜ix˜^t_i(yi−µi)²|yi

)−∑

i

Eθ(˜xi(yi−µi)|yi)Eθ

(x˜^t_i(yi−µi)|yi

)

≈ ∑

i

∑

k

˜

xikx˜^t_ik(yi−µik)²wik−∑

i

∑

k

∑

l

˜

xik˜x^t_il(yi−µik)(yi−µil)wikwil. Hence, theobserved informationis

I(y;θ) = ∑

i

E_θ(˜x_ix^t_iV_i|y_i)

−∑

i

E_θ(

˜

x_ix˜^t_i(y_i−µ_i)²|y_i)

+∑

i

E_θ(˜x_i(y_i−µ_i)|y_i)E_θ(

˜

x^t_i(y_i−µ_i)|y_i)

≈ ∑

i

∑

k

˜

xikx˜^t_ikVikwik

−∑

i

∑

k

˜

x_ikx˜^t_ik(y_i−µ_ik)²w_ik+∑

i

∑

k

∑

l

˜

x_ikx˜^t_il(y_i−µ_ik)(y_i−µ_il)w_ikw_il.

It is also worth to consider the a priori expected informationin (7) I(θ) =Eθ(I(y;θ)) =Ic(θ)−Eθ(Im(y;θ)).

Because of the assumed canonical link model, the matrix Ic(y;θ) =Eθ(Ic(y, z;θ)|y) = ˜X^tVX˜ is not a function of the observed data. Therefore,

Ic(θ) =Ic(y;θ) =∑

i

E_θ(˜x_ix˜^t_iV_i|y_i) =∑

i

E_θ(˜x_ix˜^t_iV_i).

Moreover, this also holds for Eθ

(∑

i

Eθ

(x˜ix˜^t_i(yi−µi)²|yi

))

= ∑

i

∫ ∫

˜

xix˜^t_i(yi−µi)²f(yi|zi;θ)f(zi)

f(y_i;θ) f(yi;θ)dy dz

= ∑

i

∫

˜ xix˜^t_i

(∫

(yi−µi)²f(yi|zi;θ)dy )

f(zi)dz

= ∑

i

∫

˜

x_ix˜^t_iV_if(z_i)dz=∑

i

E_θ(˜x_ix˜^t_iV_i).

Together this gives I(θ) = Eθ

(∑

i

Eθ(˜xi(yi−µi)|yi)Eθ(˜x^t_i(yi−µi)|yi) )

=∑

i

Eθ(s(yi;θ)s^t(yi;θ))

= ∑

i

V ar_θ(s(y_i;θ)) (10)

the variance-covariance matrix of the observed total score. For i.i.d. variates Meilijson (1989) uses such an empirical Fisher information. Here we like to suggest to estimate the variance contribu- tion of each individual score in an similar empirical manner for non-i.i.d. scores. Estimating the individual variance by ’samples of size one’ at a time, e.g.V ar(s(yd i;θ)) =s(yi; ˆθ)s^t(yi; ˆθ), results in ˆI =∑

is(yi; ˆθ)s^t(yi; ˆθ) with respective approximation I ≈ˆ ∑

i

∑

k

∑

l

˜

xikx˜^t_il(yi−µˆik)(yi−µˆil) ˆwikwˆil.

2.2 Unspecified Random Effect Distribution

Let us again consider the modelg(µ_i) =x^t_iβ+z_i. Like in the previous part,µ_i=E(y_i|z_i) denotes the conditional mean but z_i are now i.i.d. from an unknown distribution, which should be estimated nonparametrically. This can be done by an estimate, that is defined by giving masses π= (π₁, . . . , π_K)^ton a finite numberK of masspointsζ= (ζ₁, . . . , ζ_K)^t; both vectors are treated as unknown, whereas the numberK itself is assumed to bea priori known. In the following the explanatory vectorxshould not include an intercept term! Hence, we rewrite the model as

g(µi) =x^t_iβ+ϵ^t_iζ= ˜x^t_iθ

where ˜xi= (x^t_i, ϵ^t_i)^t,θ= (β^t, ζ^t)^t andϵi= (ϵi1, . . . , ϵiK)^{t iid}∼ M N(1;π), i.e. multinomial with P(z_i=ζ_k;π) =P(ϵ_i =e_i;π) =

∏K k=1

π_k^eik, with

∑K k=1

π_k= 1, π_k>0,

whereei = (ei1, . . . , eiK)^tis any admissible realization ofϵi. Hence, the unknown parameters can be separated into two distinct sets of parameters, whereθ corresponds tof(y|z) (again a member of the exponential family) and πonly belongs to the discrete estimate off(z) (the random effect density).

Let ˜X = (x^t_i, e^t_i)^t denote the matrix built up by all rows of the design matrixX extended by the rows of the respective row vectorse_i. The completeθscore and its negative derivative under this model is as before

S_c^θ(y, z;θ) = ∂

∂θlog (

f(y|z;θ)f(z;π) )

= f^′(y|z;θ) f(y|z;θ) = ∂

∂θlogf(y|z;θ) = ˜X^t(y−µ), I_c^θ(y, z;θ) = − ∂²

∂θ∂θ^tlogf(y|z;θ) = ˜X^tVX.˜

To determine theπscore under the multinomial model subject to the above constraint, we note that logf(z;π) is, aside from constants,

∑n i=1

(_K−1

∑

k=1

eiklogπk+eiKlog(1−π1−. . .−πK−1) )

.

Therefore, the completeπ-score vector is S_c^π(y, z;π) = ∂

∂πlog (

f(y|z;θ)f(z;π) )

= f^′(z;π) f(z;π) = ∂

∂πlogf(z;π), and equals its marginal analogue. Hence, fork= 1, . . . , K−1 we getS_c^πk(y, z;π) =∑n

i=1

(eik

π_k −^e_π^iK_K) giving the vector

S_c^π(y, z;π) =

∑n i=1

K∑−1 k=1

e_k (eik

π_k −eiK

π_K )

,

where ek = (ek1, . . . , ek,K−1)^t is a K −1 indicator vector with ekk = 1 and zeros otherwise.

Equating the complete score to zero results in the ML estimates ˆ

πc= 1 n

∑n i=1

K∑−1 k=1

ekeik.

Its respective negative derivative is therefore the (K−1)×(K−1) matrix I_c^π(y, z;π) =

(

− ∂

∂πj

S_c^πk(y, z;π) )

k,j

= { ∑n

i=1

(e_ik π²_k +^e_π^iK2

K

)

, k=j,

∑n i=1

eiK

π²_K, k̸=j.

=

∑n i=1

K∑−1 k=1

K∑−1 l=1

ek

e_iK π_K² e^t_l+

∑n i=1

K∑−1 k=1

ek

e_ik π²_ke^t_k. Becauseθ andπare orthogonal parameter sets, the full matrix is

I_c(y, z;θ, π) =

( I_c^θ(y, z;θ) 0 0 I_c^π(y, z;π)

) .

Note that the conditional expectations can be approximated by E(e_ik|y_i) =

∑1 j=0

jP(e_ik=j|y_i) =P(e_ik= 1|y_i) =πkf(yi|eik= 1) f(y_i) ≈w_ik.

By definition, the marginal score is the conditionally expected complete score. Hence, for k = 1, . . . , K−1

S^π(y;θ, π) =E (

S_c^π(y, z;π)y )

=

∑n i=1

K∑−1 k=1

e_kE_π (eik

π_k −eiK

π_K y_i

)

≈

∑n i=1

K∑−1 k=1

e_k (wik

π_k −wiK

π_K )

with respective (approximate) marginal estimates ˆ

π= 1 n

∑n i=1

K∑−1 k=1

e_kwˆ_ik.

The observed information is again described by the difference between the complete and the missing part. The observed θ information matrix part is handled in the same way as before.

Because theθscore does not depend onπ(and vice versa) the complete (θ, π) information matrix Ic(y;θ, π) consists of two blocks. Theπblock is calculated by means ofIc^π(y;π), andIm^π(y;π), the conditional variance of the completeπscores. Now we have thecomplete Fisher information

Ic^π(y;π) = E (

I_c^π(y, z;π)y )

=

∑n i=1

K∑−1 k=1

K∑−1 l=1

e_kE (e_iK

π_K² y_i

) e^t_l+

∑n i=1

K∑−1 k=1

e_kE (e_ik

π_k² y_i

) e^t_k

≈

∑n i=1

K∑−1 k=1

K∑−1 l=1

e_kwiK

π_K² e^t_l+

∑n i=1

K∑−1 k=1

e_kwik

π²_k e^t_k

and themissing information Im^π(y;π) =

∑n i=1

K∑−1 k=1

K∑−1 l=1

ekE ((eik

πk −eiK

πK

)(eil

πl −eiK

πK

)yi

) e^t_l

−

∑n i=1

K∑−1 k=1

K∑−1 l=1

ekE ((e_ik

πk −e_iK πK

)yi

) E

((e_il πl −e_iK

πK

)yi

) e^t_l.

Notice thatE(eikeil|yi) = 0 fork ̸=l and E(e²_ik|yi) =E(eik|yi). Hence, the missing information equals

Im^π(y;π) =

∑n i=1

K∑−1 k=1

ekE (eik

π²_k yi

) e^t_k+

∑n i=1

K∑−1 k=1

K∑−1 l=1

ekE (eiK

π_K² yi

) e^t_l

−

∑n i=1

K∑−1 k=1

K∑−1 l=1

ekE ((eik

πk −eiK

πK

)yi

) E

((eil

πl −eiK

πK

)yi

) e^t_l

≈

∑n i=1

K∑−1 k=1

ek

w_ik π_k² e^t_k+

∑n i=1

K∑−1 k=1

K∑−1 l=1

ek

w_iK π_K² e^t_l

−

∑n i=1

K∑−1 k=1

K∑−1 l=1

e_k (w_ik

πk

−w_iK πK

)(w_il πl

−w_iK πK

) e^t_l

giving asobserved information I^π(y;π) =

∑n i=1

K∑−1 k=1

K∑−1 l=1

ekE ((eik

π_k −eiK

π_K )yi

) E

((eil

π_l −eiK

π_K )yi

) e^t_l

≈

∑n i=1

K∑−1 k=1

K∑−1 l=1

ek

(wik

πk −wiK

πK

)(wil

πl −wiK

πK

) e^t_l.

3 The direct way – a reparameterized approach

We again assume that the conditional density is a member of the exponential family, i.e. for canonical link models with (extended) linear predictorηi=x^t_iγ+zi

f(yi|zi;θ)∝exp (

yiηi−b(ηi) )

(11) this gives conditional moments

E(yi|zi) =µi= ∂b(ηi)

∂ηi

var(y_i|z_i) =V(µ_i) =∂²b(ηi)

∂η²_i

and first derivatives

∂logf(y_i|z_i;θ)

∂ηi

= (y_i−µ_i) ∂logf(y_i|z_i;θ)

∂γ = (y_i−µ_i)x_i

∂²logf(yi|zi;θ)

∂η²_i =−V(µi) ∂²logf(yi|zi;θ)

∂γ∂γ^t =−V(µi)xix^t_i

If f(zi) does not depend on parameters, then the above results also hold for the log-likelihood based on the joint density logf(yi, zi;θ) = logf(yi|zi;θ) + logf(zi). If f(zi) is totally unknown, we can estimate through the discreteK-point distribution given by (ζk, πk),k= 1, . . . , K. We first reparametrize the massesπ=π(ϑ) as

π_k= exp(ϑ_k−κ(ϑ)), where ∂κ(ϑ)/∂ϑ_l=π_l (12) to ensureπ_k >0. The derivatives ofπw.r.tϑare therefore

∂π_k

∂ϑl

=

{ (1−π_k)π_k if k=l

−π_kπ_l if k̸=l. (13) This can be also written as

∂π_k

∂ϑl

=π_k (

I_(k=l)−π_l )

.

Now we study the derivative of the weights

wik=w(yi, ζk, γ, ϑ) = f(y_i|z_i;ζ_k, γ)π_k

∑K

l=1f(yi|zi;ζl, γ)πl

withf(yi|zi;ζk, γ) = exp (

yiηik−b(ηik) )

, where ηik=x^t_iγ+ζk. Therefore,

∂f(yi|zi;ζk, γ)

∂γ = xi(yi−µik)f(yi|zi;ζk, γ)

∂f(y_i|z_i;ζ_k, γ)

∂ζl

= I_(k=l)(y_i−µ_ik)f(y_i|z_i;ζ_k, γ) DefinefK(yi) =fK(yi;ζ, γ, ϑ) =∑K

k=1f(yi|zi;ζk, γ)πk. This gives

∂fK(yi)

∂γ =

∑K k=1

xi(yi−µik)f(yi|zi;ζk, γ)πk

∂f_K(y_i)

∂ζl

= (yi−µil)f(yi|zi;ζl, γ)πl

and

∂fK(yi)

∂ϑ_l =

∑K k=1

f(y_i|z_i;ζ_k, γ)∂πk

∂ϑ_l

= −

∑K k=1

f(y_i|z_i;ζ_k, γ)π_kπ_l+f(y_i|z_i;ζ_l, γ)π_lπ_l+f(y_i|z_i;ζ_l, γ)π_l(1−π_l)

= −πl

(

fK(yi)−f(yi|zi;ζl, γ) )

With these results we get

∂wik

∂γ = 1

fK(yi)

∂f(yi|zi;ζk, γ)

∂γ πk− 1

f_K²(yi)f(yi|zi;ζk, γ)πk

∂fK(yi)

∂γ

= x_i(y_i−µ_ik)f(y_i|z_i;ζ_k, γ)π_k fK(yi) −w_ik

∑K

l=1xi(yi−µil)f(yi|zi;ζl, γ)πl

fK(yi)

= xi(yi−µik)wik−wik

∑K l=1

xi(yi−µil)wil

∂wik

∂ζl

= 1

fK(yi)

∂f(yi|zi;ζk, γ)

∂ζl

πk− 1

f_K²(yi)f(yi|zi;ζk, γ)πk

∂fK(yi)

∂ζl

= I_(k=l)(y_i−µ_ik)f(y_i|z_i;ζ_k, γ)π_k fK(yi) −wik

(y_i−µ_il)f(y_i|z_i;ζ_l, γ)π_l fK(yi)

= I_(k=l)(yi−µik)wik−wikwil(yi−µil)

∂w_ik

∂ϑl

= 1

fK(yi)f(yi|zi;ζk, γ)∂π_k

∂ϑl − 1

f_K²(yi)f(yi|zi;ζk, γ)πk

∂f_K(y_i)

∂ϑl

= 1

f_K(y_i)f(yi|zi;ζk, γ)πk

(

I_(k=l)−πl

)

+ 1

f_K(y_i)wikπl

(

fK(yi)−f(yi|zi;ζl, γ) )

= wik(I_(k=l)−πl) +wikπl−wikwil

= wik

(

I(k=l)−wil

)

Letekbe a vector of zero with 1 at thekth position. The respective derivatives of the considered estimating equations

gγ(θ, ϑ) =

∑n i=1

∑K k=1

xi(yi−µik)wik

g_ζ(θ, ϑ) =

∑n i=1

∑K k=1

e_k(y_i−µ_ik)w_ik are

∂gγ(θ, ϑ)

∂γ^t =

∑n i=1

∑K k=1

x_i

(−V_ikw_ikx^t_i+ (y_i−µ_ik)∂wik

∂γ^t )