$\boldsymbol{y}_{hij}$, is given by the PLR model $$\pi _{hir}\left( \boldsymbol{\beta }\right) =\mathrm{E}\left[ Y_{hijr}|\boldsymbol{x}_{hi}\right] =\Pr \left( Y_{hijr}=1|\boldsymbol{x}_{hi}\right)
=\left\{ 
\begin{array}{ll}
\dfrac{\exp \{\boldsymbol{x}_{hi}^{T}\boldsymbol{\beta }_{r}\}}{1+{\sum_{l=1}^{d}}\exp \{\boldsymbol{x}_{hi}^{T}\boldsymbol{\beta }_{l}\}}, & 
r=1,...,d, \\ 
\dfrac{1}{1+{\sum_{l=1}^{d}}\exp \{\boldsymbol{x}_{hi}^{T}\boldsymbol{\beta }_{l}\}}, & r=d+1,\end{array}\right. ,  \label{2.1.0}$$with $\boldsymbol{\beta }_{r}=\left( \beta _{r0},\beta _{r1},...,\beta
_{rk}\right) ^{T}\in \mathbb{R}^{k+1}$, $r=1,...,d$. Note that, the expectation of $\boldsymbol{Y}_{hij}$ does not depends on the unit number $j$ (homogeneity), which is not a strong assumption as we generally have random sampling with the clusters in each stratum. Let $\boldsymbol{\pi }_{hi}\left( \boldsymbol{\beta }\right) $ denote the $(d+1)$-dimensional probability vector with the elements given in (\[2.1.0\]), i.e., $$\boldsymbol{\pi }_{hi}\left( \boldsymbol{\beta }\right) =\left( \pi
_{hi1}\left( \boldsymbol{\beta }\right) ,...,\pi _{hi,d+1}\left( \boldsymbol{\beta }\right) \right) ^{T}.  \label{2.2}$$Then, the associated parameter space for (\[2.1.0\]) is given by $$\Theta =\{\boldsymbol{\beta }=(\boldsymbol{\beta }_{1}^{T},...,\boldsymbol{\beta }_{d}^{T})^{T},\text{ }\boldsymbol{\beta }_{r}=\left( \beta
_{r0},...,\beta _{rk}\right) ^{T}\in \mathbb{R}^{k+1},\text{ }r=1,...,d\}=\mathbb{R}^{d(k+1)}.$$

For modeling data exhibiting overdispersion, the quasi-likelihood method is an widely used method, originally defined by Wedderburn (1974). The quasi-loglikelihood function is constructed without complete distributional knowledge but through the mean and the variance of independent sampling units, the total of clusters in all the strata $n=\sum_{h=1}^{H}n_{h}$ in the PLR model with complex sampling as given in (\[2.1.0\]). It is semiparametric method, since it only specifies the first two 