Hypotheses Testing from Complex Survey Data Using Bootstrap Weights: A Unified Approach

ABSTRACT

Standard statistical methods without taking proper account of the complexity of a survey design can lead to erroneous inferences when applied to survey data due to unequal selection probabilities, clustering, and other design features. In particular, the Type I error rates of hypotheses tests using standard methods can be much larger than the nominal significance level. Methods incorporating design features in testing hypotheses have been proposed, including Wald tests and quasi-score tests that involve estimated covariance matrices of parameter estimates. In this article, we present a unified approach to hypothesis testing without requiring estimated covariance matrices or design effects, by constructing bootstrap approximations to quasi-likelihood ratio statistics and quasi-score statistics and establishing its asymptotic validity. The proposed method can be easily implemented without a specific software designed for complex survey sampling. We also consider hypothesis testing for categorical data and present a bootstrap procedure for testing simple goodness of fit and independence in a two-way table. In simulation studies, the Type I error rates of the proposed approach are much closer to their nominal significance level compared with the naive likelihood ratio test and quasi-score test. An application to an educational survey under a logistic regression model is also presented. Supplementary materials for this article are available online.

KEYWORDS:

• Quasi-likelihood-ratio test
• Quasi-score test
• Wald test
• Wilks’ theorem

1 Introduction

Testing statistical hypotheses is one of the fundamental problems in statistics. In the parametric model approach, hypothesis testing can be implemented using Wald test, likelihood ratio test, or score test. In each case, a test statistic is calculated and then compared with the $100(1-\alpha )\%$-quantile of the reference distribution, which is the limiting distribution of the test statistic under the null hypothesis, where α is the nominal significance level. The limiting distribution is often a Chi-squared distribution (Shao 2003).

However, statistical inference with survey data involves additional steps incorporating the sampling design features. Korn and Graubard (1999) and Chambers and Skinner (2003) provided comprehensive overviews of the methods for analyzing survey data from complex sampling designs. In hypothesis testing with survey data, the limiting distribution of the test statistic is generally not a Chi-squared distribution. Rather, it can be expressed as a weighted sum of several independent random variables from ${\chi }^2(1)$, which is a Chi-square distribution with one degree of freedom, and the weights depend on unknown parameters. To handle such problems, one may consider some corrections to the test statistics to obtain a Chi-square limiting distribution approximately. Such an approach usually involves computing “design effects” associated with test statistics (Lumley and Scott 2014).

In this article, we use a different approach to computing the limiting distribution using a bootstrap approximation. Beaumont and Bocci (2009) investigated a general weighted bootstrap method for hypothesis testing under complex sampling designs in the context of Wald tests for linear regression analysis. In this article, we generalize the bootstrap testing idea of Beaumont and Bocci (2009) and present a unified bootstrap weight approach to obtain the limiting distribution of test statistics under complex sampling designs, including stratified multi-stage sampling with clusters selected with replacement.

The proposed method is a bootstrap weight method (Mashreghi, Haziza, and Léger 2016) under complex survey sampling. Although Praestgaard and Wellner (1993) and Chatterjee and Bose (2005) rigorously investigate bootstrap weight methods under regularity conditions, both may provide erroneous inferences under complex survey sampling. Praestgaard and Wellner (1993) established central limit theorems for the weighted bootstrap empirical process, but they assumed that the sample is independently and identically distributed (iid). Although it is common to assume that the finite population is iid, elements in a sample may not be due to complex sampling designs; see Pfeffermann (1993) for details. In addition, Praestgaard and Wellner (1993) only focused on mean estimation, but we are interested in hypothesis testing for parameters that solve a certain estimating equation. Chatterjee and Bose (2005) proposed a bootstrap weight method for estimating functions, and their theoretical results are valid even for correlated random variables. However, the bootstrap weight method of Chatterjee and Bose (2005) is not applicable to informative sampling designs, since survey weights are not taken into consideration; see Pfeffermann (1993) for details.

There exist bootstrap weight methods under complex survey sampling, but the proposed bootstrap method differs from them in the following aspects. Rao, Wu, and Yue (1992) proposed a bootstrap weight method to estimate the variance of a function of population total estimators under stratified random sampling, but did not investigate the theoretical properties of their bootstrap method. Chipperfield and Preston (2007) proposed a without replacement scaled bootstrap to achieve the same goal as Rao, Wu, and Yue (1992) under stratified random sampling, but their method is only applicable when the parameter of interest is a smooth function of population totals. Moreover, neither method is applicable to other complex survey sampling. Bertail and Combris (1997) proposed a bootstrap weight method to make inferences for population means under general sampling designs, and the corresponding theoretical properties were investigated. However, their results are not applicable to parameters obtained by solving an estimating equation. Although Beaumont and Patak (2012) developed a general bootstrap weight method for parameter solving an estimating equation in general sampling designs, they focused only on variance estimation without a rigorous discussion about hypothesis testing. Antal and Tillé (2011) proposed a one-one bootstrap method for variance estimation under general sampling designs. Similar to the one-one bootstrap method, Antal and Tillé (2014) developed a doubled half-bootstrap to estimate the variance of the estimator for population means. However, neither Antal and Tillé (2011) nor Antal and Tillé (2014) considered parameters solving an estimating equation. Moreover, existing work did not establish limiting distributions for the corresponding bootstrap estimator, especially for the one solving an estimating equation. Thus, up to our knowledge, there is no theoretical guarantee for existing methods when used for hypothesis testing.

Since the use of bootstrap for hypothesis testing is not fully investigated in the literature, our goal is to fill this important gap under general sampling designs. To do this, we first establish a bootstrap central limit theorem under regularity conditions. The sampling design is allowed to be informative in the sense that the sampling mechanism depends on the response variables being sampled (Pfeffermann 1993). Once the bootstrap central limit theorem is established, the proposed bootstrap method can be applied to the quasi-likelihood-ratio test and the quasi-score test. As long as the bootstrap central limit theorem holds for the sampling design, the bootstrap replicates of the test statistics can be used to approximate their distributions under the null hypothesis. Besides, this article also has a practical contribution to survey sampling, since bootstrap weights are commonly provided in agencies, including Statistics Canada.

The article is organized as follows. In Section 2, the basic setup is introduced. The proposed bootstrap method is presented in Section 3. In Section 4, the quasi-likelihood ratio test and the quasi-score test are introduced. Test for categorical survey data is briefly discussed in Section 5. The results of two simulation studies are presented in Section 5. An application to data from an educational survey under a logistic regression model is presented in Section 6. Concluding remarks are made in Section 7.

2 Basic Setup

Suppose that a finite population ${\mathcal{F}}_N=\{({\mathit{z}}_i,{\mathit{x}}_i,y_i):i=1,\dots ,N\}$ is a random sample of size N from a super-population model (Isaki and Fuller 1982; Pfeffermann 1993; Sverchkov and Pfeffermann 2004; Rubin-Bleuer and Schiopu-Kratina 2005), where y_i is the response of interest, and ${\mathit{x}}_i$ is the corresponding auxiliary vector of length d, and ${\mathit{z}}_i$ is the vector of design variables such that the first-order inclusion probability for unit i can be obtained by ${\pi }_i=\pi ({\mathit{z}}_i;{\varphi }_N)$, where ${\varphi }_N$ is a design parameter determined by $\{{\mathit{z}}_1,\dots ,{\mathit{z}}_N\}$. In other literature, the inclusion probability is expressed as ${\pi }_i=\pi ({\mathcal{F}}_N,i)$; see Sverchkov and Pfeffermann (2004) for details.

We are interested in hypothesis testing for a parameter ${\mathit{\theta }}_0=\arg {\max }_{\mathit{\theta }\in \Theta }E\{l(\mathit{\theta };\mathit{X},Y)\}$, where $\Theta \subset {\mathbb{R}}^p$ is the parameter space, $(\mathit{X},Y)$ is the random variable associated with elements of ${\mathcal{F}}_N$, the expectation is taken with respect to the super-population model, and $l(\mathit{\theta };\mathit{x},y)$ is a pre-defined objective function. For example, if the finite population is a random sample generated from a super-population model with a density function $f(\mathit{x},y;{\mathit{\theta }}_0)$ and if ${\mathit{\theta }}_0$ is the parameter of interest, then $l(\mathit{\theta };\mathit{x},y)=\hspace{0.17em}\log \hspace{0.17em}f(\mathit{x},y;\mathit{\theta })$. For regression problems, ${\mathit{\theta }}_0$ can be the parameter in the conditional expectation of the response of interest given the auxiliary vectors.

From the finite population ${\mathcal{F}}_N$, a probability sample $\{({\mathit{x}}_i,y_i):i\in A\}$ of size n is selected by a probability sampling design with the first-order inclusion probability π_i, where A is the set of sample indexes. The design variables ${\mathit{z}}_i$ and the design parameter ${\varphi }_N$ are not available at the sample level. From the sample, we can compute a pseudo M-estimator of ${\mathit{\theta }}_0$ by $$\widehat{\mathit{\theta }}=\arg \underset{\mathit{\theta }\in \Theta }{\max }{l}_w(\mathit{\theta }),$$where $l_w(\mathit{\theta })=N^{-1}\sum _{i\in A}{w}_{i}l(\mathit{\theta };{\mathit{x}}_i,y_i)$ is the survey-weighted objective function and $w_i={\pi }_i^{-1}$ is the sampling weight for $i\in A$. Note that $l_w(\mathit{\theta })$ is design-unbiased for $l_N(\mathit{\theta })=N^{-1}{\sum }_{i=1}^{N}l(\mathit{\theta };{\mathit{x}}_i,y_i)$ in the sense that $E\{l_w(\mathit{\theta })|{\mathcal{F}}_N\}=l_N(\mathit{\theta })$, where the conditional expectation conditioning on ${\mathcal{F}}_N$ is with respect to the sampling mechanism. To be rigorous, a subscript N should also be used to index the sample, the sampling weights, the estimator and other quantities, but it is suppressed in the sequel for simplicity.

