Robust semiparametric inference for polytomous logistic regression with complex survey design

TL;DR

This paper introduces robust semiparametric estimators for polytomous logistic regression models under complex survey designs, improving inference robustness in socio-economic research.

Contribution

It develops a new class of minimum quasi weighted density power divergence estimators that enhance robustness over traditional methods for complex survey data.

Findings

Estimators exhibit strong robustness properties.

The methods outperform traditional likelihood-based approaches.

Empirical validation confirms theoretical advantages.

Abstract

Analyzing polytomous response from a complex survey scheme, like stratified or cluster sampling is very crucial in several socio-economics applications. We present a class of minimum quasi weighted density power divergence estimators for the polytomous logistic regression model with such a complex survey. This family of semiparametric estimators is a robust generalization of the maximum quasi weighted likelihood estimator exploiting the advantages of the popular density power divergence measure. Accordingly robust estimators for the design effects are also derived. Robust testing of general linear hypotheses on the regression coefficients are proposed using the new estimators. Their asymptotic distributions and robustness properties are theoretically studied and also empirically validated through a numerical example and an extensive Monte Carlo study.

Figures

Tables

Table 1. Table 2: Mean standardized deviations under contamination for BMI data

	$\lambda$
	0	0.2	0.4	0.6	0.8	1
contamination of Men (45–64 years)
$masd$(${\widehat{𝜷}}_{\lambda }^{\ast }$,${\widehat{𝜷}}_{\lambda }$)	0.24396	0.23057	0.21731	0.20441	0.19187	0.17969
$masd$(${\widehat{𝝅}}_{\lambda }^{\ast }$,${\widehat{𝝅}}_{\lambda }$)	0.10170	0.09700	0.09220	0.08750	0.08280	0.07810
$asd$(${\hat{\rho }}_{\lambda }^{2\ast ,E},{\hat{\rho }}_{\lambda }^{2,E}$)	0.64785	0.58974	0.52955	0.46868	0.40858	0.35044
$asd$(${\hat{\rho }}_{\lambda }^{2\ast ,M},{\hat{\rho }}_{\lambda }^{2,M}$)	1.57103	1.41605	1.25700	1.09724	0.94068	0.79037
contamination of Women (45–64 years)
$masd$(${\widehat{𝜷}}_{\lambda }^{\ast }$,${\widehat{𝜷}}_{\lambda }$)	0.10516	0.09484	0.08533	0.07665	0.0687	0.06148
$masd$(${\widehat{𝝅}}_{\lambda }^{\ast }$,${\widehat{𝝅}}_{\lambda }$)	0.03250	0.03030	0.02810	0.02600	0.0240	0.02210
$asd$(${\hat{\rho }}_{\lambda }^{2\ast ,E},{\hat{\rho }}_{\lambda }^{2,E}$)	0.07038	0.05132	0.03358	0.01744	0.00321	0.00893
$asd$(${\hat{\rho }}_{\lambda }^{2\ast ,M},{\hat{\rho }}_{\lambda }^{2,M}$)	0.03358	0.01169	0.05244	0.08819	0.11821	0.14206

Equations

$${\mathbf{y}}_{hij}={\left(y_{hij1},....,y_{hij,d+1}\right)}^T,h=1,...,H,i=1,...,n_h,j=1,...,m_{hi},$$

$$\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.,$$

$${\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)={\left({\pi }_{hi1}\left(\mathbf{\beta }\right),...,{\pi }_{hi,d+1}\left(\mathbf{\beta }\right)\right)}^T.$$

$$\Theta =\{\mathbf{\beta }=({\mathbf{\beta }}_1^T,...,{\mathbf{\beta }}_d^T{)}^T,{\mathbf{\beta }}_r={\left({\beta }_{r0},...,{\beta }_{rk}\right)}^T\in {\mathbb{R}}^{k+1},r=1,...,d\}={\mathbb{R}}^{d(k+1)}.$$

$$f_{\mathbf{\beta }}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})$$

$$\log f_{\mathbf{\beta }}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})$$

$$\log f_{\mathbf{\beta }}({\underline{\mathbf{y}}}_{hi}|{\mathbf{x}}_{hi})\stackrel{def}{=}\sum _{j=1}^{m_{hi}}\log f_{\mathbf{\beta }}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi}),$$

$$\ell (\mathbf{\beta })=\sum _{h=1}^H\sum _{i=1}^{n_h}\sum _{j=1}^{m_{hi}}\log {\mathbf{\pi }}_{hi}^T\left(\mathbf{\beta }\right){\mathbf{y}}_{hij}=\sum _{h=1}^H\sum _{i=1}^{n_h}\log {\mathbf{\pi }}_{hi}^T\left(\mathbf{\beta }\right){\hat{\mathbf{y}}}_{hi},$$

$${\hat{\mathbf{y}}}_{hi}$$

$${\mathbf{\Sigma }}_{hi}={\mathbf{\Sigma }}_{hi}\left(\mathbf{\beta }\right)=m_{hi}\mathbf{\Delta }({\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)),$$

$${\mathbf{\Sigma }}_{hi}={\mathbf{\Sigma }}_{hi}({\nu }_{hi},\mathbf{\beta })={\nu }_{hi}{m}_{hi}\mathbf{\Delta }({\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)),$$

$$\ell (\mathbf{\beta },w)=\sum _{h=1}^H\sum _{i=1}^{n_h}{w}_{hi}\log {\mathbf{\pi }}_{hi}^T\left(\mathbf{\beta }\right){\hat{\mathbf{Y}}}_{hi}.$$

$$\sum _{h=1}^H\sum _{i=1}^{n_h}{w}_{hi}\frac{\partial {\mathbf{\pi }}_{hi}^T\left(\mathbf{\beta }\right)}{\partial \mathbf{\beta }}{\mathrm{diag}}^{-1}({\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right))\left[{\hat{\mathbf{Y}}}_{hi}-m_{hi}{\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)\right]=0_{d(k+1)},$$

$$\frac{\partial {\mathbf{\pi }}_{hi}^T\left(\mathbf{\beta }\right)}{\partial \mathbf{\beta }}$$

$$\mathbf{u}\left(\mathbf{\beta }\right)$$

$$d_{\lambda }\left(g({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi}),f_{\mathbf{\beta }}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})\right)$$

$$={\int }_{{\mathbb{Y}}_{hi}}{f}_{\mathbf{\beta }}^{\lambda }(\mathbf{y}|{\mathbf{x}}_{hi})dF(\mathbf{y}|{\mathbf{x}}_{hi})-\frac{\lambda +1}{\lambda }{\int }_{{\mathbb{Y}}_{hi}}{f}_{\mathbf{\beta }}^{\lambda }(\mathbf{y}|{\mathbf{x}}_{hi})dG(\mathbf{y}|{\mathbf{x}}_{hi})+K$$

$$=d_{\lambda }^{\ast }\left(g({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi}),f_{\mathbf{\beta }}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})\right)+K,$$

$$d_{\lambda }^{\ast }\left(g({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi}),f_{\mathbf{\beta }}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})\right)=E[f_{\mathbf{\beta }}^{\lambda }({\mathbf{Y}}_{hij}|{\mathbf{x}}_{hi})]-\frac{\lambda +1}{\lambda }{\int }_{{\mathbb{Y}}_{hi}}{f}_{\mathbf{\beta }}^{\lambda }(\mathbf{y}|{\mathbf{x}}_{hi})dG(\mathbf{y}|{\mathbf{x}}_{hi})$$

$$d_{\lambda }^{\ast }\left(\hat{g}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi}),f_{\mathbf{\beta }}({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})\right)=E[f_{\mathbf{\beta }}^{\lambda }({\mathbf{Y}}_{hij}|{\mathbf{x}}_{hi})]-\frac{\lambda +1}{\lambda }{f}_{\mathbf{\beta }}^{\lambda }({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi}).$$

$$d_{\lambda }^{\ast }(\hat{g},f_{\mathbf{\beta }})=\sum _{h=1}^H\sum _{i=1}^{n_h}\sum _{j=1}^{m_{hi}}\left(E\left[f_{\mathbf{\beta }}^{\lambda }({\mathbf{Y}}_{hij}|{\mathbf{x}}_{hi})\right]-\frac{\lambda +1}{\lambda }{f}_{\mathbf{\beta }}^{\lambda }({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})\right),\text{for }\lambda >0.$$

$$E\left[f_{\mathbf{\beta }}^{\lambda }({\mathbf{Y}}_{hij}|{\mathbf{x}}_{hi})\right]=E\left[{\mathbf{\pi }}_{hi}^{\lambda ,T}\left(\mathbf{\beta }\right){\mathbf{Y}}_{hij}\right]=\sum _{l=1}^{d+1}{\pi }_{hil}^{\lambda +1}\left(\mathbf{\beta }\right).$$

