$\boldsymbol{u}\left( \boldsymbol{\beta }\right) =\boldsymbol{0}_{d(k+1)}$, where $$\begin{aligned}
\boldsymbol{u}\left( \boldsymbol{\beta }\right) &
=\sum\limits_{h=1}^{H}\sum\limits_{i=1}^{n_{h}}\boldsymbol{u}_{hi}\left( 
\boldsymbol{\beta },\boldsymbol{x}_{hi}\right) ,  \label{Un} \quad
\boldsymbol{u}\left( \boldsymbol{\beta },\boldsymbol{x}_{hi}\right)  =w_{hi}\left[ \widehat{\boldsymbol{Y}}_{hi}^{\ast }-m_{hi}\boldsymbol{\pi }_{hi}^{\ast }\left( \boldsymbol{\beta }\right) \right] \otimes \boldsymbol{x}_{hi},  \end{aligned}$$with superscript $^{\ast }$ here, the vector (matrix) obtained by deletingthe last row from the initial vector (matrix) is denoted; thus $\boldsymbol{\pi }_{hi}^{\ast }\left( \boldsymbol{\beta }\right) =\left( \pi _{hi1}\left( 
\boldsymbol{\beta }\right) ,...,\pi _{hid}\left( \boldsymbol{\beta }\right)
\right) ^{T}$ and $\widehat{\boldsymbol{Y}}_{hi}^{\ast }=\left( \widehat{Y}_{hi1}^{\ast
},...,\widehat{Y}_{hid}^{\ast }\right) ^{T}$. The derivation from (\[2.4\]) to (\[Un\]) is in Appendix and some additional details can be found in Morel (1989).

The Minimum Quasi Weighted Density Power Divergence estimators \[sec3\]
=======================================================================

Let $f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ be the probability mass function of $\boldsymbol{Y}_{hij}|\boldsymbol{x}_{hi}$ as defined in the previous section, $g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ an unknown and true probability mass function of $\boldsymbol{Y}_{hij}|\boldsymbol{x}_{hi}$ and $\mathbb{Y}_{hi}$ the support. The DPD based on the probability mass functions of a single observation of the sample, between $f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ and $g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$, is given by, for $\lambda >0$, $$\begin{aligned}
d_{\lambda }\left( g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi}),f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})\right) &
=\int\nolimits_{\mathbb{Y}_{hi}}\left( f_{\boldsymbol{\beta }}^{\lambda +1}(\boldsymbol{y}|\boldsymbol{x}_{hi})-\frac{\lambda +1}{\lambda }f_{\boldsymbol{\beta }}^{\lambda }(\boldsymbol{y}|\boldsymbol{x}_{hi})g(\boldsymbol{y}|\boldsymbol{x}_{hi})+\frac{1}{\lambda }g^{\lambda +1}(\boldsymbol{y}|\boldsymbol{x}_{hi})\right) d\boldsymbol{y} \\
& =\int\nolimits_{\mathbb{Y}_{hi}}f_{\boldsymbol{\beta }}^{\lambda }(\boldsymbol{y}|\boldsymbol{x}_{hi})dF(\boldsymbol{y}|\boldsymbol{x}_{hi})-\frac{\lambda +1}{\lambda }\int\nolimits_{\mathbb{Y}_{hi}}f_{\boldsymbol{\beta }}^{\lambda }(\boldsymbol{y}|\boldsymbol{x}_{hi})dG(\boldsymbol{y}|\boldsymbol{x}_{hi})+K \\
& =d_{\lambda }^{\ast }\left( g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi}),f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})\right) +K,\end{aligned}$$where $K$ is a constant not depending on $\boldsymbol{\beta 