Remark 1.

When discussing a finite population in the preceding paragraph, we implicitly assume that elements of ${\mathcal{F}}_N$ are iid. Such an assumption is widely adopted in survey sampling; see Corollary 1.3.2.1 and Section 2.2.1 of Fuller (2009) for details. Besides, under a linear regression model, (1) $$Y={\mathit{X}}^{\mathrm{T}}{\mathit{\beta }}_0+ϵ$$(1) with $E(ϵ|\mathit{x})=0$, the finite population $\{({\mathit{x}}_i,y_i):i=1,\dots ,N\}$ satisfies the iid assumption; see Särndal, Swensson, and Wretman (1989, sec. 5.1), Breidt and Opsomer (2000), Särndal, Swensson, and Wretman (2003, sec. 6.4), Kim et al. (2006, sec. 5), Toth and Eltinge (2011), Wu and Thompson (2020, sec. 5.1) and Han and Wellner (2021) for more general setups about the super-population model. If we are interested in the regression parameter ${\mathit{\beta }}_0$ in (1), an objective function is $l(\mathit{\theta };\mathit{x},y)=-{(y-{\mathit{x}}^{\mathrm{T}}\mathit{\beta })}^2$.

Although we implicitly make an iid assumption for the finite population, elements in the sample A may not be iid due to complex survey designs. That is, the informative sampling no longer preserves the iid structure in the sample; see Pfeffermann (1993) for details. Under certain sampling designs, we can relax the iid assumption, and the super-population model can involve cluster random effects (Rubin-Bleuer and Schiopu-Kratina 2005); refer to Section S4 of supplementary material for an extension to stratified multi-stage sampling with cluster random effects. We focus on element sampling mainly in this article.

Since the parameter of interest is the superpopulation parameter, the final sample can be viewed as a realization of the two-phase sampling design, where the finite population is the first-phase sample and the final sample is the second-phase sample. Chen and Rao (2007) and Rubin-Bleuer and Schiopu-Kratina (2005) present rigorous asymptotic setups for two-phase sampling.

Often, the pseudo M-estimator $\widehat{\mathit{\theta }}$ can be obtained by solving (2) $${\widehat{\mathit{S}}}_w(\mathit{\theta })=\frac{\partial }{\partial \mathit{\theta }}{l}_w(\mathit{\theta })=\frac{1}{N}\sum _{i\in A}{w}_i\mathit{S}(\mathit{\theta };{\mathit{x}}_i,y_i)=\mathbf{0},$$(2) where $\mathit{S}(\mathit{\theta };\mathit{x},y)=\partial l(\mathit{\theta };\mathit{x},y)/\partial \mathit{\theta }$ is the tangent function of $\mathit{\theta }$. To discuss the asymptotic properties of the pseudo M-estimator, we assume the following regularity conditions for a sequence of finite populations and samples (Isaki and Fuller 1982).

C1. The parameter space Θ is an open and convex set containing ${\mathit{\theta }}_0$ as an interior point. For $\mathit{\theta }\in \Theta ,l(\mathit{\theta };\mathit{x},y)$ is concave and twice-continuously differentiable with respect to $\mathit{\theta },E\left\{\right|l(\mathit{\theta };\mathit{X},Y)\left|\right\}<\infty$, and $E\{l(\mathit{\theta };\mathit{X},Y)\}$ is uniquely maximized at ${\mathit{\theta }}_0$.

C2. For $\mathit{\theta }\in \Theta ,V\{l_w(\mathit{\theta })|{\mathcal{F}}_N\}\to 0$ in probability with respect to the super-population model as $n\to \infty$, where $V\{l_w(\mathit{\theta })|{\mathcal{F}}_N\}$ is the design variance of $l_w(\mathit{\theta })$ conditional on the finite population ${\mathcal{F}}_N$. In addition, there exists a nonnegative design variance estimator $\widehat{V}\{l_w(\mathit{\theta })|{\mathcal{F}}_N\}$, such that $E[\widehat{V}\{l_w(\mathit{\theta })\left|{\mathcal{F}}_N\right\}|{\mathcal{F}}_N]=V\{l_w(\mathit{\theta })|{\mathcal{F}}_N\}$ almost surely.

C3. There exists a compact set $\mathcal{B}\subset \Theta$ containing ${\mathit{\theta }}_0$ as an interior point, such that $\underset{\mathit{\theta }\in \mathcal{B}}{\sup }\Vert n\mathit{V}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })|{\mathcal{F}}_N\}-{\mathbf{\Sigma }}_S(\mathit{\theta })\Vert \to 0$ in probability with respect to the super-population model as $n\to \infty$, where $\mathit{V}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })|{\mathcal{F}}_N\}$ is the design variance of ${\widehat{\mathit{S}}}_w(\mathit{\theta })$ conditional on the finite population, $\Vert \mathit{M}\Vert =\underset{\mathit{x}\in {\mathbb{R}}^m,\Vert \mathit{x}\Vert =1}{\sup }\Vert \mathit{Mx}\Vert$ is the operator norm of an m × m matrix $\mathit{M}$ associated with the Euclidean vector norm $\Vert ·\Vert$, and ${\mathbf{\Sigma }}_S(\mathit{\theta })$ is a positive-definitive nonstochastic matrix.

C4. A central limit theorem holds for the survey-weighted score function ${\widehat{\mathit{S}}}_w({\mathit{\theta }}_0)$ in (2). Specifically, $$n^{1/2}{\widehat{\mathit{S}}}_w({\mathit{\theta }}_0)\to N(\mathbf{0},{\mathbf{\Sigma }}_S({\mathit{\theta }}_0))$$in distribution as $n\to \infty$ with respect to the super-population model and the sampling mechanism.

C5. The weak convergence of ${\widehat{\mathit{I}}}_w(\mathit{\theta })=N^{-1}\sum _{i\in A}{w}_i\mathit{I}(\mathit{\theta };{\mathit{x}}_i,y_i)$ is uniform in $\mathit{\theta }\in \mathcal{B}$, where $\mathit{I}(\mathit{\theta };\mathit{x},y)=-\partial \mathbf{S}(\mathit{\theta };\mathit{x},y)/\partial {\mathit{\theta }}^{\mathrm{T}}$. Specifically, $\underset{\mathit{\theta }\in \mathcal{B}}{\sup }\Vert {\widehat{\mathit{I}}}_w(\mathit{\theta })-\mathcal{I}(\mathit{\theta })\Vert \to 0$ in probability with respect to the super-population model and the sampling mechanism as $n\to \infty$, where $\mathcal{I}(\mathit{\theta })=E\{\mathit{I}(\mathit{\theta };\mathit{X},Y)\}$, and $\mathcal{I}({\mathit{\theta }}_0)$ is invertible.

Conditions C1–C2 are sufficient conditions for the weak consistency of $\widehat{\mathit{\theta }}$. Conditions C3–C5 are used to establish the asymptotic normality of $\widehat{\mathit{\theta }}$. In Conditions C3–C4, we consider the case where $V[E\{{\widehat{S}}_w(\theta )\left|F_N\right\}]$ is asymptotically negligible compared with $E[\mathit{V}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })\left|{\mathcal{F}}_N\right\}]$, and it is often the case when $n/N\to 0$ as $n\to \infty$. That is, we can safely ignore the variability for generating the finite population and focus on the design variance of ${\widehat{\mathit{S}}}_w(\mathit{\theta })$ given the finite population. This assumption is a building block of the proposed bootstrap method in the next section, and a similar assumption is also proposed by Beaumont and Bocci (2009); see Assumption 5 and its discussion of Beaumont and Bocci (2009) for details.

As mentioned by Han and Wellner (2021), “establishing a finite-dimensional CLT … can be a nontrivial problem for general sampling designs,” and it is beyond our scope to prove Condition C4 under general sampling designs. In fact, there is no guarantee that the central limit theorem holds for general sampling designs. Instead, we simply consider the set of sampling designs such that the central limit theorem holds; see Hájek (1964), Rosén (1972a), Rosén (1972b), Berger (1998), Rubin-Bleuer and Schiopu-Kratina (2005), sections 1.3.3–1.3.4 of Fuller (2009), Chauvet (2015), and Han and Wellner (2021) for details. For simplicity, we use the same notation, “$\Vert ·\Vert$”, to denote the matrix norm and the vector norm. Detailed comments on these regularity conditions are presented in Section S1 of supplementary material. For the proposed conditions, we implicitly assume that $N\to \infty$ as $n\to \infty$. The following theorem establishes the asymptotic normality of $\widehat{\mathit{\theta }}$. Its proof can be found in Section S2 of supplementary material.

Theorem 1.