$$d_{\lambda }^{\ast }(\hat{g},f_{\mathbf{\beta }},w)=\sum _{h=1}^H\sum _{i=1}^{n_h}\sum _{j=1}^{m_{hi}}{w}_{hi}\left(E\left[f_{\mathbf{\beta }}^{\lambda }({\mathbf{Y}}_{hij}|{\mathbf{x}}_{hi})\right]-\frac{\lambda +1}{\lambda }{f}_{\mathbf{\beta }}^{\lambda }({\mathbf{y}}_{hij}|{\mathbf{x}}_{hi})\right),\text{for }\lambda >0,$$

$$d_{\lambda }^{\ast }(\hat{g},f_{\mathbf{\beta }},w)=\sum _{h=1}^H\sum _{i=1}^{n_h}{w}_{hi}{\mathbf{\pi }}_{hi}^{\lambda ,T}\left(\mathbf{\beta }\right)\left(m_{hi}{\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)-\frac{\lambda +1}{\lambda }{\hat{\mathbf{y}}}_{hi}\right),\text{for }\lambda >0.$$

$${\hat{\mathbf{\beta }}}_{\lambda ,Q}=\arg \underset{\mathbf{\beta }\in {\mathbb{R}}^{d(k+1)}}{\min }{d}_{\lambda }^{\ast }(\hat{g},f_{\mathbf{\beta }},w).$$

$${\mathbf{u}}_{\lambda }(\mathbf{\beta })$$

$${\mathbf{u}}_{\lambda }(\mathbf{\beta },{\mathbf{x}}_{hi})$$

$${\mathbf{U}}_{\lambda }(\mathbf{\beta },{\mathbf{x}}_{hi})=\left[w_{hi}{\mathbf{\Delta }}^{\ast }({\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)){\mathrm{diag}}^{\lambda -1}\{{\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)\}\{{\hat{\mathbf{Y}}}_{hi}-m_{hi}{\mathbf{\pi }}_{hi}\left(\mathbf{\beta }\right)\}\right]\otimes {\mathbf{x}}_{hi}$$

$$\sqrt{n}({\hat{\mathbf{\beta }}}_{\lambda ,Q}-{\mathbf{\beta }}_0)\stackrel{\mathcal{L}}{\underset{n\to \infty }{⟶}}\mathcal{N}\left(0_{d(k+1)},{\mathbf{\Psi }}_{\lambda }^{-1}\left({\mathbf{\beta }}_0\right){\mathbf{\Omega }}_{\lambda }\left({\mathbf{\beta }}_0\right){\mathbf{\Psi }}_{\lambda }^{-1}\left({\mathbf{\beta }}_0\right)\right),$$

$${\mathbf{\Omega }}_{\lambda }\left(\mathbf{\beta }\right)$$

$${\mathbf{\Psi }}_{\lambda }\left(\mathbf{\beta }\right)$$

$${\mathbf{\Omega }}_{hi,\lambda }\left(\mathbf{\beta },{\mathbf{\Sigma }}_{hi}\right)$$

$${\mathbf{\Sigma }}_{hi}$$

Full text

Keywords: Cluster sampling; Design effect; Minimum quasi weighted DPD estimator; Polytomous logistic regression model; Pseudo minimum phi-divergence estimator; Quasi-likelihood; Robustness

1 Introduction

In many real-life applications, we come across data that have been collected through a complex survey scheme, like stratified sampling or cluster sampling, etc., rather than the simple random sampling. Such situations commonly arise in large scale data collection, for example, within several states of a country or even among different countries. Suitable statistical methods are required to analyze these data by taking care of the stratified structure of the data; this is because there often exist several inter and intra-class correlations within such stratification and ignoring them often lead to erroneous inference. Further, in many such complex surveys, stratified observations are collected on some categorical responses having two or more mutually exclusive unordered categories along with some related covariates and inference about their relationship is of up-most interest for insight generation and policy making. Polytomous logistic regression (PLR) model is a useful and popular tool in such situations to model categorical responses with associated covariates. However, most of classical literature deals with the cases of simple random sampling scheme (e.g. McCullagh, 1980; Lesaffre and Albert, 1989; Agresti, 2002; Gupta et al., 2006, 2008). The application of PLR model under complex survey setting can be found, for example, in Binder (1983), Roberts et al. (1987), Morel (1989), Morel and Neerchal (2012) and Castilla et al. (2018); most of them, except the last one, are based on the quasi maximum likelihood approach.

Even though the maximum quasi weighted likelihood estimator is the base of most of the existing literature on logistic models under complex survey designs, it is known to be non-robust in the presence of possible outliers in the data. In practice, with such a complex survey design, it is quite natural to have some outlying observations that make the likelihood based inference highly unstable. So, we often may need to make additional efforts to find and discard the outliers from the data before their analyses. A robust method providing stable solution even in presence of the outliers will be really helpful and more efficient in practice. The cited work by Castilla et al (2018) has developed an alternative minimum divergence estimator based on $\phi$-divergences (Pardo, 2006), but the important issue of robustness is still ignored there.

In this paper, we develop a robust estimator under the PLR model with a complex survey based on a minimum quasi weighted divergence approach. In particular, we exploit the nice properties of the density power divergence (DPD) of Basu et al. (1998). This measure become very popular in recent literature for yielding highly robust and efficient estimators under various statistical models; see, for example, Ghosh and Basu (2013, 2016, 2018) and, in particular, the recent paper by Castilla et al (2019) which discussed a PLR model but under the simple random sampling scheme.

We first start with the mathematical description of the PLR model with complex survey set-up and a brief discussion about the maximum quasi weighted likelihood estimator of the underlying parameters in Section 2. Then we introduce a class of new robust parameter estimates for the PLR model with complex survey by minimizing a suitably defined DPD measure in Section 3. The asymptotic distribution of the resulting estimator and the design effect for the PLR model with complex survey are also described there. In Section 4, a new family of Wald-type tests is introduced based on our new estimators for testing linear hypotheses about the parameters of the PLR model. In Section 5 we theoretically study the robustness of the proposed estimators and Wald-type tests through the influence function analysis. After presenting an illustrative example in Section 6, an extensive simulation study is presented in Section 7. The paper ends with a brief concluding remark in Section 8. For brevity in presentation, proofs of all the results are given in Appendix A.

2 Maximum Quasi Weighted Likelihood Estimator

Let us assume that the whole population is partitioned into $H$ distinct strata and the data consist of $n_{h}$ clusters in stratum $h$ for each $h=1,\ldots,H$. Further, for each cluster $i=1,\ldots,n_{h}$ in the stratum $h$, we have observed the values of a categorical response variable ($Y$) for $m_{hi}$ units. Assuming $Y$ has $(d+1)$ categories, we denote these observed responses by a $(d+1)$-dimensional classification vector

[TABLE]

with $y_{hijr}$ $=1$ and $y_{hijl}$ $=0$ for $l\in\{1,...,d+1\}-\{r\}$ if the $j$-th unit selected from the $i$-th cluster of the $h$-th stratum falls in the $r$-th category. We also have data on $(k+1)$ explanatory variables which are common for all the individuals in the $i$-th cluster of the $h$-th stratum (very common with dummy or qualitative explanatory variables) to be denoted as $\boldsymbol{x}_{hi}=\left(x_{hi0},x_{hi1},....,x_{hik}\right)^{T}$; the first one $x_{hi0}=1$ is associated with the intercept. Let us denote the sampling weight from the $i$-th cluster of the $h$-th stratum by $w_{hi}$. For each $i$, $h$ and $j$, the expectation of the $r$-th element of the random variable $\boldsymbol{Y}_{hij}=(Y_{hij1},...,Y_{hij,d+1})^{T}$, corresponding to the realization $\boldsymbol{y}_{hij}$, is given by the PLR model

[TABLE]

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), i.e.,

[TABLE]

Then, the associated parameter space for (2) is given by

[TABLE]

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). It is semiparametric method, since it only specifies the first two 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

[TABLE]

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,

[TABLE]

and thus

[TABLE]

where

[TABLE]

denotes the realization value of counts in the $i$-th cluster of the $h$-th stratum.

An important feature of the quasi loglikelihood is that the marginal distributions of $\boldsymbol{Y}_{hij}$ are completely known but the components of $\underline{\boldsymbol{Y}}_{hi}$, jointly, might be correlated. This means that distribution of their total, $\widehat{\boldsymbol{Y}}_{hi}$, might be also unknown, but the expectation is obtained as the total of the expectation of $\boldsymbol{Y}_{hij}$, $j=1,...,m_{hi}$. The most common assumption is to consider that $\widehat{\boldsymbol{Y}}_{hi}$ has a multinomial sampling scheme, which means that $\boldsymbol{Y}_{hij}$, $j=1,...,m_{hi}$ are independent random variables and

[TABLE]

