multivariate moments of $\boldsymbol{Y}_{hij}$. Let $f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ be the probability mass function of $\boldsymbol{Y}_{hij}$ such that $\boldsymbol{\pi }_{hi}\left( \boldsymbol{\beta }\right) $ is modeled by the PLR model with complex sampling, i.e. since the support of $\boldsymbol{Y}_{hij}$, say $\mathbb{Y}_{hi}$, is the set of column vectors of identity matrix $\boldsymbol{I}_{d+1}$ it holds $$\begin{aligned}
f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})& =\boldsymbol{\pi }_{hi}^{T}\left( \boldsymbol{\beta }\right) \boldsymbol{y}_{hij}={\sum_{l=1}^{d}}\pi _{hil}\left( \boldsymbol{\beta }\right) y_{hijl},
\\
\log f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})&
=\log \left( \boldsymbol{\pi }_{hi}^{T}\left( \boldsymbol{\beta }\right) 
\boldsymbol{y}_{hij}\right) ={\sum_{l=1}^{d}\log }\pi _{hil}\left( 
\boldsymbol{\beta }\right) y_{hijl}=\log \boldsymbol{\pi }_{hi}^{T}\left( 
\boldsymbol{\beta }\right) \boldsymbol{y}_{hij},\end{aligned}$$where we use the notation $\log \boldsymbol{\pi }_{hi}\left( \boldsymbol{\beta }\right) =\left( \log \pi _{hi1}\left( \boldsymbol{\beta }\right)
,...,\log \pi _{hi,d+1}\left( \boldsymbol{\beta }\right) \right) ^{T}$. It is important to be aware that the correct loglikelihood should be $
\mathcal{\ell (}\boldsymbol{\beta })={\sum\limits_{h=1}^{H}}{\sum\limits_{i=1}^{n_{h}}}\log f_{\boldsymbol{\beta }}(\underline{\boldsymbol{y}}_{hi}|\boldsymbol{x}_{hi}),
$ where $\underline{\boldsymbol{y}}_{hi}=(\boldsymbol{y}_{hi1},...,\boldsymbol{y}_{him_{hi}})^{T}$. Under homogeneity assumption within the clusters, since $f_{\boldsymbol{\beta }}(\underline{\boldsymbol{y}}_{hi}|\boldsymbol{x}_{hi}) $ is unknown, the quasi loglikelihood, $\mathcal{\ell (}\boldsymbol{\beta })$, considers the likelihood within each cluster the same as independent case as an approximation, $$\log f_{\boldsymbol{\beta }}(\underline{\boldsymbol{y}}_{hi}|\boldsymbol{x}_{hi})\overset{def}{=}{\sum\limits_{j=1}^{m_{hi}}}\log f_{\boldsymbol{\beta }}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi}),  \label{def0}$$and thus$$\mathcal{\ell (}\boldsymbol{\beta })={\sum\limits_{h=1}^{H}}{\sum\limits_{i=1}^{n_{h}}}{\sum\limits_{j=1}^{m_{hi}}}\log \boldsymbol{\pi }_{hi}^{T}\left( \boldsymbol{\beta }\right) \boldsymbol{y}_{hij}={\sum\limits_{h=1}^{H}}{\sum\limits_{i=1}^{n_{h}}}\log \boldsymbol{\pi }_{hi}^{T}\left( \boldsymbol{\beta }\right) \widehat{\boldsymbol{y}}_{hi},$$where$$\begin{aligned}
\widehat{\boldsymbol{y}}_{hi}& =\boldsymbol{1}_{m_{hi}}^{T}\underline{\boldsymbol{y}}_{hi}=\sum\limits_{j=1}^{m_{hi}}\boldsymbol{y}_{hij},
\label{2.201} \quad
\widehat{\boldsymbol{y}}_{hi} =(\widehat{y}_{hi1},...,\widehat{y}_{hi,d+1})^{T}=\left(
\sum\limits_{j=1}^{m_{hi}}y_{hij1},...,\sum\limits_{j=1}^{m_{hi}}y_{hij,d+1}\right) ^{T},  \end{aligned}$$denotes the realization value of counts in the $i$-th cluster of the $h$-th stratum.