Under Conditions C1–C5, we have (3) $$\sqrt{n}(\widehat{\mathit{\theta }}-{\mathit{\theta }}_0)\to N(\mathbf{0},{\mathbf{\Sigma }}_{\mathit{\theta }})$$(3) in distribution as $n\to \infty$ with respect to the joint distribution of the super-population model and the sampling mechanism, where ${\mathbf{\Sigma }}_{\mathit{\theta }}=\mathbf{ℐ}{({\mathit{\theta }}_0)}^{-1}{\mathbf{\Sigma }}_S({\mathit{\theta }}_0){\left\{\mathbf{ℐ}{({\mathit{\theta }}_0)}^{-1}\right\}}^{\mathrm{T}}.$

Now, by the second-order Taylor expansion, we obtain (4) $$\begin{array}{cc}{l}_w({\mathit{\theta }}_0)\hfill & =l_w(\widehat{\mathit{\theta }})+{\widehat{\mathit{S}}}_w^{\mathrm{T}}(\widehat{\mathit{\theta }})({\mathit{\theta }}_0-\widehat{\mathit{\theta }})-\frac{1}{2}{(\widehat{\mathit{\theta }}-{\mathit{\theta }}_0)}^{\mathrm{T}}{\widehat{\mathit{I}}}_w(\widehat{\mathit{\theta }})(\widehat{\mathit{\theta }}-{\mathit{\theta }}_0)\hfill \\ \hfill & \hspace{1em}+o_p(n^{-1})\hfill \\ \hfill & =l_w(\widehat{\mathit{\theta }})-\frac{1}{2}{(\widehat{\mathit{\theta }}-{\mathit{\theta }}_0)}^{\mathrm{T}}{\widehat{\mathit{I}}}_w(\widehat{\mathit{\theta }})(\widehat{\mathit{\theta }}-{\mathit{\theta }}_0)+o_p(n^{-1}),\hfill \end{array}$$(4) where ${\widehat{\mathit{I}}}_w(\mathit{\theta })=N^{-1}\sum _{i\in A}{w}_i\mathit{I}(\mathit{\theta };{\mathit{x}}_i,y_i).$ Define (5) $$W({\mathit{\theta }}_0)=-2n\{l_w({\mathit{\theta }}_0)-l_w(\widehat{\mathit{\theta }})\}.$$(5)

By (4), we obtain $$W({\mathit{\theta }}_0)=n{(\widehat{\mathit{\theta }}-{\mathit{\theta }}_0)}^{\mathrm{T}}{\widehat{\mathit{I}}}_w(\widehat{\mathit{\theta }})(\widehat{\mathit{\theta }}-{\mathit{\theta }}_0)+o_p(1).$$

Thus, by Theorem 1, we have (6) $$W({\mathit{\theta }}_0)\to \mathcal{G}=\sum _{i=1}^p{c}_i{Z}_i^2$$(6) in distribution as $n\to \infty$, and the reference distribution is the joint distribution of the super-population model and the sampling mechanism, where $c_1,\dots ,c_p$ are the eigenvalues of $\mathit{D}={\mathbf{\Sigma }}_{\mathit{\theta }}\mathbf{ℐ}({\mathit{\theta }}_0)$, and $Z_1,\dots ,Z_p$ are p independent random variables from the standard normal distribution. Result (6) was established by Rao and Scott (1984) and Lumley and Scott (2014), and it can be regarded as a version of the Wilks’ theorem for survey sampling. For p = 1, we can use $c_1^{-1}W({\mathit{\theta }}_0)$ as the test statistic with ${\chi }^2(1)$ distribution as the limiting distribution under the null hypothesis, where $c_1={\sigma }_S^2({\theta }_0)/I({\theta }_0)$, and ${\sigma }_S^2(\theta )$ and $I(\theta )$ are the counterparts of ${\mathbf{\Sigma }}_s(\mathit{\theta })$ and $\mathit{I}(\mathit{\theta })$, respectively. Unless the sampling design is simple random sampling and the sampling fraction is negligible, the limiting distribution does not reduce to the standard Chi-squared distribution with p degrees of freedom.

3 Bootstrap Calibration

We propose using a bootstrap weight method to approximate the limiting distribution in (6) and other test statistics under complex survey sampling. Such a bootstrap calibration is very attractive in practice since there is no need to derive the analytic form of the limiting distribution of the test statistic. That is, standard software that takes account of the weights can be used to implement bootstrap hypothesis testing from data files reporting weights and associated sets of bootstrap weights, as done in Statistics Canada and some other agencies. Given the sample A, the proposed bootstrap method is as follows:

• Step 1 Generate a rescaling vector ${\mathit{r}}^{*}$ of size $n=\left|A\right|$ from a bootstrap distribution with $E_{*}({\mathit{r}}^{*})={\mathbf{1}}_n$, and the covariance of ${\mathit{r}}^{*}$ is determined by the sampling design, where $\left|A\right|$ is the cardinality of set A, $E_{*}(·)$ is the expectation with respect to the bootstrap distribution conditional on the sample A, and ${\mathbf{1}}_n$ is a vector of 1 with length n. The bootstrap weights are chosen in such a way that $E_{*}\{{\mathit{S}}_w^{*}(\mathit{\theta })\}={\widehat{\mathit{S}}}_w(\mathit{\theta })$ and ${\mathit{V}}_{*}\{{\mathit{S}}_w^{*}(\mathit{\theta })\}=\widehat{\mathit{V}}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })|{\mathcal{F}}_N\}$, where ${\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })=\sum _{i\in A}{w}_i^{*}\mathit{S}(\mathit{\theta };{\mathit{x}}_i,y_i),w_i^{*}=r_i^{*}{w}_i,{\mathit{V}}_{*}(·)$ is the bootstrap variance conditional on the sample A, and $\widehat{\mathit{V}}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })|{\mathcal{F}}_N\}$ is a nonnegative design variance estimator of ${\widehat{\mathit{S}}}_w(\mathit{\theta })$ conditional on the finite population ${\mathcal{F}}_N$.

• Step 2 Obtain a bootstrap replicate of $\widehat{\mathit{\theta }}$, denoted as ${\widehat{\mathit{\theta }}}^{*}$, by solving ${\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })=\mathbf{0}$.

Different distributions may be used to generate ${\mathit{r}}^{*}$ under different sampling designs, but the generation of ${\mathit{r}}^{*}$ only depends on the inclusion probabilities and not on the specific form of the estimating function for parameter estimation. Thus, once the bootstrap weights are generated, they can be used for various inference problems based on the sample A. However, the commonly used linearization technique is based on the specific form of the estimating function, so different linearization is required for different hypothesis testing problems. Thus, compared with the linearization technique, the proposed bootstrap method is more appealing to the practitioners.

We now present additional regularity conditions to validate the proposed bootstrap method.

C6. For $\mathit{\theta }\in \Theta$, the bootstrap weights are nonnegative and satisfy $E_{*}\{{\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })\}={\widehat{\mathit{S}}}_w(\mathit{\theta })$ and ${\mathit{V}}_{*}\{{\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })\}=\widehat{\mathit{V}}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })\}$ conditional on the sample A, where $\widehat{\mathit{V}}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })\}$ satisfies $$n\underset{\mathit{\theta }\in \mathcal{B}}{\sup }\Vert \widehat{\mathit{V}}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })\}-\mathit{V}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })\}\Vert \to 0$$in probability with respect to the super-population model and the sampling mechanism as $n\to \infty$.

C7. A central limit theorem holds for the bootstrap weighted score function ${\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })$ uniformly over $\mathcal{B}$. Specifically, (7) $$n^{1/2}\{{\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })-{\widehat{\mathit{S}}}_w(\mathit{\theta })\}|A \to N(\mathbf{0},{\mathbf{\Sigma }}_S(\mathit{\theta }))$$(7) in distribution uniformly over $\mathcal{B}$, and the reference distribution is the bootstrap sampling distribution conditional on the sample A.

C8. $\underset{\mathit{\theta }\in \mathcal{B}}{\sup }\Vert {\widehat{\mathit{I}}}_w^{*}(\mathit{\theta })-{\widehat{\mathit{I}}}_w(\mathit{\theta })\Vert \to 0$ in probability with respect to the bootstrap sampling distribution conditional on the sample A as $n\to \infty$, where ${\widehat{\mathit{I}}}_w^{*}(\mathit{\theta })=N^{-1}\sum _{i\in A}{w}_i^{*}\mathit{I}(\mathit{\theta };{\mathit{x}}_i,y_i)$.

Condition C6 guarantees that the scaled bootstrap variance of ${\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })$ converges to ${\mathbf{\Sigma }}_S(\mathit{\theta })$ in Condition C3, and they are also used to show the bootstrap central limit theorem in Condition C7. Condition C8 is needed to establish the asymptotic normality of the bootstrap estimator ${\widehat{\mathit{\theta }}}^{*}$. To check Condition C7, we may use a bootstrap version of the Lyapounov condition. As discussed in Condition C4, we do not claim that Condition C7 holds under general sampling designs, but we can validate both under certain sampling designs. For example, in Section S4 of the supplementary material, we use a Lyapounov-type condition to prove the two central limit theorems under stratified multi-stage sampling. Detailed comments on these regularity conditions are presented in Section S1 of the supplementary material. The high-level conditions are also verifiable under other commonly used sampling designs, including Poisson sampling.

Theorem 2.

Let the assumptions for Theorem 1 hold. Then, under Conditions C6–C8, we have (8) $$\sqrt{n}({\widehat{\mathit{\theta }}}^{*}-\widehat{\mathit{\theta }})|A \to N(\mathbf{0},{\mathbf{\Sigma }}_{\mathit{\theta }})$$(8) in distribution as $n\to \infty$, and the reference distribution on the left side of (8) is the bootstrap sampling distribution conditional on the sample A, where ${\mathbf{\Sigma }}_{\mathit{\theta }}$ is defined in Theorem 1.