where $\boldsymbol{\Delta}(\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right))=\mathrm{diag}(\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right))-\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right)\boldsymbol{\pi}_{hi}^{T}\left(\boldsymbol{\beta}\right)$. Since (4) is not an approximation, the term “quasi” should be dropped. A weaker assumption is to consider that $\widehat{\boldsymbol{Y}}_{hi}$ has a multinomial sampling scheme with a overdispersion parameter $\nu_{hi}=1+\rho_{hi}^{2}(m_{hi}-1)$, which means that $\boldsymbol{Y}_{hij}$, $j=1,...,m_{hi}$ are correlated random variables ($Cor[\boldsymbol{Y}_{hia},\boldsymbol{Y}_{hib}]=\rho_{hi}^{2}$, $a\neq b$, $a,b\in\{1,...,m_{hi}\}$) and

[TABLE]

but the distribution of $\widehat{\boldsymbol{Y}}_{hi}$ is not in principle used for the estimators. Distributions such as Dirichlet Multinomial, Random Clumped and $m$-inflated belong this family (see Alonso et al. (2017), Morel and Neerchal (2012) and Raim et al. (2015) for details). The weakest assumption is to consider that $\widehat{\boldsymbol{Y}}_{hi}$ has an unknown distribution, with $\boldsymbol{Y}_{hij}$, $j=1,...,m_{hi}$ being possibly correlated but with no specific pattern. It is worth of mentioning that $\nu_{hi}$ plays here a role of nuisance parameter and it is possible to consider a model with additional nuisance parameters, which are more complex than (7) but simpler than the option of completely unknown distribution for $\widehat{\boldsymbol{Y}}_{hi}$ (see Morel and Koehler (1995) for details).

Taking into account weights for each cluster, the quasi weighted loglikelihood is defined as

[TABLE]

Then, the maximum quasi weighted likelihood estimator of $\boldsymbol{\beta}$, say $\widehat{\boldsymbol{\beta}}_{P}$, is obtained by maximizing the quasi weighted loglikelihood, $\mathcal{\ell(}\boldsymbol{\beta},w)$, with respect to $\boldsymbol{\beta}$. The corresponding estimating equation is then given by

[TABLE]

with

[TABLE]

The system of equations (9) can be written as $\boldsymbol{u}\left(\boldsymbol{\beta}\right)=\boldsymbol{0}_{d(k+1)}$, where

[TABLE]

with superscript ∗ here, the vector (matrix) obtained by deleting the 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 (9) to (10) is in Appendix and some additional details can be found in Morel (1989).

3 The Minimum Quasi Weighted Density Power Divergence estimators

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$,

[TABLE]

where $K$ is a constant not depending on $\boldsymbol{\beta}$, $F(\boldsymbol{y}|\boldsymbol{x}_{hi})$ and $G(\boldsymbol{y}|\boldsymbol{x}_{hi})$ are the distribution functions corresponding to the densities $f(\boldsymbol{y}|\boldsymbol{x}_{hi})$ and $g(\boldsymbol{y}|\boldsymbol{x}_{hi})$, respectively, and

[TABLE]

is the kernel of $d_{\lambda}\left(g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})\text{, }f_{\boldsymbol{\beta}}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})\right)$. In practice, since $G$ is unknown, it must be estimated from the sample, which is in this case a single individual, so

[TABLE]

Based on Ghosh and Basu (2013), the “kernel of the ordinary DPD” between the probability mass functions $\widehat{g}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ and $f_{\boldsymbol{\beta}}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ for the whole sample is defined as a total discrepancy given by

[TABLE]

Since the support of $\boldsymbol{Y}_{hij}|\boldsymbol{x}_{hi}$, is the set of column vectors of identity matrix $\boldsymbol{I}_{d+1}$, it holds

[TABLE]

In the current paper it is defined for the first time the “kernel of the quasi weighted DPD” between the probability mass functions $\widehat{g}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ and $f_{\boldsymbol{\beta}}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ in the whole sample as a weighted sum

[TABLE]

whose expression for the PLR model with complex sampling is

[TABLE]

Based on $d_{\lambda}^{\ast}(\widehat{g},f_{\boldsymbol{\beta}},w)$ the minimum quasi weighted DPD estimator is formally defined as follows.

Definition 1

The minimum quasi weighted DPD estimator, say $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$, of $\boldsymbol{\beta}$ is defined as

[TABLE]

At the particular choice $\lambda\rightarrow 0$, the DPD measure, defined as the limit of (11) coincides (in limit) with the weighted quasi-loglikelihood, $\mathcal{\ell(}\boldsymbol{\beta},w)$, given in (8); thus the minimum quasi weighted DPD estimator of $\boldsymbol{\beta}$ at $\lambda=0$ coincides with the maximum weighted quasi likelihood estimator. With the same philosophy, the following result generalizes $\boldsymbol{u}(\boldsymbol{\beta},\boldsymbol{x}_{hi})$ given in (10) which plays an important role for the derivation of the asymptotic distribution of $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$.

Theorem 2

The minimum quasi weighted DPD estimate of $\boldsymbol{\beta}$, say $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$, can be obtained by solving the system of equations $\boldsymbol{u}_{\lambda}(\boldsymbol{\beta})=\boldsymbol{0}_{d(k+1)}$, where

[TABLE]

Let $\boldsymbol{U}_{\lambda}(\boldsymbol{\beta},\boldsymbol{X})$ be a random variable generator of (13) associated with a generic random explanatory variable $\boldsymbol{X}$, with no stratum and cluster assignment. In what is to follow,

[TABLE]

denotes $\boldsymbol{U}_{\lambda}(\boldsymbol{\beta},\boldsymbol{X})|\boldsymbol{X=x}_{hi}$. An important property is that the estimating equation given in Theorem 2 is unbiased. This property is very important to consider these estimators similar to the maximum quasi weighted likelihood ones in the construction of the asymptotic properties, which was not the case for other distance based estimators such as the one proposed in Castilla et al (2018).

3.1 Asymptotic distribution and estimates of the design effect

The following results are generalized for the PLR model with complex sampling and random explanatory variables. Without any loss of generality it is assumed that $\widehat{\Pr}(\boldsymbol{X=x}_{hi})=\frac{1}{n}$ is estimated from the sample of strata in order to get assymptotic results for the fixed explanatory variables (as a particular case of random explanatory variables).

Theorem 3

Let $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$ the minimum quasi weighted DPD estimate of $\boldsymbol{\beta}$ in the PLR model (2) under a complex survey. Then we have

[TABLE]

where $\boldsymbol{\beta}_{0}$ is the true parameter value and

[TABLE]

where

[TABLE]

Notice that the expression of $\mathbf{\Psi}_{\lambda=0}\left(\boldsymbol{\beta}\right)$ is the same as the so called Fisher information matrix for multinomial sampling. In addition, if $\widehat{\boldsymbol{Y}}_{hi}$ has a multinomial sampling scheme $\mathbf{\Omega}_{\lambda=0}\left(\boldsymbol{\beta}\right)=\mathbf{\Psi}_{\lambda=0}\left(\boldsymbol{\beta}\right)$ and $\mathbf{\Omega}_{\lambda=0}^{-1}\left(\boldsymbol{\beta}\right)\mathbf{\Psi}_{\lambda=0}\left(\boldsymbol{\beta}\right)\mathbf{\Omega}_{\lambda=0}^{-1}\left(\boldsymbol{\beta}\right)$ is the inverse of the Fisher information matrix.

Consistency is usually considered as a minimal requirement for an inference procedure. The following result is useful as a tool for estimating $\mathbf{\Omega}_{\lambda}\left(\boldsymbol{\beta}\right)$ and $\mathbf{\Psi}_{\lambda}\left(\boldsymbol{\beta}\right)$ consistently plugging a consistent estimator into $\boldsymbol{\beta}$, and hence also $\mathbf{\Omega}_{\lambda}^{-1}\left(\boldsymbol{\beta}\right)\mathbf{\Psi}_{\lambda}\left(\boldsymbol{\beta}\right)\mathbf{\Omega}_{\lambda}^{-1}\left(\boldsymbol{\beta}\right)$ again as a (double) plug-in estimator.

Corollary 4

The following ones are (weak) consistent estimators as $n$ goes to infinity:

a)

$\widehat{\boldsymbol{\beta}}_{\lambda,Q}$* is a consistent estimator of the true regression coefficient $\boldsymbol{\beta}_{0}$.*

b)

$\mathbf{\Psi}_{n,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})=\frac{1}{n}\sum_{h=1}^{H}\sum_{i=1}^{n_{h}}\mathbf{\Psi}_{hi,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})$* is a consistent estimator of $\mathbf{\Psi}_{\lambda}\left(\boldsymbol{\beta}_{0}\right)$.*

c)