Theorem 2 is proved under general sampling designs; see Section S3 of the supplementary material for details. Theorem 2 establishes the bootstrap central limit theorem for ${\widehat{\mathit{\theta }}}^{*}$, and the limiting distribution coincides with the one in (3) for $\widehat{\mathit{\theta }}$. Thus, the proposed bootstrap method can be used to approximate the sampling distribution of $\widehat{\mathit{\theta }}$.

For estimating the mean of a finite population, Chao and Lo (1985) first proved a bootstrap central limit theorem under simple random sampling without replacement, and Bickel and Freedman (1984) proved the bootstrap central limit theorem under stratified random sampling without replacement. We now present stratified multi-stage sampling with clusters selected with replacement as an example to validate the bootstrap central limit theorem in Theorem 2. Stratified multi-stage sampling with clusters selected with replacement is extensively used in socio-economic surveys when the first-stage sampling fractions within strata are negligible, as is often the case with large-scale sample surveys.

Example 1.

Under stratified multi-stage sampling, a popular approach is to assume that $n_h(\ge 2)$ clusters are selected with replacement from the hth stratum, with the selection probabilities p_hi satisfying ${\sum }_{i=1}^{N_h}{p}_{hi}=1$ for $h=1,\dots ,H$ (Rao and Wu 1988), where N_h is the number of clusters in the hth stratum. As mentioned in Section 2, the selection probabilities can be random with respect to the super-population model. Within each sampled cluster, denote ${\widehat{\mathit{S}}}_{w,hi}(\mathit{\theta })$ as a design-unbiased estimator of the cluster mean under sampling at the second and subsequent stages. For simplicity, we still use ${\widehat{\mathit{S}}}_{w,hi}(\mathit{\theta })$ for the case when the (hi)th cluster is fully observed. Then, the weighted score function is ${\widehat{\mathit{S}}}_w(\mathit{\theta })={\sum }_{h=1}^H{W}_h{\widehat{\mathit{S}}}_{w,h}(\mathit{\theta })$, where $W_h=M_h/M,M={\sum }_{h=1}^H{M}_h,M_h={\sum }_{i=1}^{N_h}{M}_{hi}$, M_hi is the size of the (hi)th cluster, ${\widehat{\mathit{S}}}_{w,h}(\mathit{\theta })=n_h^{-1}{\sum }_{i=1}^{n_h}{\tilde{\mathit{S}}}_{w,hi}(\mathit{\theta }),{\tilde{\mathit{S}}}_{w,hi}(\mathit{\theta })=d_{a(hi)}{p}_{a(hi)}^{-1}{\widehat{\mathit{S}}}_{w,a(hi)}(\mathit{\theta }),d_{a(hi)}=M_{a(hi)}/M_h$, and a(hi) is the index of the ith selected cluster in the hth stratum.

The bootstrap replicate of ${\widehat{\mathit{S}}}_w(\mathit{\theta })$ can be obtained in the following way. For the hth stratum, generate ${\mathit{n}}_h^{*}={(n_{h1}^{*},\dots ,n_{h{n}_h}^{*})}^{\mathrm{T}}$ by a multinomial distribution with $n_h-1$ trials and a success probability vector $n_h^{-1}{(1,\dots ,1)}^{\mathrm{T}}$ of length n_h. Then, the bootstrap weighted score function is ${\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })={\sum }_{h=1}^H{W}_h{\widehat{\mathit{S}}}_{w,h}^{*}(\mathit{\theta })$, where ${\widehat{\mathit{S}}}_{w,h}^{*}(\mathit{\theta })=n_h^{-1}{\sum }_{i=1}^{n_h}{r}_{hi}^{*}{\tilde{\mathit{S}}}_{w,hi}(\mathit{\theta })$, $r_{hi}^{*}=k_h{n}_{hi}^{*}$, and $k_h=n_h/(n_h-1)$; Rao, Wu, and Yue (1992) proposed a similar rescaling bootstrap method under stratified multi-stage sampling, but they did not investigate the corresponding theoretical properties with respect to hypothesis testing of super-population model parameters. By Lemma S7 in Section S4 of supplementary material, we conclude that $E_{*}\{{\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })\}={\widehat{\mathit{S}}}_w(\mathit{\theta })$ and ${\mathit{V}}_{*}\{{\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })\}=\widehat{V}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })|{\mathcal{F}}_N\}$, where $\widehat{V}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })|{\mathcal{F}}_N\}={\sum }_{h=1}^H{W}_h^2{\{n_h(n_h-1)\}}^{-1}{\tilde{\mathit{S}}}_{h,\text{vec}}(\theta )({\mathit{I}}_{n_h}-{\mathcal{P}}_{\mathbf{1},n_h}){\tilde{\mathit{S}}}_{h,\text{vec}}^{\mathrm{T}}(\theta ),{\mathcal{P}}_{\mathbf{1},n_h}={\mathbf{1}}_{n_h}{({\mathbf{1}}_{n_h}^{\mathrm{T}}{\mathbf{1}}_{n_h})}^{-1}{\mathbf{1}}_{n_h}^{\mathrm{T}}$ is the projection matrix to the linear space spanned by ${\mathbf{1}}_{n_h}$, ${\mathbf{1}}_{n_h}={(1,\dots ,1)}^{\mathrm{T}}$ is a vector consisting of n_h 1’s, ${\mathit{I}}_{n_h}$ is an $n_h\times n_h$ identity matrix, and $${\tilde{\mathit{S}}}_{h,\text{vec}}(\theta )=(d_{a(h1)}{p}_{a(h1)}^{-1}{\widehat{\mathit{S}}}_{w,a(h1)}(\mathit{\theta }),\dots ,d_{a(h{n}_h)}{p}_{a(h{n}_h)}^{-1}{\widehat{\mathit{S}}}_{w,a(h{n}_h)}(\mathit{\theta })).$$

In Section S4 of supplementary material, we show that the bootstrap central limit theorem in Theorem 2 holds with $n={\sum }_{h=1}^H{n}_h$.

Under stratified simple random sampling without replacement (Bickel and Freedman 1984), a bootstrap method is proposed in Section S7 of the supplementary material, and it is essentially the same as the one considered by Rao and Wu (1988) when estimating the population mean; see Rao and Wu (1988) and Beaumont and Patak (2012) for details.

To establish a bootstrap version of Wilks’ Theorem in (6), we use $$W^{*}(\widehat{\mathit{\theta }})=-2n\{l_w^{*}(\widehat{\mathit{\theta }})-l_w^{*}({\widehat{\mathit{\theta }}}^{*})\}$$

as the bootstrap version of $W({\mathit{\theta }}_0)$ in (5). Under the assumptions for Theorem 2, we can show that (9) $$W^{*}(\widehat{\mathit{\theta }})|A \to \mathcal{G}$$(9) in distribution as $n\to \infty$, and the reference distribution is the bootstrap sampling distribution conditional on the sample A, where $\mathcal{G}$ is the limiting distribution (6) of the quasi-likelihood ratio test statistic $W({\mathit{\theta }}_0)$ in (5) under the null hypothesis.

By (9), we conclude that the corresponding test statistic generated by the proposed bootstrap method has the same asymptotic distribution as in (5). Thus, we can use the bootstrap distribution of $W^{*}(\widehat{\mathit{\theta }})$ to approximate the sampling distribution of $W({\mathit{\theta }}_0)$.

Remark 2. A

s mentioned in Section 1, most existing bootstrap methods were proposed for variance estimation under survey sampling. There indeed exist some papers discussing bootstrap confidence intervals under survey sampling, but these methods do not apply to general hypothesis testing problems. For example, Rao, Wu, and Yue (1992) proposed a bootstrap-t confidence interval for the parameter of interest, but their bootstrap confidence interval does not apply when the parameter of interest is multi-dimensional. In addition, their bootstrap method is only valid under the stratified simple random sampling. Beaumont and Patak (2012) and Bertail and Combris (1997) discussed bootstrap confidence intervals in their simulation studies, but it is not clear how the corresponding intervals are constructed. Also, neither Bertail and Combris (1997) nor Beaumont and Patak (2012) investigated their bootstrap confidence intervals theoretically. Theorem 2, on the other hand, is a building block to show that the proposed bootstrap method can be used for hypothesis testing under general sampling designs. Although we adopt almost the same technique to prove Theorems 1–2, we are the first to establish that the bootstrap estimator ${\widehat{\mathit{\theta }}}^{*}$ has the same asymptotic distribution as the original $\widehat{\mathit{\theta }}$. Then, we can show that the corresponding bootstrap statistic $W^{*}(\widehat{\mathit{\theta }})$ shares the same asymptotic distribution as the original one $W({\mathit{\theta }}_0)$ under the null hypothesis, so the proposed bootstrap method can be used for hypothesis testing without deriving the corresponding variance estimator as in Lumley and Scott (2014). In addition, the proposed bootstrap method can be widely applied in commonly used sampling designs; see Sections S4–S5 of the supplementary material for details.

Remark 3.