$\mathbf{\Omega}_{n,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\{\widehat{\boldsymbol{\Sigma}}_{hi}\}_{h=1,...,H;i=1,...,n_{h}})=\frac{1}{n}\sum_{h=1}^{H}\sum_{i=1}^{n_{h}}\mathbf{\Omega}_{hi,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\widehat{\boldsymbol{\Sigma}}_{hi})$* with*

[TABLE]

is a consistent estimator of $\mathbf{\Omega}_{\lambda}\left(\boldsymbol{\beta}_{0}\right)$, whenever $\Sigma_{hi}=\mathrm{Var}[\widehat{\boldsymbol{Y}}_{hi}]$ is consistently estimated through $\widehat{\boldsymbol{\Sigma}}_{hi}$ for all $(h,i)\in\{1,...,H\}\times\{1,...,m_{hi}\}$.

The following two cases of $\mathbf{\Omega}_{n,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\{\widehat{\boldsymbol{\Sigma}}_{hi}\}_{h=1,...,H;i=1,...,n_{h}})$ are taken into account for two cases of $\boldsymbol{\Sigma}_{hi}$ for which there exists a consistent estimator:

• •

if $\widehat{\boldsymbol{Y}}_{hi}$ has (an ordinary) multinomial sampling scheme then

[TABLE]

where $\mathbf{\Omega}_{hi,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})$ is $\mathbf{\Omega}_{hi,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\widehat{\boldsymbol{\Sigma}}_{hi})$ given in (20) with

[TABLE]

• •

if $\widehat{\boldsymbol{Y}}_{hi}$ has an overdispersed multinomial sampling scheme then

[TABLE]

where $\widetilde{\nu}_{hi}$ is a consistent estimator of the overdispersion parameter $\nu_{hi}$ which is established later in Corollary 7, and $\mathbf{\Omega}_{n,\lambda}(\widetilde{\nu}_{hi},\widehat{\boldsymbol{\beta}}_{\lambda,Q})$ is $\mathbf{\Omega}_{hi,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\widehat{\boldsymbol{\Sigma}}_{hi})$ given in (20) with

[TABLE]

Remark 5

If $\widehat{\boldsymbol{Y}}_{hi}$ has a multinomial sampling scheme, Theorem 3 can be similarly formulated taking the assumption that $\widehat{\Pr}(\boldsymbol{X=x}_{hi})=\frac{1}{m}$, where $m=\sum_{h=1}^{H}\sum_{i=1}^{n_{h}}m_{hi}$, for each individual of all the clusters in the sample, and taking for the estimation equation $m$ summands rather than $n$, i.e. plugging ${\sum\limits_{j=1}^{m_{hi}}}\boldsymbol{y}_{hij}=\widehat{\boldsymbol{y}}_{hi}$ into (12) and considering the system of equations $\boldsymbol{u}_{\lambda}(\boldsymbol{\beta})=\boldsymbol{0}_{d(k+1)}$, where

[TABLE]

Hence, it holds

[TABLE]

with

[TABLE]

The formal proof is omitted, but the derivation of the expressions is almost the same considering $\boldsymbol{U}_{\lambda}(\boldsymbol{\beta},\boldsymbol{x}_{hi})={\sum\limits_{j=1}^{m_{hi}}}\boldsymbol{V}_{j,\lambda}(\boldsymbol{\beta},\boldsymbol{x}_{hi})$, with $\boldsymbol{V}_{j,\lambda}(\boldsymbol{\beta},\boldsymbol{x}_{hi})$ i.i.d. random variables $j=1,...,m_{hi}$. This idea matches the philosophy of the asymptotic result developed in Castilla et al. (2019), where for $H=1$ and $w_{1i}=1$, $i=1,...,m_{1i}$.

The following result is useful for any sample of polytomous logistic regression with complex sample design, more general in comparison with Corollary 4, since it is not necessary to get any consistent estimators for $\boldsymbol{\Sigma}_{hi}$.

Theorem 6

The estimator $\widehat{\mathbf{\Omega}}_{n,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})=\frac{1}{n}\sum_{h=1}^{H}\sum_{i=1}^{n_{h}}\widehat{\mathbf{\Omega}}_{hi,\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})$, with

[TABLE]

is consistent for $\mathbf{\Omega}_{\lambda}\left(\boldsymbol{\beta}\right)$ as $n$ goes to infinity.

Corollary 7

Let $\widehat{\boldsymbol{Y}}_{hi}$ be a random variable with overdisepersed multinomial sampling scheme with a common overdispersion parameter $\nu$ and $m_{hi}=\overline{m}$,

[TABLE]

then, for $\nu$ and $\rho^{2}$:

a)

“robust and consistent estimators based on the estimating equation” are given respectively by

[TABLE]

where the matrices of interest are as follows, the one associated with multinomial sampling

[TABLE]

and Theorem 6 for overdispersed multinomial sampling.

b)

“robust and consistent estimators based on the method of moments” are given by

[TABLE]

The preceding results were established for $n=\sum\nolimits_{h=1}^{H}n_{h}$ tending to infinity which implies in practice that $H$ is fixed and there exists $\eta_{h}=\lim_{n\rightarrow\infty}\frac{n_{h}}{n}\in(0,1)$, $\sum\nolimits_{h=1}^{H}\eta_{h}=1$. This means that it is approriate to consider in the place of Theorem 3, the following one, assuming that $\widehat{\Pr}(\boldsymbol{X}_{h}\boldsymbol{=x}_{hi})=\frac{1}{n_{h}}$.

Theorem 8

Let $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$ the quasi MDPDE of $\boldsymbol{\beta}$ in the PLR model (2) under a complex survey. Then we have

[TABLE]

where $\boldsymbol{\beta}_{0}$ is the true parameter value and

[TABLE]

where $\mathbf{\Omega}_{i,\lambda}^{(h)}\left(\boldsymbol{\beta},\boldsymbol{\Sigma}_{hi}\right)$ and $\mathbf{\Psi}_{i,\lambda}^{(h)}\left(\boldsymbol{\beta}\right)$ have the same expressions (not meaning) as $\mathbf{\Omega}_{hi,\lambda}\left(\boldsymbol{\beta},\boldsymbol{\Sigma}_{hi}\right)$ and $\mathbf{\Psi}_{hi,\lambda}\left(\boldsymbol{\beta}\right)$ respectively, given in (16), (19).

The idea of the previous theorem matches the philosophy of the asymptotic result developed in Castilla et al. (2018), where other family of divergences was considered and the derivations of the results were quite different. It is also appropriate to consider the following results in the place of Corollary 4, b and c.

Corollary 9

The following ones are (weak) consistent estimators as $n_{h}$ goes to infinity, for every $h=1,...,H$:

a)

$\mathbf{\Psi}_{\lambda}^{(h)}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})=\frac{1}{n_{h}}\sum_{i=1}^{n_{h}}\mathbf{\Psi}_{i,\lambda}^{(h)}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})$* is a consistent estimator of $\mathbf{\Psi}_{\lambda}^{(h)}\left(\boldsymbol{\beta}_{0}\right)$, for every $h=1,...,H$.*

b)

$\mathbf{\Omega}_{\lambda}^{(h)}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\{\widehat{\boldsymbol{\Sigma}}_{hi}\}_{h=1,...,H:i=1,...,n_{h}})=\frac{1}{n_{h}}\sum_{i=1}^{n_{h}}\mathbf{\Omega}_{i,\lambda}^{(h)}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\widehat{\boldsymbol{\Sigma}}_{hi})$* is a consistent estimator of $\mathbf{\Omega}_{\lambda}^{(h)}\left(\boldsymbol{\beta}_{0}\right)$, for every $h=1,...,H$, whenever $\Sigma_{hi}=\mathrm{Var}[\widehat{\boldsymbol{Y}}_{hi}]$ is consistently estimated through $\widehat{\boldsymbol{\Sigma}}_{hi}$ for all $(h,i)\in\{1,...,H\}\times\{1,...,m_{hi}\}$.*

The following result is useful for any sample of polytomous logistic regression with complex sample design, more general in comparison with Corollary 9, since it is not neccessary to get any consistent estimators for $\boldsymbol{\Sigma}_{hi}$.

Theorem 10

The estimator

[TABLE]

is consistent of $\mathbf{\Omega}_{\lambda}^{(h)}\left(\boldsymbol{\beta}\right)$, as $n_{h}$ goes to infinity, for every $h=1,...,H$.

Corollary 11

Let $\widehat{\boldsymbol{Y}}_{hi}$ be a random variable with overdisepersed multinomial sampling scheme with an overdispersion parameter $\nu_{h}$, specific for each stratum, and $m_{hi}=\overline{m}_{h}$,

[TABLE]

then, for $\nu_{h}$ and $\rho_{h}^{2}$:

a)

“robust and consistent estimators based on the estimating equation” are given respectively by

[TABLE]

where the matrices of interest are the one associated with multinomial sampling for the $h$-th stratum,

[TABLE]

and Theorem 10 for overdispersed multinomial sampling.

b)

“robust and consistent estimators based on the method of moments” are given by

[TABLE]

4 Testing linear hypotheses for the PLR coefficients with complex survey design

Based on the asymptotic distribution of the minimum quasi weighted DPD estimator, $\widehat{\boldsymbol{\beta}}_{\lambda,Q},$ presented in Theorem 8, we now introduce and study a family of Wald-type test statistics for testing

[TABLE]

where $\boldsymbol{M}$ is a $d(k+1)\times r$ full row-rank matrix with $r\leq d(k+1)$ and the right-hand-side $r-$vector consists of constants that in many situations are $\boldsymbol{l=0}_{r}.$ For solving the problem of testing given in (23) we define a family of Wald-type test statistics based on quasi MDPDE as follows.

Definition 12

Let $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$ the minimum quasi weighted DPD estimator of $\boldsymbol{\beta}$ in the PLR model (2) under a complex survey and we denote

[TABLE]

where

[TABLE]

The family of Wald-type test statistics for testing the null hypothesis given in (23) is given by

[TABLE]

For $\lambda=0,$ $\widehat{\boldsymbol{\beta}}_{\lambda=0,Q},$ is the maximum quasi weighted likelihood estimator of $\boldsymbol{\beta}$, with estimating equations given in (9). It is not difficult to see that $\mathbf{Q}_{\lambda=0}(\widehat{\boldsymbol{\beta}}_{\lambda=0,Q})$ is the Fisher information matrix and $W_{n}(\widehat{\boldsymbol{\beta}}_{\lambda=0,Q})$ will be the classical Wald test statistic.

Theorem 13

Under the null hypothesis given in (23) the asymptotic distribution of the Wald-type test statistics $W_{n}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})$, defined in (24), is chi-square with $r$ degrees of freedom.

The proof is immediate using the asymptotic distribution of the minimum quasi weighted DPD estimator, presented in Theorem 8, and taking into account the consistence of the matrix $\mathbf{Q}_{\lambda}\left(\widehat{\boldsymbol{\beta}}_{\lambda,Q}\right)$ presented in Corollary 9. Based on the previous theorem the null hypothesis given in (23) will be rejected if

[TABLE]

Next we present a result in order to give an approximation of the power function for the test statisitcs given in (25). Let $\boldsymbol{\beta}^{\ast}$ $\in\Theta$ be such that $\boldsymbol{M}^{T}\boldsymbol{\beta}^{\ast}\neq\boldsymbol{c,}$ i.e., $\boldsymbol{\beta}^{\ast}$ does not belong to the null parameter space. Let us denote

[TABLE]

Then we have the following result.

Theorem 14

Let $\boldsymbol{\beta}^{\ast}\in\Theta$ with $\boldsymbol{M}^{T}\boldsymbol{\beta}^{\ast}\neq\boldsymbol{l}$ the true value of the parameter so that $\widehat{\boldsymbol{\beta}}_{\lambda,Q}\underset{n\rightarrow\infty}{\overset{P}{\longrightarrow}}\boldsymbol{\beta}^{\ast}$. The power function of the test statistic given in (25), at $\boldsymbol{\beta}^{\ast}$, is given by

[TABLE]

where $\phi_{n}\left(x\right)$ tends uniformly to the standard normal distribution $\phi\left(x\right)$ and $\sigma\left(\boldsymbol{\beta}^{\ast}\right)$ is given by

[TABLE]

It is clear that

[TABLE]

for all $\alpha\in\left(0,1\right)$. Therefore, the Wald-type test statistics are consistent in the sense of Fraser.

Remark 15

Based on the previous theorem we can obtain the sample size necessary to get a fix power $Po_{W_{n}\left(\widehat{\boldsymbol{\beta}}_{\lambda,Q}\right)}\left(\boldsymbol{\beta}^{\ast}\right)=\pi_{0}.$ By (26) we must solve the equation

[TABLE]

and we get that $n=\left[n^{\ast}\right]+1$ with

[TABLE]

being $A=\sigma^{2}\left(\boldsymbol{\beta}^{\ast}\right)\left(\phi^{-1}\left(1-\pi_{0}\right)\right)^{2}\text{ and }B=\frac{1}{2}\chi_{r,\alpha}^{2}l_{\boldsymbol{\beta}^{\ast}}(\boldsymbol{\beta}^{\ast}).$

We may also find an approximation of the power of the Wald-type tests given in (25) at an alternative close to the null hypothesis. Let $\boldsymbol{\beta}_{n}$ such that $\boldsymbol{M}^{T}\boldsymbol{\mathbf{\beta}_{n}\neq l}$ be a given alternative and let $\boldsymbol{\mathbf{\beta}}_{0}$ be the element such that $\boldsymbol{M}^{T}\boldsymbol{\mathbf{\beta}_{0}}=\boldsymbol{l}$ closest to $\boldsymbol{\beta}_{n}$ in the Euclidean distance sense. A first possibility to introduce contiguous alternative hypotheses is to consider a fixed $\boldsymbol{d}\in\mathbb{R}^{d(k+1)}$ and to permit $\boldsymbol{\beta}_{n}$ moving towards $\boldsymbol{\beta}_{0}$ as $n$ increases in the following way

[TABLE]

A second approach is to relax the condition $\boldsymbol{M}^{T}\boldsymbol{\mathbf{\beta}_{0}}=\boldsymbol{l}$ defining the null hypothesis. Let $\boldsymbol{\delta}\in\mathbb{R}^{r}$ and consider the following sequence, $\boldsymbol{\beta}_{n},$ of parameters moving towards $\boldsymbol{\beta}_{0}$ according to

[TABLE]

Note that

[TABLE]

Then the equivalence between the two hypotheses is given by $M\boldsymbol{{}^{T}\boldsymbol{d}}=\boldsymbol{\delta}.$

If we denote by $\chi_{r}^{2}(\Delta)$ the non central chi-square distribution with $r$ degrees of freedom and noncentrality parameter $\Delta,$ we can state the following theorem.

Theorem 16

We have,

1. i)

$W_{n}\left(\widehat{\boldsymbol{\beta}}_{\lambda,Q}\right)\underset{n\longrightarrow\infty}{\overset{\mathcal{L}}{\longrightarrow}}\chi_{r}^{2}\left(\Delta_{1}\right)$under $H_{1,n}$ given in (27), with

[TABLE] 2. ii)

$W_{n}\left(\widehat{\boldsymbol{\beta}}_{\lambda,Q}\right)\underset{n\longrightarrow\infty}{\overset{\mathcal{L}}{\longrightarrow}}\chi_{r}^{2}\left(\Delta_{2}\right)$* under $H_{1,n}^{\ast}$ given in (28), with*

[TABLE]

5 Influence Function

Influence function is a classical tool to measure robustness of an estimator (Hampel et al., 1968). However, the present set-up of complex survey is not as simple as the iid set-up; in fact the observations within a cluster of a stratum are iid but the observations in different cluster and stratum are independent non-homogeneous. So we need to modify the definition of the influence function accordingly. Recently, Ghosh and Basu (2013, 2016, 2018) have discussed the extended definition of the influence function for the independent but non-homogeneous observations; we will extend their approach to define the influence function in the present case of PLR model under complex design.

5.1 Influence Function of the minimum quasi weighted DPD estimator

We first need to define the statistical functional corresponding to the minimum quasi weighted DPD estimator as the minimizer of the DPD between the true and model densities. Assume the set-up and notations of Section 3 with the DPD kernel being given by $d_{\lambda}^{\ast}\left(g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi}),f_{\boldsymbol{\beta}}(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})\right)$ for individual densities $g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ and $f_{\boldsymbol{\beta}}(\boldsymbol{y}_{hij})$. Then, following Ghosh and Basu (2013), the minimum quasi weighted DPD estimator functional is to be defined by the minimizer of the total weighted DPD measure given by

[TABLE]

Definition 17

We consider the PLR model with complex survey defined in (2). The minimum quasi weighted DPD estimator functional, ${\boldsymbol{T}}_{\lambda,Q}(\boldsymbol{G})$, of $\boldsymbol{\beta}$ at $\boldsymbol{G}=(G(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi}),~{}h=1,\ldots,H,i=1,\ldots,n_{h})$ is defined as

[TABLE]

where $d_{\lambda}(g,f_{\boldsymbol{\beta}},w)$ is as defined above in (LABEL:def_func).

Note that, by the property of the DPD measure, it is immediate that the minimum quasi weighted DPD estimator functional ${\boldsymbol{T}}_{\lambda,Q}(\boldsymbol{G})$ is Fisher consistent at the assumed PLR model (2), i.e., ${\boldsymbol{T}}_{\lambda,Q}(\boldsymbol{\pi}(\boldsymbol{\beta}))=\boldsymbol{\beta}$ for all parameter values $\boldsymbol{\beta}$. Also, following the calculations for Theorem 2, one can see that the minimum quasi weighted DPD estimator functional ${\boldsymbol{T}}_{\lambda,Q}(\boldsymbol{G})$ can also be derived as a solution to the estimating equations