In many practical situations, the design weights $w_i={\pi }_i^{-1}$ are further modified to incorporate population-level constraints (as in calibration weighting) or to handle unit nonresponse (as in nonresponse weighting adjustment). In this case, the bootstrap design weights can be modified in the same way to obtain the bootstrap final weights. See Bessonneau et al. (2021) for a detailed illustration of constructing bootstrap final weights in household surveys. When the test statistic is computed using the final weights rather than the design weights, the bootstrap final weights can be applied to the same test statistics to obtain the bootstrap distribution of the test statistic.

4 Quasi-Likelihood-Ratio Test

Without loss of generality, we denote $\mathit{\theta }$ as the true parameter ${\mathit{\theta }}_0$ in this section. Let ${\Theta }_0\subset \Theta$ be the parameter space under the null hypothesis, so the null hypothesis can be written as $H_0:\mathit{\theta }\in {\Theta }_0$. In this section, we consider $H_0:{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}$ for some known vector ${\mathit{\theta }}_2^{(0)}$, where $\mathit{\theta }={({\mathit{\theta }}_1^{\mathrm{T}},{\mathit{\theta }}_2^{\mathrm{T}})}^{\mathrm{T}}$. Thus, we have ${\Theta }_0=\{\mathit{\theta }\in \Theta ;{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}\}.$ Let ${\widehat{\mathit{\theta }}}_1^{(0)}$ be the profile pseudo M-estimator of ${\mathit{\theta }}_1$ under $H_0:{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}$, which is obtained by maximizing $l_w({\mathit{\theta }}_1,{\mathit{\theta }}_2^{(0)})$ with respect to ${\mathit{\theta }}_1$.

The quasi-likelihood ratio test statistic for testing $H_0:{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}$ is defined as (10) $$W({\mathit{\theta }}_2^{(0)})=-2n\{l_w({\widehat{\mathit{\theta }}}^{(0)})-l_w(\widehat{\mathit{\theta }})\},$$(10) where ${\widehat{\mathit{\theta }}}^{(0)}={({\widehat{\mathit{\theta }}}_1^{(0)}{}^{\mathrm{T}},{\mathit{\theta }}_2^{(0)}{}^{\mathrm{T}})}^{\mathrm{T}}$. Under simple random sampling with replacement, $W({\mathit{\theta }}_2^{(0)})$ in (10) is asymptotically distributed as ${\chi }^2(q)$ with $q=p-p_0$ and $p_0=\text{dim}({\Theta }_0)$, and the reference distribution is the joint distribution of the super-population model and the sampling mechanism, where $\text{dim}(\Theta )$ is the dimension of Θ. The following lemma, proved by Lumley and Scott (2014), presents the limiting distribution of the quasi-likelihood ratio test statistic in (10) for general sampling designs.

Lemma 1.

Under $H_0:{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}$ and the assumptions for Theorem 2, (11) $$W({\mathit{\theta }}_2^{(0)}) \to {\mathcal{G}}_1=\sum _{i=1}^q{c}_i{Z}_i^2$$(11) in distribution as $n\to \infty$, and the reference distribution is the joint distribution of the super-population model and the sampling mechanism, where $c_1\ge c_2\ge \cdots \ge c_q>0$ are the eigenvalues of $\mathit{P}={\mathbf{\Sigma }}_{\mathit{\theta },2}{\mathbf{ℐ}}_{22·1}({\mathit{\theta }}_1,{\mathit{\theta }}_2^{(0)})$, $Z_1,\dots ,Z_q$ are q independent random variables from the standard normal distribution, ${\mathbf{\Sigma }}_{\mathit{\theta },2}$ is the asymptotic variance of $n^{1/2}{\widehat{\mathit{\theta }}}_2$ in (3), ${\mathbf{ℐ}}_{22·1}({\mathit{\theta }}_1,{\mathit{\theta }}_2)={\mathbf{ℐ}}_{22}({\mathit{\theta }}_1,{\mathit{\theta }}_2)-{\mathbf{ℐ}}_{21}({\mathit{\theta }}_1,{\mathit{\theta }}_2){\{{\mathbf{ℐ}}_{11}({\mathit{\theta }}_1,{\mathit{\theta }}_2)\}}^{-1}{\mathbf{ℐ}}_{12}({\mathit{\theta }}_1,{\mathit{\theta }}_2),$ and ${\mathbf{ℐ}}_{ij}({\mathit{\theta }}_1,{\mathit{\theta }}_2)=E[-{\partial }^{2}l\{({\mathit{\theta }}_1,{\mathit{\theta }}_2);\mathit{X},Y\}/(\partial {\mathit{\theta }}_i\partial {\mathit{\theta }}_j^{\mathrm{T}})]$ for $i,j=1,2$.

Lumley and Scott (2014) proposed to estimate the limiting distribution in (11) using a design-based estimator of $\mathit{P}={\mathbf{\Sigma }}_{\mathit{\theta },2}{\mathbf{ℐ}}_{22·1}({\mathit{\theta }}_1,{\mathit{\theta }}_2^{(0)})$.

We consider an alternative test using a novel application of the bootstrap method discussed in Section 3. To do this, a bootstrap version of the quasi-likelihood ratio test statistic in (10) is obtained as (12) $$W^{*}({\widehat{\mathit{\theta }}}_2)=-2n\{l_w^{*}({\widehat{\mathit{\theta }}}_1^{*(0)},{\widehat{\mathit{\theta }}}_2)-l_w^{*}({\widehat{\mathit{\theta }}}^{*})\},$$(12) where ${\widehat{\mathit{\theta }}}_1^{*(0)}=\arg {\max }_{{\mathit{\theta }}_1}{l}_w^{*}({\mathit{\theta }}_1,{\widehat{\mathit{\theta }}}_2)$, and ${\widehat{\mathit{\theta }}}_2$ is the second component of $\widehat{\mathit{\theta }}={({\widehat{\mathit{\theta }}}_1^{\mathrm{T}},{\widehat{\mathit{\theta }}}_2^{\mathrm{T}})}^{\mathrm{T}}$. Under the conditions of Theorem 2, we can show that the limiting distribution for $W^{*}({\widehat{\mathit{\theta }}}_2)$ is the same as that in (11). Thus, the bootstrap distribution using $W^{*}({\widehat{\mathit{\theta }}}_2)$ in (12) can be used to obtain the reference distribution of the test statistic $W({\mathit{\theta }}_2^{(0)})$ in (10) under $H_0:{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}$.

Remark 4.

In addition to the quasi-likelihood ratio test, the proposed bootstrap method is also applicable to quasi-score tests (Rao, Scott, and Skinner 1998). To develop a quasi-score test for $H_0:{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}$, we first decompose the survey-weighted score function as (13) $${\widehat{\mathit{S}}}_w(\mathit{\theta })=\left(\begin{array}{c}{\widehat{\mathit{S}}}_{w1}(\mathit{\theta })\\ {\widehat{\mathit{S}}}_{w2}(\mathit{\theta })\end{array}\right)=\left(\begin{array}{c}{N}^{-1}\sum _{i\in A}{w}_i{\mathit{u}}_1(\mathit{\theta };{\mathit{x}}_i,y_i)\hfill \\ N^{-1}\sum _{i\in A}{w}_i{\mathit{u}}_2(\mathit{\theta };{\mathit{x}}_i,y_i)\hfill \end{array}\right)$$(13) with ${\mathit{u}}_j(\mathit{\theta };{\mathit{x}}_i,y_i)=(y_i-{\mu }_i){\{V_0({\mu }_i)\}}^{-1}(\partial {\mu }_i/\partial {\mathit{\theta }}_j)$ for j = 1, 2, where ${\mu }_i=\mu ({\mathit{x}}_i;\mathit{\theta }),\mu (\mathit{x};\mathit{\theta })=E(Y|\mathit{x})$ is a known function with an unknown parameter $\mathit{\theta }$, and $V_0({\mu }_i)=V(Y|\mathit{x})$ is a known function of μ_i. Note that quasi-score tests clearly specify the model without requiring a correctly specified density function, so they are more important in practice. Let ${\widehat{\mathit{\theta }}}_1^{(0)}$ be the solution to ${\widehat{\mathit{S}}}_{w1}({\mathit{\theta }}_1,{\mathit{\theta }}_2^{(0)})=\mathbf{0},$ where ${\widehat{\mathit{S}}}_{w1}(\mathit{\theta })$ is defined in (13). The quasi-score test statistic for $H_0:{\mathit{\theta }}_2={\mathit{\theta }}_2^{(0)}$ is (14) $$X_{QS}^2\left({\mathit{\theta }}_2^{(0)}\right)={\widehat{\mathit{S}}}_{w2}^{\mathrm{T}}\left({\widehat{\mathit{\theta }}}^{(0)}\right){\left\{{\widehat{\mathit{I}}}_{w22·1}\left({\widehat{\mathit{\theta }}}^{(0)}\right)\right\}}^{-1}{\widehat{\mathit{S}}}_{w2}\left({\widehat{\mathit{\theta }}}^{(0)}\right),$$(14) where (15) $${\widehat{\mathit{I}}}_{w22·1}(\mathit{\theta })={\widehat{\mathit{I}}}_{w22}(\mathit{\theta })-{\widehat{\mathit{I}}}_{w21}(\mathit{\theta }){\{{\widehat{\mathit{I}}}_{w11}(\mathit{\theta })\}}^{-1}{\widehat{\mathit{I}}}_{w12}(\mathit{\theta })$$(15)