[TABLE]

where ${g}_{hi}=g(\boldsymbol{y}_{hij}|\boldsymbol{x}_{hi})$ and

[TABLE]

Note that, these estimating equations are unbiased at the model probability $\boldsymbol{g}=\boldsymbol{\pi}(\boldsymbol{\beta})$.

Now, in order to define the influence function of the minimum DPD functional ${\boldsymbol{T}}_{\lambda,Q}(\boldsymbol{g})$, we note that the functional itself depends on the sample sizes and cluster weights and so it IF will have the same dependence in analogue to the non-homogeneous case of Ghosh and Basu (2013, 2016, 2018). Also note that, the contamination can be in any one particular cluster within one stratum or simultaneously in many of them (or all). For simplicity, let us first assume that the contamination is only in one cluster probability ${g}_{h_{0}i_{0}}$ for some fix $h_{0}$ and $i_{0}$. Consider the contaminated probability vector ${g}_{h_{0}i_{0},\epsilon}=(1-\epsilon){g}_{h_{0}i_{0}}+\epsilon\delta_{\boldsymbol{t}}$, where $\epsilon$ is the contamination proportion and $\delta_{\boldsymbol{t}}$ is the degenerate probability at the outlier point $\boldsymbol{t}=(t_{1},\ldots,t_{d+1})^{T}\in\{0,1\}^{d+1}$ with $\sum_{s=1}^{d+1}t_{s}=1$. Denote the corresponding contaminated full probability vector as $\boldsymbol{g}_{\epsilon}$ which is same as $\boldsymbol{g}$ except ${g}_{h_{0}i_{0}}$ being replaced by ${g}_{h_{0}i_{0},\epsilon}$ and let the corresponding contaminated distribution vector be $\boldsymbol{G}_{\epsilon}$. Then, the corresponding influence function is defined as

[TABLE]

In order to calculate this influence function, we start with the estimating equation for ${\boldsymbol{T}}_{\lambda,Q}$. Note that, ${\boldsymbol{T}}_{\lambda,Q}(\boldsymbol{G}_{\epsilon})$ satisfies the equations

[TABLE]

Now, differentiating it with respect to $\epsilon$ art $\epsilon=0$ and simplifying, we get the required influence function as given by

[TABLE]

where

[TABLE]

Note that, at the model probability $\boldsymbol{g}=\boldsymbol{\pi}(\boldsymbol{\beta})$, we have $\boldsymbol{u}_{\lambda,h_{0}i_{0}}^{*}\left(\boldsymbol{\pi}_{h_{0}i_{0}}(\boldsymbol{\beta}),\boldsymbol{\beta}\right)=\boldsymbol{0}$ and $\boldsymbol{\Psi}_{n,\lambda}^{*}(\boldsymbol{\pi}(\boldsymbol{\beta}),\boldsymbol{\beta})=\boldsymbol{\Psi}_{n,\lambda}(\boldsymbol{\beta})$, as defined in Theorem 3. Hence, the influence function of the proposed minimum quasi weighted DPD estimator at the model probability simplifies to

[TABLE]

Clearly, by the assumed form of the model probability $\boldsymbol{\pi}(\boldsymbol{\beta})$ this influence function is bounded for all $\lambda>0$ but unbounded at $\lambda=0$. So, the proposed minimum quasi weighted DPD estimators are expected to be robust for any $\lambda>0$ but the corresponding maximum quasi weighted likelihood estimator (at $\lambda=0$) is non-robust against data contamination.

Similarly, one can show that, if there is contamination in some of the clusters within some stratum indexed by ${h,i}\in\Gamma\subseteq\{h=1,\ldots,H;i=1,\ldots,n_{h}\}$ at the contamination points $\boldsymbol{t}_{hi}$, the corresponding influence function of the proposed minimum quasi weighted DPD estimator at $\boldsymbol{g}$ and the model $\pi(\boldsymbol{\beta})$ are respectively given by

[TABLE]

and

[TABLE]

The boundedness and robustness implications for these influence functions are exactly the same as before.

5.2 Influence Function of the Wald-type Tests in PLRM

Let us now study the robustness of proposed Wald-type test through the influence function of the corresponding test statistics in (24). Considering the notation and set-up of Section 2 and 3, and define

[TABLE]

Then the Wald-type test functional $W_{\lambda}(\boldsymbol{G})$ corresponding to (24) for testing the null hypothesis given in (23) at the true distribution vector $\boldsymbol{G}$ is defined as

[TABLE]

where ${\boldsymbol{T}}_{\lambda,Q}$ is the minimum quasi weighted DPD estimator functional as defined in the previous subsection. Also, let $\boldsymbol{\beta}_{0}$ denote the true null parameter value for the hypothesis in (23).

Let us now derive the influence function of the test functional $W_{\lambda}$. As before, here also the contamination can be any particular cluster and strata (a given $h_{0},i_{0}$) combination or in many (or all) of them. The influence function of general Wald-type tests under such non-homogeneous set-up has been extensively studied in Basu et al. (2018). Here, we follow the general theory of Basu et al. (2018) to conclude that the first order influence functions of $W_{\lambda}$, defined as the first order derivative of its value at the contaminated distribution with respect to $\epsilon$ at $\epsilon=0$, in either case of contamination become identically zero at the null distribution $\boldsymbol{g}=\boldsymbol{\pi}(\boldsymbol{\beta}_{0})$. Therefore, the first order influence function is not informative in this case of Wald-type tests, and we need to investigate the second order influence function of $W_{\lambda}$.

The second order influence function $IF^{(2)}$, which measures the second order approximation to the asymptotic bias due to infinitesimal contamination, is defined as the second order derivative of the value of $W_{\lambda}$ at the contaminated distribution with respect to the contamination proportion $\epsilon$ at $\epsilon=0$. Again, following Basu et al. (2018), we can derive these second order influence functions of the Wald-type tests in either case of contaminations; at the null distribution $\boldsymbol{g}=\boldsymbol{\pi}(\boldsymbol{\beta}_{0})$, they simplifies to

[TABLE]

for the two types of contamination as before in the case of estimator. Note that, the second order influence functions of the proposed Wald-type tests are a quadratic function of the corresponding influence functions of the minimum quasi weighted DPD estimator for any type of contamination. Therefore, the boundedness of the influence functions of minimum quasi weighted DPD estimator at $\lambda>0$ also indicates the boundedness of the influence functions of the Wald-type test functional $W_{\lambda}$ indicating their robustness against contamination in any cluster or stratum of the sample data.

6 Illustrative Example: BMI data set

Let us consider an illustrative real-life dataset on BMI which was previously studied in Castilla et al. (2018). This data set, obtained from CANSIM Canada’s database and presented in Table 1, shows the body mass indexes of population in Canada in the year 1994. Each person in the sample is divided by their body mass index category under the international standard: acceptable weight, overweight or obese. The data set consists on a stratified sample design with clusters nested on them, with the strata being three different age groups (20–34 years, 35–44 years and 45–64 years) and the genders (male or female) as the clusters. The qualitative explanatory variables are valid to distinguish the clusters within the strata. They are given by $\boldsymbol{x}_{h1}^{T}=(1,0)$, and $\boldsymbol{x}_{h2}^{T}=(0,1)$, $h=1,\dots,5$ for men and women, respectively.

To illustrate the robustness of minimum quasi weighted DPD estimators, we contaminate the BMI data by permuting the categories overweight and obese in the Men with age range 45–64 years. After obtaining the correpsonding minimum quasi weighted DPD estimators estimates, we compute the mean absolute standardized deviations ($mssd$) between the estimated parameters and corresponding estimated probabilities obtained for the modified and original data, as given by

[TABLE]

where, with superscript ∗ we denote the contaminated case. Assuming common overdispersion parameter, absolute standardized deviation ($asd$) of the two versions of the intra-cluster correlation estimator are also computed as in Corollary 7. Their values, as presented in Table 2, clearly show that the minimum quasi weighted DPD estimators become more robust as $\lambda$ increases. Similar results are obtained when permuting the categories for Women instead of Men, as can also be seen in Table 2.

7 Monte Carlo Simulation Study

7.1 Simulation Scheme

We develop a complex design extension of the simulation scheme previously studied by Castilla et al. (2019), where a simple sample design was considered. $H$ strata consisting of $n_{h}=n$ clusters having $m_{hi}=m$ units each, are taken. Three overdispersed multinomial distributions: the random-clumped (RC), the $m$-inflated ($m$-I) and the Dirichlet Multinomial (DM) distributions (Alonso et al., 2017) having the same parameters $\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}_{0}\right)$ and $\rho$, are then considered for $\widehat{\boldsymbol{Y}}_{hi}$; which are characterized by

[TABLE]

As in Castilla et al. (2019), we further assume that the outcome nominal variable $\boldsymbol{Y}$ has $d+1=3$ categories depending on $k=2$ explanatory variables through the PLR model probabilities

[TABLE]

with $\boldsymbol{\beta}=(\beta_{01},\beta_{11},\beta_{21},\beta_{02},\beta_{12},\beta_{22})^{T}=(0,-0.9,0.1,0.6,-1.2,0.8)^{T}$ and $\boldsymbol{x}_{hi}\overset{iid}{\sim}\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$ for all $i=1,\ldots,n$, $h=1,\dots,H$.

Considering different values of the intra-cluster correlation parameter ($\rho^{2}$), the number of clusters in each stratum ($n$) and the clusters sizes ($m$), we simulate data from different scenarios.

1. Scenario 1: $H=2$, $n\in\{10i\}_{i=4}^{15}$, $m=21$, $\rho^{2}=0.25$, RC distribution.
2. Scenario 1b: $H=2$, $n\in\{10i\}_{i=4}^{15}$, $m=21$, $\rho^{2}=0.50$, RC distribution.
3. Scenario 2: $H=2$, $n=60$, $m\in\{10i\}_{i=2}^{12}$, $\rho^{2}=0.25$, RC m-I and DM distributions.
4. Scenario 3: $H=2$, $n=60$, $m=21$, $\rho^{2}\in\{0.1i\}_{i=0}^{9}$, RC, m-I and DM distributions.

In order to study the robustness issue, these simulations are repeated under contaminated data having $7\%$ outliers. These outliers are generated by permuting the elements of the outcome variable, such that categories 1, 2, 3 are classified as categories 3, 1, 2 for the outlying observations.

7.2 Performance of minimum quasi weighted DPD estimators and $\rho^{2}$ estimates

For the above scenarios, we compute the minimum quasi weighted DPD estimator of $\boldsymbol{\beta}$, for different tuning parameters $\lambda\in\{0,0.2,0.4,0.6,0.8\}$, and the corresponding estimate of $\rho^{2}$, both for the methods of moments and the estimating equations. The root of the mean square error (RMSE) are then computed based on $1000$ replications (Figures 5–8).

While classical maximum quasi weighted likelihood estimator presents the best behaviour under pure data, minimum quasi weighted DPD estimators with $\lambda>0$ are much more robust. In particular, as $\lambda$ increases, the change on their behaviour is accentuated. This is independent to the sample size (Figures 5 and 6), the intra-cluster correlation and the distribution (Figure 7) considered.

Estimator of $\rho^{2}$ by the method of moments seems less precise than the method of estimating equations, with independence of the tuning parameter chosen (Figures 5 and 6). Best estimators of $\rho^{2}$ by the estimating equations are obtained from minimum quasi weighted DPD estimators with $\lambda>0$, both for pure and contaminated data and for any of the distributions considered (Figure 7). Error of the estimators of $\rho^{2}$ tends to be smaller with the DM and RC distributions in comparison with the mI distribution, as can be seen in Figure 1, where density plots based on $1,000$ replications are shown for $\rho^{2}=0.5,n=60,m=21$ and tuning parameter $\lambda=0.4$. These results on the design effect parameter are consistent with the previous work of Alonso et al. (2017).

Notice that $\boldsymbol{\beta}$ estimates are obtained through $nlm()$ procedure with tolerance $10^{-6}$ in the software R, used for the whole Monte Carlo study. In Figure 1 iterations needed to compute the corresponding minimum quasi weighted DPD estimators were also obtained, with not a significant difference among the different distributions.

7.3 Comparison of minimum quasi weighted DPD estimators with pseudo minimum phi-divergence estimators

In the previous Section we have illustrated how the quasi minimum quasi weighted DPD estimators present a much more robust behavior than the maximum quasi weighted likelihood estimator for the estimation of $\boldsymbol{\beta}$ when a moderate/large sample size was considered. In Castilla et al. (2018), another family of estimators, the pseudo minimum phi-divergence estimators, was defined in the PLR model with complex survey. In particular, pseudo minimum Cressie-Read divergence estimators for positive tuning parameter $\lambda$ were studied. Results of the simulation study suggested these estimators as an interesting alternative to maximum quasi weighted likelihood estimator in terms of efficiency when a small sample size and a large intra cluster correlation was considered. Nevertheless, the robustness of these estimators was not studied.

In this section we want to make a general empirical comparison between minimum quasi weighted DPD estimator and pseudo minimum phi-divergence estimators (with positive and negative tuning parameters in the Cressie-Read subfamily) of $\boldsymbol{\beta}$, in terms of efficiency and robustness. Although they were not studied in the cited paper of Castilla et al. (2018), phi-divergences with negative tuning parameter and, in particular, the Hellinger distance ($\lambda=-0.5$), are known to have good robustness properties in some models. See for example Beran (1977) Lindsay (1994) and Toma (2007). This seems to happen in the context of PLR model with complex survey too, when a moderate/large sample size is considered. Nevertheless, pseudo minimum phi-divergence estimators with positive tuning parameter present an important lack of robustness. This behaviour can be summarized in Figure 2.

Let us now consider Scenario 3 of the previous section. A comparison between quasi minimum quasi weighted DPD estimators with $\lambda\in\{0.4,0.8\}$, classical maximum quasi weighted likelihood estimator, and pseudo minimum phi-divergence estimators with $\lambda\in\{-0.5,-0.3,0.66\}$ is made. Hellinger distance is shown to be, by far, the best choice when a contaminated scheme with low intra-cluster correlation is considered. With medium/high intra-correlation the behaviour turns to be the opposite. Minimum quasi weighted DPD estimators present a more stable behaviour, both in pure and contaminated schemes, with independence to the correlation parameter.

The same divergences are considered now in Scenario 1 for the estimation of the parameter $\rho^{2}$. The estimation of $\rho^{2}$ with pseudo minimum phi-divergence estimators is made by the method of Binder (Castilla et al., 2018) while the estimation of $\rho^{2}$ with minimum quasi weighted DPD estimators is made by the method of the estimating equations. Pseudo minimum phi-divergence estimators with negative tuning parameter are not competitive neither in the pure nor in the contaminated schemes. Pseudo minimum phi-divergence estimator with tuning parameter $\lambda=2/3$ is a good alternative to the minimum quasi weighted DPD estimators in contaminated data. This estimator showed also a good behaviour in the paper of Castilla et al. (2018).

7.4 Performance of the Wald-type tests

With the same model as in Section 7.1 we now empirically study the robustness of the minimum quasi weighted DPD estimator based Wald-type tests for the PLR model. As in Castilla et al. (2018) we first study the level under the true null hypothesis $H_{0}:\beta_{02}=0.6$. For studying the power robustness, the true data generating parameter value is considered as $\beta_{02}=1.08$.

Under Scenario 1, and both for pure and contaminated data, we compute observed levels and powers (measured as the proportion of test statistics exceeding the corresponding chi-square critical value), as can be seen in Figure 4. Under pure data, all levels present a very similar behaviour, while power attains their best value for classical maximum quasi weighted likelihood estimator. Under contamination, minimum quasi weighted DPD estimator based Wald-type tests for $\lambda>0$ present a more robust behaviour than maximum quasi weighted likelihood estimator, in concordance with previous results.

8 Concluding Remarks

We have presented minimum quasi weighted DPD estimators, as a robust alternative to classical approaches, in the modeling of categorical responses with associated covariates through polytomous logistic regression models. A Wald-type family of tests based on these estimators is also presented, for the problem of testing linear hypotheses on regression coefficients. The robustness of both estimators and tests are theoretically justified in terms of the influence function, which is shown to be bounded for the new procedures. An extensive simulation study is also provided, to empirically illustrate their robustness along with a comparison with the pseudo minimum phi-divergence estimators of Castilla et al. (2018). The results clearly show how the proposed minimum quasi weighted DPD estimators seem to be the best choice for dealing with the robustness issue for a moderate sample size.

Appendix A Appendix

Derivation of (10) from (9)

Proof. From (9) it holds ${\sum\limits_{h=1}^{H}}{\sum\limits_{i=1}^{n_{h}}}\boldsymbol{u}_{hi}(\boldsymbol{\beta})=\boldsymbol{0}_{d(k+1)}$, with

[TABLE]

The first term of $\boldsymbol{u}_{hi}(\boldsymbol{\beta})$ is$w_{hi}\boldsymbol{\Delta}^{\ast}(\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right))\mathrm{diag}^{-1}\{\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right)\}\{\widehat{\boldsymbol{y}}_{hi}-m_{hi}\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right)\},$ where

[TABLE]