and $$\left(\begin{array}{cc}{\widehat{\mathit{I}}}_{w11}(\mathit{\theta })\hfill & {\widehat{\mathit{I}}}_{w12}(\mathit{\theta })\hfill \\ {\widehat{\mathit{I}}}_{w21}(\mathit{\theta })\hfill & {\widehat{\mathit{I}}}_{w22}(\mathit{\theta })\hfill \end{array}\right)=-\left(\begin{array}{cc}{\widehat{\mathit{S}}}_{w1}(\mathit{\theta })/\partial {\mathit{\theta }}^{\prime }{}_1\hfill & {\widehat{\mathit{S}}}_{w1}(\mathit{\theta })/\partial {\mathit{\theta }}^{\prime }{}_2\hfill \\ {\widehat{\mathit{S}}}_{w2}(\mathit{\theta })/\partial {\mathit{\theta }}^{\prime }{}_1\hfill & {\widehat{\mathit{S}}}_{w2}(\mathit{\theta })/\partial {\mathit{\theta }}^{\prime }{}_2\hfill \end{array}\right).$$

The bootstrap replicates of $X_{QS}^2({\mathit{\theta }}_2^{(0)})$ can be constructed by replacing the sampling weights w_i and ${\widehat{\mathit{\theta }}}^{(0)}$, respectively, by the bootstrap weights $w_i^{*}$ and the bootstrap estimator that solves ${\widehat{\mathit{S}}}_w^{*}(\mathit{\theta })=N^{-1}\sum _{i\in A}{w}_i^{*}\{y_i-p_i(\mathit{\theta })\}{(1,{\mathit{x}}_i^{\mathrm{T}})}^{\mathrm{T}}=\mathbf{0}$ subject to ${\mathit{\theta }}_2={\widehat{\mathit{\theta }}}_2$. Then, the bootstrap version of the quasi-score test statistic $X_{QS}^2({\mathit{\theta }}_2^{(0)})$ is (16) $$X_{QS}^{2*}({\widehat{\mathit{\theta }}}_2)={\widehat{\mathit{S}}}_{w2}^{*\mathrm{T}}\left({\widehat{\mathit{\theta }}}^{*(0)}\right){\left\{{\widehat{\mathit{I}}}_{w22·1}^{*}\left({\widehat{\mathit{\theta }}}^{*(0)}\right)\right\}}^{-1}{\widehat{\mathit{S}}}_{w2}^{*}\left({\widehat{\mathit{\theta }}}^{*(0)}\right),$$(16) where ${\widehat{\mathit{\theta }}}^{*(0)}={({\widehat{\mathit{\theta }}}_1^{*(0)}{}^{\mathrm{T}},{\widehat{\mathit{\theta }}}_2{}^{\mathrm{T}})}^{\mathrm{T}}$, and ${\widehat{\mathit{I}}}_{w22·1}^{*}({\widehat{\mathit{\theta }}}^{*(0)})$ is computed similarly to (15) using bootstrap weights.

5 Test for Categorical Survey Data

5.1 Goodness-of-Fit Test

Suppose that a finite population U is partitioned into K categories with $U=U^{(1)}\cup \cdots \cup U^{(K)}$. Define $p_k=N_k/N$ to be the population proportion of the kth category, where N_k is the size of $U^{(k)}$. Thus, we implicitly assume that the finite population is an iid realization of the super-population model following a multinomial distribution with K categories.

Suppose that we are interested in the simple goodness-of-fit test $H_0:p_k=p_k^{(0)}$ for $k=1,\dots ,K$, where ${(p_1^{(0)},\dots ,p_K^{(0)})}^{\mathrm{T}}$ is a pre-specified vector satisfying ${\sum }_{k=1}^K{p}_k^{(0)}=1$. From the sample A, we compute ${\widehat{p}}_k={\widehat{N}}_k/\widehat{N}$ as an estimator of p_k, where ${\widehat{N}}_k=\sum _{i\in A_k}{w}_i$ is a design-unbiased estimator of N_k, $A_k=A\cap U^{(k)}$, and $\widehat{N}={\sum }_{k=1}^K{\widehat{N}}_k=\sum _{i\in A}{w}_i$. The estimated proportions can be obtained by $\widehat{\mathit{p}}=\arg {\max }_{\mathit{p}\in \Theta }\sum _{i\in A}{w}_{i}l(\mathit{p};{\mathit{x}}_i)$, where $l(\mathit{p};{\mathit{x}}_i)={\sum }_{k=1}^K{x}_{ik}\hspace{0.17em}\log \hspace{0.17em}(p_k)$, with x_ik = 1 if the ith element belongs to the kth category and 0 otherwise, $\mathit{p}={(p_1,\dots ,p_K)}^{\mathrm{T}}$.

Then, the Pearson Chi-squared goodness-of-fit test statistic for H₀ is $$X^2({\mathit{p}}^{(0)})=n\sum _{k=1}^K{\left({\widehat{p}}_k-p_k^{(0)}\right)}^2/p_k^{(0)},$$

where ${\mathit{p}}^{(0)}={(p_1^{(0)},\dots ,p_{K-1}^{(0)})}^{\mathrm{T}}$. Under regularity conditions, according to Rao and Scott (1981), we have (17) $$X^2({\mathit{p}}^{(0)}) \to {\mathcal{G}}_2=\sum _{k=1}^{K-1}{\lambda }_k{Z}_k^2$$(17) in distribution as $n\to \infty$ under H₀, and the reference distribution is the joint distribution of the super-population model and the sampling mechanism, where ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_{K-1}>0$ are the eigenvalues of the design effect matrix $\mathit{D}={\mathit{P}}_0^{-1}{\mathbf{\Sigma }}_p,{\mathit{P}}_0=\text{diag}({\mathit{p}}^{(0)})-{\mathit{p}}^{(0)}{({\mathit{p}}^{(0)})}^{\mathrm{T}},{\mathbf{\Sigma }}_p$ is the asymptotic variance of $\sqrt{n}\widehat{\mathit{p}}$, and $Z_1,\dots ,Z_{K-1}$ are K – 1 independent random variables from the standard normal distribution. Under simple random sampling with replacement, the limiting distribution in (17) is a Chi-squared distribution with K – 1 degrees of freedom since ${\mathit{P}}_0={\mathbf{\Sigma }}_p$.

We now apply the proposed bootstrap method to approximate the limiting distribution in (17), without computing the estimated covariance matrix ${\widehat{\mathbf{\Sigma }}}_p$. To describe the proposed method, let ${\widehat{\mathit{p}}}^{*}$ be the estimator of the population proportion $\mathit{p}$ based on the bootstrap weights $w_i^{*}$. The proposed bootstrap statistic of $X^2({\mathit{p}}^{(0)})$ is $$X^{2*}(\widehat{\mathit{p}})=n\sum _{i=1}^K{({\widehat{p}}_i^{*}-{\widehat{p}}_i)}^2/{\widehat{p}}_i.$$

We use ${\widehat{p}}_i$ in place of $p_i^{(0)}$ in the bootstrap test statistics. Under the conditions of Theorem 2, we can show that the limiting distribution for $X^{2*}(\widehat{\mathit{p}})$ is the same as that in (17).

5.2 Test of Independence in a Two-Way Table

The proposed bootstrap can also be used to test independence in a two-way count table. Let $p_{ij}=N_{ij}/N$ be the population proportion for cell (i, j) with margins $p_{i+}$ and $p_{+j}$ for $i=1,\dots ,R$ and $j=1,\dots ,C$, where R and C are the numbers of rows and columns, and $\{N_{ij};i=1,\dots ,R,j=1,\dots ,C\}$ is the set of population counts with margins $N_{i+}$ and $N_{+j}$. Let ${\widehat{N}}_{ij}$ be a design-unbiased estimator of N_ij and ${\widehat{p}}_{ij}={\widehat{N}}_{ij}/\widehat{N}$. By assuming a multinomial distribution for the super-population model, the Chi-squared statistic for testing independence $H_0:p_{ij}=p_{i+}{p}_{+j}$ for all i and j is (18) $$X_I^2=n\sum _{i=1}^R\sum _{j=1}^C\frac{{({\widehat{p}}_{ij}-{\widehat{p}}_{i+}{\widehat{p}}_{+j})}^2}{{\widehat{p}}_{i+}{\widehat{p}}_{+j}}.$$(18)

The limiting distribution of $X_I^2$ is a mixture of ${\chi }^2(1)$ distributions (Rao and Scott 1981).

To develop bootstrap replicates of the test statistics under H₀, let ${\widehat{p}}_{ij}^{*}$ be the bootstrap cell proportion computed using the bootstrap weights, ${\widehat{p}}_{i+}^{*}={\sum }_{j=1}^C{\widehat{p}}_{ij}^{*}$, and ${\widehat{p}}_{+j}^{*}={\sum }_{i=1}^R{\widehat{p}}_{ij}^{*}$. The proposed bootstrap version of $X_I^2$ is (19) $$X_I^{2*}=n\sum _{i=1}^R\sum _{j=1}^C\frac{{\{({\widehat{p}}_{ij}^{*}-{\widehat{p}}_{i+}^{*}{\widehat{p}}_{+j}^{*})-({\widehat{p}}_{ij}-{\widehat{p}}_{i+}{\widehat{p}}_{+j})\}}^2}{({\widehat{p}}_{i+}{\widehat{p}}_{+j})}.$$(19)

Under H₀, the terms in the numerator of $X_I^2$ are identical to ${\{({\widehat{p}}_{ij}-{\widehat{p}}_{i+}{\widehat{p}}_{+j})-(p_{ij}-p_{i+}{p}_{+j})\}}^2$. That is, the bootstrap test statistic is calculated by simply replacing $\{{\widehat{p}}_{ij},{\widehat{p}}_{i+},{\widehat{p}}_{+j}\}$ and $\{p_{ij},p_{i+},p_{+j}\}$ with $\{{\widehat{p}}_{ij}^{*},{\widehat{p}}_{i+}^{*},{\widehat{p}}_{+j}^{*}\}$ and $\{{\widehat{p}}_{ij},{\widehat{p}}_{i+},{\widehat{p}}_{+j}\}$, respectively. It can be shown that the limiting distribution of $X_I^{2*}$ conditional on the sample is the same as that of $X_I^2$.

6 Simulation Study

6.1 Single-Stage Sampling

In this section, we test the performance of the proposed bootstrap method under the PPS (probability proportional to size) sampling with replacement. A finite population of size N is generated as follows. For $i=1,\dots ,N$, $x_i\sim \mathrm{U}(-5,5)$ and $y_i={\theta }_1+{\theta }_2{x}_i+{ϵ}_i$, where $\mathrm{U}(a,b)$ is a uniform distribution on the interval (a, b), $({\theta }_1,{\theta }_2)=(1,1)$, and ${ϵ}_i\sim N(0,2^2)$. To make the sampling design informative, the PPS sampling with replacement is carried out to generate a sample of size n with selection probabilities $p_i\propto (6I_i+1)$ and ${\sum }_{i=1}^N{p}_i=1$, where I_i = 1 if $x_i{ϵ}_i>0$ and 0 otherwise; see Section S6 of the supplementary material for details about the probability proportional to size sampling and the bootstrap weights. The sampling design is informative (Pfeffermann 1993), since the selection probabilities are partially determined by the residuals $\{{ϵ}_i:i=1,\dots ,N\}$. A similar model setup was also used by Verret, Rao, and Hidiroglou (2015) in the context of small area estimation. We consider two scenarios: $(N,n)=(2,000,200)$ and $(N,n)=(15,000,500)$.

We are interested in testing $H_0:{\theta }_2={\theta }_2^{(0)}$ with the nominal significance level of $\alpha =0.05$ and consider three different values for ${\theta }_2^{(0)}\in \{1,1.1,1.2\}$. The following testing methods are compared:

1. Naive likelihood ratio method with $W({\theta }_2^{(0)})$ in (10) and ${\chi }^2(1)$ as test statistic and reference distribution, respectively.

2. Naive quasi-score method with test statistic (14) and ${\chi }^2(1)$ as the reference distribution.

3. Lumley and Scott (2014) method. The test statistic is $W({\theta }_2^{(0)})/\widehat{\delta }$, where $\widehat{\delta }=n\widehat{V}({\widehat{\theta }}_2){\widehat{I}}_{w,22·1},{\widehat{I}}_{w,22·1}$ is defined in (15), and $\widehat{V}({\widehat{\theta }}_2)$ is obtained from ${\widehat{\mathit{I}}}_w^{-1}(\widehat{\mathit{\theta }})\widehat{V}\{{\widehat{\mathit{S}}}_w(\mathit{\theta })\}{\{{\widehat{\mathit{I}}}_w^{-1}(\widehat{\mathit{\theta }})\}}^{\mathrm{T}}$. The reference distribution is $F_{1,k}$, where $F_{{\nu }_1,{\nu }_2}$ is an F distribution with parameters ν₁ and ν₂, k is the degrees of freedom of the variance estimator based on the sampling design, and k is obtained by subtracting the number of parameters from the effective sample size associated with the sampling design.

4. Bootstrap quasi-likelihood-ratio method with $W({\theta }_2^{(0)})$ in (10) being the test statistic, and the reference distribution is approximated by the empirical distribution of $W^{*}({\widehat{\theta }}_2)$ in (12).

5. Bootstrap quasi-score method with $X_{QS}^2({\theta }_2^{(0)})$ in (14) being the test statistic, and the reference distribution is approximated by the empirical distribution of the bootstrap test statistics,

When calculating the likelihood-ratio test statistics, we assume a normal distribution with a constant but unknown variance for the conditional distribution y_i given x_i. For each scenario, we generate 1000 Monte Carlo samples, and we consider M = 500 and M = 1000 iterations for both bootstrap methods. Table 1 summarizes the simulation results. For both scenarios, the naive likelihood-ratio method and the naive quasi-score method have a significantly inflated Type I error rate when the null hypothesis is true. The corresponding design effects (Wu and Thompson 2020) are more than 2 for estimating θ₂ under both scenarios; see Section S8.1 of the supplementary material for details. Thus, the distributions of $W({\theta }_2^{(0)})$ and the quasi-score test statistic (14) are more “dispersed” than ${\chi }^2(1)$ under the null hypothesis, and due to the lack of design-effect corrections, the two naive methods are not theoretically valid. The Type I error rates of the Lumley and Scott (2014) method and the two proposed bootstrap testing methods are similar under H₀, regardless of the population and sample sizes. The Type I error rates of the three methods come closer to the nominal truth as the population and sample sizes increase. Furthermore, the powers of the two bootstrap testing methods are reasonable compared to the Lumley-Scott method. The naive methods have a slightly larger power at the expense of elevated Type I error rates. In addition, we get similar test results regardless of the number of bootstrap repetitions.

Table 1 Test power for the hypothesis test $H_0:{\theta }_2={\theta }_2^{(0)}$ based on 1 000 Monte Carlo simulations, the true parameter is ${\theta }_2=1.00$, and the nominal significance level is 0.05.

(N, n)	${\theta }_2^{(0)}$	NLR	NQS	LS	M = 500		M = 1000
					BLR	BQS	BLR	BQS
$(2 000,200)$	1.00	0.20	0.20	0.07	0.07	0.07	0.07	0.07
	1.10	0.53	0.53	0.25	0.24	0.24	0.23	0.23
	1.20	0.89	0.89	0.73	0.72	0.72	0.72	0.72
$(15 000,500)$	1.00	0.22	0.23	0.06	0.06	0.06	0.06	0.06
	1.10	0.76	0.76	0.52	0.53	0.53	0.53	0.53
	1.20	1.00	1.00	0.98	0.98	0.98	0.98	0.98

NOTE: NLR, naive likelihood-ratio method; NQS, naive quasi-score method; LS, Lumley and Scott (2014) method; BLR, bootstrap quasi-likelihood-ratio method; BQS: bootstrap quasi-score method. The numbers of bootstrap repetitions are set to be M = 500 and M = 1000.

Besides, Section S9 of the supplementary material reports results of an additional simulation study for a stratified two-stage sampling design.

6.2 Test of Independence

In this section, we consider test of independence in a 3 × 3 table of counts to check the performance of the proposed bootstrap testing method. For a finite population $U=\{{\mathit{y}}_i:i=1,\dots ,N\}$, ${\mathit{y}}_i$ is generated by a multinomial distribution with one trial and a success probability vector $\mathit{p}$, where $\mathit{p}={(p_{11},\dots ,p_{ij},\dots ,p_{33})}^{\mathrm{T}}$, and p_ij is the success probability for the cell in the ith row and jth column for $i,j=1,\dots ,3$. That is, ${\mathit{y}}_i$ is a dummy variable consisting of eight 0’s and a 1. For the success probability vector $\mathit{p}$, we consider three cases:

• Case I: $p_{11}=1/4,p_{12}=p_{13}=p_{21}=p_{31}=1/8,p_{22}=p_{23}=p_{32}=p_{33}=1/16$.

• Case II: $p_{11}=1/4,p_{12}=p_{13}=(1.4)/8,p_{21}=p_{31}=(0.6)/8,p_{22}=p_{33}=1/16,p_{23}=(1.4)/16,p_{32}=(0.6)/16$.

• Case III: $p_{11}=p_{2,3}=p_{3,2}=1/6,p_{12}=p_{13}=p_{21}=p_{31}=p_{22}=p_{33}=1/12$.

Case I satisfies independence for the two-way table of counts, but Cases II–III do not. The level of nonindependence can be expressed by $\gamma ={\sum }_{i=1}^3{\sum }_{j=1}^3{(p_{ij}-p_{i+}{p}_{+j})}^2/(p_{i+}{p}_{+j})$, and γ = 0 corresponds to independence. The values of γ are 0, 0.017 and 0.125 for Cases I–III, respectively.

For each ${\mathit{y}}_i$, we generate an auxiliary variable $x_i={\mathit{\beta }}^{\mathrm{T}}{\mathit{y}}_i$, where $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_9)}^{\mathrm{T}},{\beta }_j=0.5+e_j$ for $j=1,\dots ,9,e_j\sim \text{Ex}(1)$. Probability proportional to size sampling with replacement is used to generate a sample of size n with selection probability proportional to x_i. We consider two scenarios for the population and sample sizes: $(N,n)=(2000,75)$ and $(N,n)=(10,000,150)$.

We are interested in testing independence in the two-way table with $\alpha =0.05$ nominal significance level. For each sample, we compare the following test methods:

1. Naive Pearson method based on $X_I^2$ in (18) with ${\chi }^2(4)$ being the reference distribution.

2. The second-order Rao-Scott correction method using the Satterthwaite approximation (Thomas and Rao 1987). This method is widely used in analyzing two-way contingency table and is implemented by the command svychisq in the R package survey (Lumley 2020).