but $\boldsymbol{1}^{T}\{\widehat{\boldsymbol{y}}_{hi}-m_{hi}\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right)\}=0$ and hence

[TABLE]

Proof of Theorem 2

Proof. The minimum quasi weighted DPD estimator of $\boldsymbol{\beta}$, is defined as

[TABLE]

which can also be obtained by solving the system of equations $\boldsymbol{u}_{\lambda}\left(\boldsymbol{\beta}\right)=\boldsymbol{0}_{d(k+1)}$ where

[TABLE]

Now, taking into account that

[TABLE]

we get

[TABLE]

and hence the system of equations becomes

[TABLE]

Proof of Theorem 3

Proof. By following formulas (3.3) and (3.4) of Ghosh and Basu (2013), it holds

[TABLE]

which can be expressed as a limit of

[TABLE]

where the summands given in (14) and (15) are

[TABLE]

With respect to $\mathbf{\Omega}_{\lambda}\left(\boldsymbol{\beta}\right)=\mathrm{Var}\left[\boldsymbol{U}_{\lambda}(\boldsymbol{\beta},\boldsymbol{X})\right]=\mathrm{E}\left[\mathrm{Var}\left[\boldsymbol{U}_{\lambda}(\boldsymbol{\beta},\boldsymbol{X})|\boldsymbol{X}\right]\right]+\mathrm{Var}\left[\mathrm{E}\left[\boldsymbol{U}_{\lambda}(\boldsymbol{\beta},\boldsymbol{X})|\boldsymbol{X}\right]\right]$, the first summand estimated from the sample is

[TABLE]

and the second one

[TABLE]

for $\boldsymbol{\beta}=\boldsymbol{\beta}_{0}$, since $\mathrm{E}\left[\boldsymbol{U}_{\lambda}(\boldsymbol{\beta}_{0},\boldsymbol{x}_{hi})\right]=\boldsymbol{0}_{d(k+1)}$ from $\mathrm{E}[\widehat{\boldsymbol{Y}}_{hi}-m_{hi}\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}_{0}\right)]=\boldsymbol{0}_{d(k+1)}$. On the other hand,

[TABLE]

With respect to $\mathbf{\Psi}_{\lambda}\left(\boldsymbol{\beta}\right)=-\mathrm{E}\left[\frac{\partial}{\partial\boldsymbol{\beta}}\boldsymbol{U}_{\lambda}^{T}(\boldsymbol{\beta},\boldsymbol{X})\right]$, estimated from the sample is given by

[TABLE]

where the third equality is again true for $\boldsymbol{\beta}=\boldsymbol{\beta}_{0}$. The derivations are omitted for $\lambda=0$, in which case similar ideas as $\lambda>0$ might be followed.

Proof of Theorem 14

Proof. We have

[TABLE]

Now we are going to get the asymptotic distribution of the random variable $\sqrt{n}\left(l_{\widehat{\boldsymbol{\beta}}_{\lambda,Q}}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})-l_{\boldsymbol{\beta}^{\ast}}(\boldsymbol{\beta}^{\ast})\right)$.

It is clear that $l_{\widehat{\boldsymbol{\beta}}_{\lambda,Q}}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})$ and $l_{\widehat{\boldsymbol{\beta}}_{\lambda,Q}}(\boldsymbol{\beta}^{\ast})$ have the same asymptotic distribution because $\widehat{\boldsymbol{\beta}}_{\lambda,Q}\underset{n\rightarrow\infty}{\overset{P}{\longrightarrow}}\boldsymbol{\beta}^{\ast}.$ A first Taylor expansion of $l_{\widehat{\boldsymbol{\beta}}_{\lambda,Q}}(\boldsymbol{\beta}^{\ast})$ at $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$ around $\boldsymbol{\beta}^{\ast}$ gives

[TABLE]

Therefore,

[TABLE]

being

[TABLE]

Proof of Corollary 4

Proof. It is based on the weak consistency of $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$ as well as in the continuity with respect to $\boldsymbol{\beta}$ of the elements of different matrices. Based on Theorem 3, it holds $\lim_{n\rightarrow\infty}\mathrm{E}[\widehat{\boldsymbol{\beta}}_{\lambda,Q}]=\boldsymbol{0}_{d(k+1)}$ and $\lim_{n\rightarrow\infty}\mathrm{Var}[\widehat{\boldsymbol{\beta}}_{\lambda,Q}]=\boldsymbol{O}_{d(k+1)\times d(k+1)}$ and hence $\widehat{\boldsymbol{\beta}}_{\lambda,Q}\underset{n\rightarrow\infty}{\overset{P}{\longrightarrow}}\boldsymbol{\beta}_{0}$.

Proof of Theorem 6

Proof. Let

[TABLE]

be another possible consistent estimator of $\mathbf{\Omega}_{\lambda}\left(\boldsymbol{\beta}\right)$, alternative to $\mathbf{\Omega}_{n,\lambda}\left(\boldsymbol{\beta}\right)$. Since

[TABLE]

and $\mathrm{E}[\boldsymbol{U}_{\lambda}(\widehat{\boldsymbol{\beta}}_{\lambda,Q},\boldsymbol{X})|\boldsymbol{x}_{hi}]=0$, for being $\mathrm{E}[\widehat{\boldsymbol{Y}}_{hi}-m_{hi}\boldsymbol{\pi}_{hi}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})]=\boldsymbol{0}_{k+1}$, for all $(h,i)\in\{1,...,H\}\times\{1,...,m_{hi}\}$ according to the PLR model with complex design, it holds

[TABLE]

is another possible consistent estimator of $\mathbf{\Omega}_{\lambda}\left(\boldsymbol{\beta}\right)$, alternative to $\mathbf{\Omega}_{n,\lambda}\left(\boldsymbol{\beta}\right)$.

Proof of Corollary 7

Proof. Section a) is straightforward considering the expressions of matrix $\mathbf{G}_{n,\lambda}(\boldsymbol{\beta})$ and the related ones. For section b) we consider the vector

[TABLE]

Taking into account that $\mathrm{Var}[\widehat{\boldsymbol{Y}}_{hi}^{\ast}]=m\nu\boldsymbol{\Delta}\left(\boldsymbol{\pi}_{hi}^{\ast}\left(\beta\right)\right),$ we have

[TABLE]

for $h=1,...,H$ and $i=1,...,n_{h}$. An unbiased estimator of $\mathrm{Var}\left[\boldsymbol{Z}_{hi}^{\ast}\left(\boldsymbol{\beta}\right)\right]$ is

[TABLE]

i.e. $\mathrm{E}\left[\widehat{\mathrm{Var}}\left[\boldsymbol{Z}_{hi}^{\ast}\left(\boldsymbol{\beta}\right)\right]\right]=\mathrm{Var}\left[\boldsymbol{Z}_{hi}^{\ast}\left(\boldsymbol{\beta}\right)\right]$ for which the trace is

[TABLE]

Since $\widehat{\boldsymbol{\beta}}_{\lambda,Q}$ is consistent of $\boldsymbol{\beta}$, it holds that $\frac{1}{nd}\sum\limits_{h=1}^{H}\sum\limits_{i=1}^{n_{h}}\boldsymbol{Z}_{hi}^{\ast}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})\boldsymbol{Z}_{hi}^{\ast T}(\widehat{\boldsymbol{\beta}}_{\lambda,Q})$ is consistent of $\nu$, but

[TABLE]

is in principle different in shape in comparison with (22) with

[TABLE]

where $\Delta^{-}\left(\boldsymbol{\pi}_{hi}^{\ast}\left(\boldsymbol{\beta}\right)\right)=\mathrm{diag}^{-1}\left(\boldsymbol{\pi}_{hi}\left(\boldsymbol{\beta}\right)\right)$. To assure that both expressions are equivalent we count with the result related to invariance of quadratic forms with generalized variances, Lemma 1a, of Moore (1977), from which is concluded that

[TABLE]

Proof of Theorem 8

Proof. The details are omitted for being very similar to Theorem 3 conditioning within each stratum, i.e.

[TABLE]

Proof of Theorem 16

Proof. We have

[TABLE]

Therefore,

[TABLE]

We know, under $H_{1,n}$ that

[TABLE]

and

[TABLE]

The Wald-type test statistics can be written as the quadratic form $W_{n}\left(\widehat{\boldsymbol{\beta}}_{\lambda,Q}\right)=\boldsymbol{Z}^{T}\boldsymbol{Z}$ with

[TABLE]

and

[TABLE]

where $\boldsymbol{I}$ is the identity $r\times r$ matrix. The application of the previous gives i). The noncentrality parameter is $\boldsymbol{d}^{T}\boldsymbol{M}\left[M\boldsymbol{{}^{T}\mathbf{\Psi}_{\lambda}\left(\boldsymbol{\beta}_{0}\right)M}\right]^{-1}\boldsymbol{M^{T}d.}$ Result ii) follows using relation (29).