3. Bootstrap Pearson method using $X_I^2$, and its distribution is approximated by that of $X_I^{2*}$ in (19).

For each scenario, we generate 1000 Monte Carlo samples, and we consider M = 500 and M = 1000 iterations for both bootstrap methods. Table 2 summarizes the simulation results. For Case I when the null hypothesis is true, the Type I error rate of the naive Pearson method is much larger than the nominal significance level 0.05 for different sample sizes, since the asymptotic distribution of the test statistic $X_I^2$ in (18) is more “dispersed” than a ${\chi }^2(4)$ distribution under the null hypothesis; see Section S8.2 of the supplementary material for details. The performance of the Rao-Scott method with second-order correction is similar to the bootstrap testing method in terms of Type I error and power. The power of the bootstrap testing method increases with the value of γ. Similar to the previous simulation, the power of the naive method is larger at the expense of elevated Type I error rates. In addition, we get similar test results regardless of the number of bootstrap repetitions.

Table 2 Power of the test procedures for independence based on 1000 Monte Carlo simulation samples, Case I corresponds to the independence scenario, and the nominal significance level is 0.05.

(N, n)	Case	NP	RS	BP_I	BP_II
$(2 000,75)$	I	0.16	0.06	0.07	0.06
	II	0.21	0.10	0.12	0.11
	III	0.71	0.53	0.56	0.55
$(10 000,150)$	I	0.14	0.05	0.07	0.06
	II	0.25	0.13	0.14	0.14
	III	0.94	0.84	0.84	0.84

NOTE: NP, naive Pearson method; NLR: naive likelihood-ratio method; RS, Rao-Scott method using second-order correction; BP_I, bootstrap Pearson method with M = 500; BP_II, bootstrap Pearson method with M = 1000. The numbers of bootstrap repetitions are set to be M = 500 and M = 1000.

7 Application

We present an analysis of the 2011 Private Education Expenditure Survey (PEES) in South Korea using the proposed bootstrap testing methods. The purpose of this survey is to study the relationship between private education expenditure and academic performance of students before entering college.

A stratified two-stage sampling design was used for the 2011 PEES, and the strata consist of 16 first-tier administrative divisions, including most provinces and metropolitan cities in South Korea. For each stratum, PPS sampling with replacement was conducted in the first stage, and the primary sampling unit was the school. The students were randomly selected from the sample school in the second stage. There are about 1000 sample schools and 45,000 students involved in this survey.

For the student i in the sample A, let y_i be the academic performance assessed by the teacher, and it takes a value from 1 through 3 corresponding to low, middle and high academic performance, respectively. Associated with y_i, let ${\mathit{x}}_i$ be the vector of covariates of interest. As discussed by Kim, Park, and Lee (2017), we consider the following covariates: after-school education, hours taking lessons provided by the school after regular classes in a month; father’s education, 1 for college or higher and 0 otherwise; gender, 1 for female and 0 for male; household income per month; mother’s education, 1 for college or higher and 0 otherwise; private education, hours taking private lessons in a month.

In this section, we study the academic performance of students in middle school and high school separately and are interested in studying the conditional probability of achieving high academic performance. Specifically, consider the following logistic model (20) $$\text{logit}\{\text{pr}(Y=1|\mathit{x})\}=(1,{\mathit{x}}^{\mathrm{T}})\mathit{\theta },$$(20) where $\text{logit}(x)=\hspace{0.17em}\log \hspace{0.17em}(x)-\hspace{0.17em}\log \hspace{0.17em}(1-x)$ for $x\in (0,1)$, Y = 1 if high academic performance is achieved and 0 otherwise, $\mathit{x}$ is a vector of six covariates, and $\mathit{\theta }={({\theta }_0,{\theta }_1,\dots ,{\theta }_6)}^{\mathrm{T}}$. We are interested in testing the null hypotheses $H_{0,i}:{\theta }_i=0$ for $i=1,\dots ,6$ with $\alpha =0.05$ nominal significance level. Since the Wald method is widely used in practice, the naive likelihood ratio method, naive quasi-score method, bootstrap quasi-likelihood ratio method, bootstrap quasi-score method with 1000 iterations and a two-sided Wald test are compared. The p-values for the two naive methods and the Wald test are obtained using reference distributions for simple random sampling. Specifically, the reference distribution for the two naive methods is ${\chi }^2(1)$, and that of the Wald test is a normal distribution with estimated variance (Fuller 2009, sec. 1.2.8) for $\widehat{\mathit{\theta }}$ using the sandwich formula. The p-values for the proposed bootstrap testing methods are obtained by bootstrap empirical distributions of the corresponding test statistics.

The results of our analysis are summarized in Table 3. The two bootstrap testing methods and the Wald test perform approximately the same in terms of p-values. The p-values of the two naive methods are approximately the same, but they differ from those of the bootstrap testing methods, especially for “after school education” at the middle school level and most covariates at the high school level, as the naive methods may not properly reflect the intra-cluster correlation due to cluster sampling.

Table 3 Estimates (Est) and the p-values (Unit: ${10}^{-2}$) for testing $H_{0,i}:{\theta }_i=0$, where $i=1,\dots ,6$, in (20) for the middle school and high school.

School level	Cov	Est	p-value (Unit: ${10}^{-2}$)
			NLR	NQS	BLR	BQS	Wald
Middle School	After-school Edu	0.03	16.9	16.6	39.4	39.0	35.6
	Father’s Edu	0.42	0.0	0.0	0.0	0.0	0.0
	Gender	0.3	0.3	0.3	0.7	0.6	0.6
	Income	0.14	0.0	0.0	0.2	0.2	0.4
	Mother’s Edu	0.37	0.2	0.2	0.0	0.0	0.0
	Private Edu	0.05	0.0	0.0	0.0	0.0	0.0
High School	After-school Edu	0.03	2.0	2.0	0.3	0.3	0.1
	Father’s Edu	0.19	10.4	10.3	14.1	14.0	13.8
	Gender	0.21	2.5	2.5	17.0	16.8	15.7
	Income	0.1	0.0	0.0	0.1	0.1	0.1
	Mother’s Edu	0.15	20.4	20.4	41.2	41.1	40.3
	Private Edu	0.04	0.0	0.0	3.5	3.0	3.2

NOTE: Cov, covariate; Est, estimate; NLR, naive likelihood-ratio method; NQS, naive quasi-score method; BLR, bootstrap quasi-likelihood-ratio method; BQS: bootstrap quasi-score method; Wald, two-sided Wald test.

Based on the two bootstrap testing methods in Table 3, we have the following conclusions under $\alpha =0.05$ nominal significance level. Controlling other covariates, the probability that female students achieve high academic performance is significantly higher than that of male students in middle school, but the gender effect is not significant in high school. Hours spent on private education and after-school education can significantly increase the probability of achieving high academic performance in both middle school and high school. The household income has a significantly positive influence on their child’s academic performance. However, the level of education of the father and mother has a significant influence during the middle school period only.

8 Concluding Remarks

Many statistical agencies provide microdata files containing survey weights and several sets of associated replication weights, in particular bootstrap weights. Standard statistical packages often permit the use of survey-weighted test statistics, and we have shown how to approximate their distributions under the null hypothesis by their bootstrap analogues computed from the bootstrap weights supplied in the data file. Using bootstrap weights, we can easily apply hypothesis testing without using specialized software designed for complex survey data, as noted by Beaumont and Bocci (2009). The proposed bootstrap method is applied to quasi-likelihood ratio tests and quasi-score tests.

Our theories depend on establishing bootstrap central limit theorems under specified sampling designs, including stratified multi-stage sampling with clusters selected with replacement. Under some sampling designs, the central limiting theorem does not necessarily hold and the proposed method is not applicable. In the case of sampling designs that allow for the central limit theorem, the bootstrap central limit theorems can be established under some additional assumptions. Under the assumption that the bootstrap central limit theorem holds in the sampling design, the proposed hypothesis testing methods are applicable. The hypotheses we have discussed are special cases of the one considered in Section 4.2 of Jiang (2019), so the theoretical properties of the proposed bootstrap method can also be investigated using the general method developed in Section 4.2 of Jiang (2019). The proposed method can be applied to the case of categorical data by developing bootstrap procedures for testing simple goodness of fit and independence in a two-way table. We plan to extend our bootstrap method for categorical data to test hypotheses from multi-way tables of weighted counts or proportions, using a log-linear model approach proposed by Rao and Scott (1984). Moreover, the bootstrap method can also be applied to nonresponse and weighting adjustments for variance estimation under regularity conditions, and extending the proposed bootstrap method for weight adjustments for hypothesis testing is also a promising project.

Supplementary Materials

Supplementary Material contains the following: comments on the regularity conditions (S1), proof of Theorem 1 (S2), proof of Theorem 2 (S3), verification of the proposed bootstrap method under stratified multi-stage sampling (S4), bootstrap method for Poisson sampling (S5) bootstrap method for probability proportional to size sampling with replacement (S6), stratified simple random sampling (S7), design effect (S8) and an additional simulation study: stratified two-stage sampling (S9).

Acknowledgments

We are grateful to the associate editor and three anonymous referees for their constructive comments, which have greatly improved the article.

Additional information

Funding

Kim is partially supported by NSF (OAC-1931380). Rao is supported by a grant from the Natural Sciences and Engineering Research Council of Canada. Wang is partially supported by NSFC (11901487, 71988101, 72033002, 12231011) and National Key R&D Program of China (2022YFA10038002).