2 Analysis of embedded K x L factorial experiments

Jan A. van den Brakel

Previous | Next

2.1 Experimental designs embedded in probability samples

In a $K\times L$ factorial design, the effects of two factors are tested simultaneously. The first factor, denoted A contains $K\ge 2$ levels. The second factor, denoted B contains $L\ge 2$ levels. The purpose of the experiment is to test the main effects of the two factors and the interactions between both factors on the main parameter estimates of the ongoing survey. To this end a probability sample $s$ of size $n$ is drawn from a finite target population U of size N according the sample design of the regular survey. This sample design can be generally complex, and is described by its first order inclusion probabilities ${\pi }_i$ for unit i and second order inclusion probabilities ${\pi }_{ii\text{'}}$ for units i and i'.

Subsequently, this sample is randomly divided into $KL$ subsamples according to a randomized experiment. In the case of a CRD, the sample $s$ of size $n$ is randomly divided into $KL$ subsamples $s_{kl}$, each with a size of $n_{kl}$ sampling units. The sampling units of each subsample are assigned to one of the $KL$ treatment combinations. Under a CRD, $n_{++}={\displaystyle {\sum }_{k=1}^K{\displaystyle {\sum }_{l=1}^L{n}_{kl}}}$ denotes the total number of sampling units in the sample $s$. The probability that sampling unit i is assigned to subsample $s_{kl}$, conditionally on the realization of $s$, equals $n_{kl}/n_{++}$. The unconditional probability that sampling unit i is selected in subsample $s_{kl}$ equals ${\pi }_i^{*}={\pi }_i(n_{kl}/n_{++})\text{.}$

The power of an experiment might be improved by using sampling structures such as strata, clusters or interviewers as block variables in an RBD since restricted randomization removes the variance between the blocks from the analysis of the experiment, (Fienberg and Tanur (1987, 1988)). In the case of an RBD, the sampling units are deterministically grouped in B more or less homogeneous blocks $s_b$. Within each block, the sampling units are randomly assigned to one of the $KL$ treatment combinations. Let $n_{bkl}$ denote the number of sampling units in block $b$ assigned to treatment combination $kl$, and $n_{b++}={\displaystyle {\sum }_{k=1}^K{\displaystyle {\sum }_{l=1}^L{n}_{bkl}}}$ the number of sampling units in block $b$. The probability that sampling unit i is assigned to subsample $s_{kl}$, conditionally on the realization of $s$ and $i\in s_b$, equals $n_{bkl}/n_{b++}$, $i\in s_b$. The unconditional probability that sampling unit i is selected in subsample $s_{kl}$ equals ${\pi }_i^{*}={\pi }_i(n_{bkl}/n_{b++})\text{.}$

In many practical applications one of the $KL$ subsamples is assigned to the regular survey and serves, besides being used to produce estimates for the regular publication, as the control group in the experiment. In such situations, the size of this subsample will be substantially larger than the other subsamples.

There are a lot of issues in the planning and design stage of embedded experiments. The field staff, for example, requires special attention, since an embedded experiment can have a large impact on their daily routine of data collection, to which they are accustomed. See van den Brakel and Renssen (1998) and van den Brakel (2008) for more details about such design issues.

Although factorial designs are efficient from a statistical point of view, there might be strong practical arguments against a factorial set-up. The number of treatment combinations increases rapidly with the number of factors in full factorial designs, which might be difficult to implement in the data collection of a survey process. A general solution, known from standard experimental design theory, is to confound higher order interactions with blocks or to apply fractional factorial designs (Hinkelmann and Kempthorne (2005); Montgomery (2001)). These balanced designs, however, are generally hard to combine with the fieldwork restrictions encountered in the daily practice of survey sampling. In many applications the factors that changed in a survey redesign are therefore combined into one treatment. The total effect of these modifications is tested against the standard alternative in a two-treatment experiment. This implies that the effects of all factors in the experiment are confounded and cannot be separately estimated.

2.2 Testing hypotheses about finite population parameters

The purpose of embedded experiments is to test whether alternative survey implementations result in significantly different estimates for finite population parameters. Such differences are the result of non-sampling errors, like measurement errors and response bias. A measurement error model is required to link systematic differences between finite population parameters due to different survey implementations or treatments. Therefore the measurement error model for single-factor experiments proposed by van den Brakel and Renssen (2005) and van den Brakel (2008) is extended to factorial designs.

Let $y_{iqkl}$ denote the observation obtained from the $i^{th}$ individual observed under the $k{l}^{th}$ treatment combination and the $q^{th}$ interviewer. It is assumed that the observations are a realization of the measurement error model

$$y_{iqkl}=u_i+{\beta }_{kl}+{\gamma }_q+{\epsilon }_{ikl}.\left(2.1\right)$$

Here $u_i$ is the true intrinsic value of the $i^{th}$ individual, ${\beta }_{kl}$ the effect of the $k{l}^{th}$ treatment combination and ${\epsilon }_{ikl}$ an error component. The model also allows for interviewer effects, i.e. ${\gamma }_q=\psi +{\xi }_q$, where $\psi$ denotes a systematic interviewer bias and ${\xi }_q$ the random effect of the $q^{th}$ interviewer, respectively. Let ${\text{E}}_m$ and ${\mathrm{cov}}_m$ denote the expectation and the covariance with respect to the measurement error model. It is assumed that ${\text{E}}_m({\epsilon }_{ikl})=0$, ${\mathrm{var}}_m({\epsilon }_{ikl})={\sigma }_{ikl}^2$, and that measurement errors between sampling units are independent. Furthermore it is assumed that ${\text{E}}_m({\xi }_q)=0$, ${\mathrm{var}}_m({\xi }_q)={\tau }_q^2$ and that random interviewer effects between interviewers are independent. As a result the model allows for correlated response between sampling units that are interviewed by the same interviewer. The measurement error model allows for separate variances for measurement errors under different treatment combinations and separate variances for interviewers.

The treatment effects ${\beta }_{kl}$ can be interpreted as the bias in the estimated population parameter if the true intrinsic population value of $u$ is measured by means of the $k{l}^{th}$ survey implementation. The treatment effect can be decomposed in the traditional way of an analysis of variance for a two-way layout:

$${\beta }_{kl}=u+A_k+B_l+A{B}_{kl},\left(2.2\right)$$

with $u$ the overall effect, $A_k$ and $B_l$ the main effects of treatment factors $A$ and $B$ and $A{B}_{kl}$ the interactions between treatment factors $A$ and $B$. If the treatment effects are defined as fixed deviations from the individuals' intrinsic value $u_i$, then the overall mean $u$ equals zero. In that case $A_k$ corresponds with the bias associated with the $k^{th}$ level of factor $A$ averaged over all levels of factor $B$, $B_l$ the bias associated with the $l^{th}$ level of factor $B$, averaged over all levels of factor $A$, and $A{B}_{kl}$ the additional bias associated with the combination of the $k^{th}$ level of factor $A$ and the $l^{th}$ level of factor $B$ on top of $A_k$ and $B_l\text{.}$

The following restrictions are required to identify model (2.2):

$${\displaystyle \sum _{k=1}^K{A}_k}=0,{\displaystyle \sum _{l=1}^L{B}_l}=0,\left(2.3\right)$$

and

$${\displaystyle \sum _{k=1}^{K}A{B}_{kl}}=0,l=1,2,\dots ,L,{\displaystyle \sum _{l=1}^{L}A{B}_{kl}}=0,k=1,2,\dots ,K.\left(2.4\right)$$

For each sampling unit, a potential response variable is defined under each of the $KL$ treatment combinations. Therefore the measurement error model can be expressed in matrix notation as:

$$y_{iq}=j_{KL}{u}_i+\beta +j_{KL}{\gamma }_q+{\epsilon }_i,\left(2.5\right)$$

where $y_{iq}={(y_{iq11},\mathrm{...},y_{iqkl},\mathrm{...},y_{iqKL})}^t$, $\beta ={({\beta }_{11},\mathrm{...},{\beta }_{kl},\mathrm{...},{\beta }_{KL})}^t$, $j_{KL}$ a vector of order $KL$ with each element equal to one and ${\epsilon }_i={({\epsilon }_{i11},\mathrm{...},{\epsilon }_{ikl},\mathrm{...},{\epsilon }_{iKL})}^t$. The sampling units are assigned to one of the treatment combinations only, so only one of the responses of $y_{iq}$ is actually observed. The model assumptions specified above are stated as:

$${\text{E}}_m\left({\epsilon }_i\right)=0,\left(2.6\right)$$

$${\text{cov}}_m\left({\epsilon }_i,{\epsilon }_{i^{\prime }}\right)=\{\begin{array}{lll}{\Sigma }_i\hfill & :\hfill & i=i^{\prime }\hfill \\ {\rm O}\hfill & :\hfill & i\ne i^{\prime }\hfill \end{array},\left(2.7\right)$$

$${\text{E}}_m\left({\xi }_q\right)=0,\left(2.8\right)$$

$${\text{cov}}_m\left({\xi }_q,{\xi }_{q^{\prime }}\right)=\{\begin{array}{lll}{\tau }_q^2\hfill & :\hfill & q=q^{\prime }\hfill \\ 0\hfill & :\hfill & q\ne q^{\prime }\hfill \end{array},\left(2.9\right)$$

$${\mathrm{cov}}_m\left({\epsilon }_{ikl},{\xi }_q\right)=0,\left(2.10\right)$$

where $0$ is a vector of order $KL$ with each element zero, ${\Sigma }_i$ a matrix of order $KL\times KL$ containing the variances of the measurement errors ${\sigma }_{ikl}^2$, and ${\rm O}$ a matrix of order $KL\times KL$ with each element zero.

Let $\overline{Y}={({\overline{Y}}_{11},\mathrm{...},{\overline{Y}}_{1L},\mathrm{...},{\overline{Y}}_{kl},\mathrm{...},{\overline{Y}}_{K1},\mathrm{...},{\overline{Y}}_{KL})}^t$ denote the $KL$ dimensional vector of population means of $y_{iq}$ defined by (2.5). These are the values obtained under a complete enumeration of the finite population under each of the treatment combinations and are defined as:

$$\overline{Y}=j_{KL}\frac{1}{N}{\displaystyle \sum _{i=1}^N{u}_i}+\beta +j_{KL}\psi +j_{KL}{\displaystyle \sum _{q=1}^Q\frac{N_q}{N}{\xi }_q}+\frac{1}{N}{\displaystyle \sum _{i=1}^N{\epsilon }_i},\left(2.11\right)$$

where $Q$ denotes the total number of interviewers available for the data collection and $N_q$ the number of units assigned to the $q-\text{th}$ interviewer in the case of a complete enumeration.

Only systematic differences between the population parameters that are reflected by the treatment effects $\beta$ should lead to a rejection of the null hypotheses of no treatment effects. This is accomplished by formulating hypotheses about $\overline{Y}$ in expectation over the measurement error model, i.e.

$${\text{E}}_m\overline{Y}=j_{KL}\frac{1}{N}{\displaystyle \sum _{i=1}^N{u}_i}+\beta +j_{KL}\psi .\left(2.12\right)$$

Consequently, hypotheses about main effects and interactions are formulated as

$$\begin{array}{l}{H}_0\text{:}C{\text{E}}_m\overline{Y}=0,\left(2.13\right)\\ H_1\text{:}C{\text{E}}_m\overline{Y}\ne 0,\end{array}$$

where $C$ denotes an appropriate contrast matrix, and $0$ a vector with elements equal to one and a dimension that is equal to the number of contrasts (rows) defined by $C$. The contrast matrix for the hypothesis about the main effects of factor $A$ is defined as

$$C_A=\frac{1}{L}\left(j_{(K-1)}|-I_{(K-1)}\right)\otimes j_L^t\equiv \frac{1}{L}{\tilde{C}}_A\otimes j_L^t\text{, }\left(2.14\right)$$

with $I_{(K-1)}$ the identity matrix of order $K-1$. Matrix ${\tilde{C}}_A$ defines the $K-1$ contrasts between the $K$ levels of factor $A$, averaged over the $L$ levels of factor $B$. From (2.12) and due to restrictions (2.3) and (2.4) it follows that the contrasts between the population parameters exactly correspond to the contrasts between the main effects of the first factor:

$${\tilde{C}}_A{\text{E}}_m\overline{Y}={\tilde{C}}_A\beta ={(A_1-A_2,\mathrm{...},A_1-A_K)}^t\text{.}$$

The contrast matrix for the hypothesis about the main effects of factor $B$ is defined as

$$C_B=\frac{1}{K}{j}_K^t\otimes \left(j_{(L-1)}|-I_{(L-1)}\right)\equiv \frac{1}{K}{j}_K^t\otimes {\tilde{C}}_B\text{.}\left(2.15\right)$$

This matrix defines the $L-1$ contrasts between the $L$ levels of factor $B$, averaged over the $K$ levels of factor $A$. From (2.12) and due to restrictions (2.3) and (2.4) it follows that the contrasts between the population parameters exactly correspond to the contrasts between the main effects of the second factor:

$${\tilde{C}}_B{\text{E}}_m\overline{Y}={\tilde{C}}_B\beta ={(B_1-B_2,\mathrm{...},B_1-B_L)}^t\text{.}$$

The contrast matrices for the main effects use the first level of factors $A$ and $B$ as the reference category. This implies that treatment combination $A_1\times B_1$ is considered as the control group in the experiment.

Interactions between the two treatment factors are defined as the $L-1$ contrasts of factor $B$ between the $K-1$ contrasts of factor $A$ or, equivalently, as the $K-1$ contrasts of factor $A$ between the $L-1$ contrasts of factor $B$, Hinkelmann and Kempthorne (1994, chapter 11). Therefore the contrast matrix for the hypothesis about the interactions between factor $A$ and $B$ can be defined as

$$C_{AB}=\left(j_{\left(K-1\right)}|-I_{\left(K-1\right)}\right)\otimes \left(j_{\left(L-1\right)}|-I_{\left(L-1\right)}\right)={\tilde{C}}_A\otimes {\tilde{C}}_B.\left(2.16\right)$$

This matrix contains the $(K-1)(L-1)$ contrasts that define the interactions between factor $A$ and $B$. The contrasts between the population parameters exactly correspond to the interactions between the first and the second factor, since

$${\tilde{C}}_{AB}{\text{E}}_m\overline{Y}={\tilde{C}}_{AB}\beta =(A{B}_{11}-A{B}_{12}-A{B}_{21}+A{B}_{22},\mathrm{...},$$

$$A{B}_{11}-A{B}_{1L}-A{B}_{21}+A{B}_{2L},\mathrm{...},A{B}_{11}-A{B}_{12}-A{B}_{K1}+A{B}_{K2},\mathrm{...},A{B}_{11}-A{B}_{1L}-A{B}_{K1}+A{B}_{KL}{)}^t$$

Each element of this $(K-1)(L-1)$ vector defines one of the $(K-1)(L-1)$ interactions, which neatly corresponds to the contrasts between the interaction effects defined by (2.2). The first element e.g. can be interpreted as the deviation of the treatment effect of the particular combination of factor $A$ at level 2 and factor $B$ at level 2 from the two main effects of these factors.

2.3 Wald test

The hypotheses specified in section 2.2, can be tested with a Wald test (Wald 1943), which is frequently applied in design-based testing procedures, see for example Skinner, Holt and Smith (1989) or Chambers and Skinner (2003). If $\widehat{\overline{Y}}$ denotes a design-unbiased estimator for $\overline{Y}$, $C$ the contrast matrix $C_A$, $C_B$, or $C_{AB}$ defined in (2.14), (2.15) and (2.16), and $\mathrm{cov}(C\widehat{\overline{Y}})$ the covariance matrix of the contrasts between $\widehat{\overline{Y}}$, then hypotheses can be tested with the Wald statistic $W={\widehat{\overline{Y}}}^t{C}^t{\{\mathrm{cov}(C\widehat{\overline{Y}})\}}^{-1}C\widehat{\overline{Y}}$. The GREG estimators, proposed by van den Brakel and Renssen (2005) and van den Brakel (2008) for single-factor experiments are extended to embedded factorial designs in this section. For notational convenience, the subscript q will be omitted in $y_{iqkl}$, since there is no need to sum explicitly over the interviewer subscript in most of the formulas developed in the rest of this paper.

To apply the model-assisted mode of inference to the analysis of embedded experiments, it is assumed for each unit in the population that the intrinsic value $u_i$ in measurement error model (2.5) is an independent realization of the following linear regression model:

$$u_i={\beta }^t{x}_i+e_i,\left(2.17\right)$$

where $x_i$ H-vector with auxiliary information, $\beta$ a H-vector with the regression coefficients and $e_i$ the residuals, which are independent random variables with variance ${\omega }_i^2$. It is required that all ${\omega }_i^2$ are known up to a common scale factor, that is ${\omega }_i^2={\omega }^2{\nu }_i$, with ${\nu }_i$ known. The GREG estimator for ${\overline{Y}}_{kl}$, based on the $n_{kl}$ observations of subsample $s_{kl}$, is defined as (Särndal et al., 1992)

$${\widehat{\overline{Y}}}_{kl;greg}={\widehat{\overline{Y}}}_{kl}+{\widehat{b}}_{kl}^t(\overline{X}-\widehat{\overline{X}})\text{, }k=1,2,\dots ,K, \text{and} l=1,2,\dots ,L,\left(2.18\right)$$

where,

$${\widehat{\overline{Y}}}_{kl}=\frac{1}{N}{\displaystyle \sum _{i=1}^{n_{kl}}\frac{y_{ikl}}{{\pi }_i^{*}}},\left(2.19\right)$$

denotes the HT estimator for ${\overline{Y}}_{kl}$, $\overline{X}$ the finite population means of the auxiliary variables $x$, and $\widehat{\overline{X}}$ the HT estimator for $\overline{X}$ based on the $n_{kl}$ sample units of subsample $s_{kl}$. Furthermore,

$${\widehat{b}}_{kl}={\left({\displaystyle \sum _{i=1}^{n_{kl}}\frac{x_i{x}_i^t}{{\omega }_i^2{\pi }_i^{*}}}\right)}^{-1}{\displaystyle \sum _{i=1}^{n_{kl}}\frac{x_i{y}_{ikl}}{{\omega }_i^2{\pi }_i^{*}}},\left(2.20\right)$$

denotes the HT-type estimator for the regression coefficients in (2.17) based on the $n_{kl}$ sampling units in subsample $s_{kl}$. In (2.19) and (2.20), ${\pi }_i^{*}$ are the first order inclusion probabilities for the sampling units in the $KL$ different subsamples, derived in subsection 2.1. Now ${\widehat{\overline{Y}}}_{GREG}={({\widehat{\overline{Y}}}_{11;greg},\mathrm{...},{\widehat{\overline{Y}}}_{KL;greg})}^t$ is an approximately design-unbiased estimator for $\overline{Y}$ and also for ${\text{E}}_m\overline{Y}$ by definition.

Under the null hypotheses that there are no treatment effects and no interactions, it follows that $b_{kl}=b_{k\text{'}l\text{'}}$. In that case, it might be efficient to substitute for ${\widehat{b}}_{kl}$ in the GREG estimator (2.18) the pooled estimator

$$\widehat{b}={\left({\displaystyle \sum _{i=1}^n\frac{x_i{x}_i^t}{{\omega }_i^2{\pi }_i^{*}}}\right)}^{-1}{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{i=1}^{n_{kl}}\frac{x_i{y}_{ikl}}{{\omega }_i^2{\pi }_i^{*}}}}}.\left(2.21\right)$$

Since $H$ instead of $KL\times H$ regression coefficients have to be estimated, the pooled estimates of the regression coefficients $\widehat{b}$ will be more precise, particularly in the case of small subsamples. Note, however, that many commonly used weighting schemes meet the condition that a constant vector $\lambda$ exists such that ${\omega }_i^2=\lambda x_i$ for all $i\in U$. In this situation the GREG estimator reduces to the simplified form ${\widehat{\overline{Y}}}_{kl;greg}={\widehat{b}}_{kl}^t\overline{X}$ (Särndal et al. 1992, section 6.5). Under this simplified form, the treatment effects are completely included in the regression coefficients. In case of the pooled estimator (2.21), the $KL$ GREG estimators are exactly equal by definition, since ${\widehat{\overline{Y}}}_{kl;greg}={\widehat{b}}_{}^t\overline{X}$ for all k and l.

An expression for the covariance matrix of the contrasts between the elements of ${\widehat{\overline{Y}}}_{GREG}$ where the covariance is taken over the sampling design, the experimental design and the measurement error model, is given by

$$\mathrm{cov}\left(C{\widehat{\overline{Y}}}_{GREG}\right)={\text{E}}_m{\text{E}}_{s}CD{C}^t,\left(2.22\right)$$

where ${\text{E}}_s$ denotes the expectation with respect to the sampling design, and $D$ a $KL\times KL$ diagonal matrix with diagonal elements

$$d_{kl}=\frac{1}{n_{kl}\left(n_{++}-1\right)}{\displaystyle \sum _{i=1}^{n_{++}}{\left(\frac{n_{++}\left(y_{ikl}-b_{kl}^t{x}_i\right)}{N{\pi }_i}-\frac{1}{n_{++}}{\displaystyle \sum _{i^{\prime }=1}^{n_{++}}\frac{n_{++}\left(y_{i^{\prime }kl}-b_{kl}^t{x}_{i^{\prime }}\right)}{N{\pi }_{i^{\prime }}}}\right)}^2},\left(2.23\right)$$

in the case of a CRD and

$$d_{kl}={\displaystyle \sum _{b=1}^B\frac{1}{n_{bkl}\left(n_{b++}-1\right)}{\displaystyle \sum _{i=1}^{n_{b++}}{\left(\frac{n_{b++}\left(y_{ikl}-b_{kl}^t{x}_i\right)}{N{\pi }_i}-\frac{1}{n_{b++}}{\displaystyle \sum _{i\text{'}=1}^{n_{b++}}\frac{n_{b++}\left(y_{i\text{'}kl}-b_{kl}^t{x}_{i\text{'}}\right)}{N{\pi }_{i\text{'}}}}\right)}^2}}\text{,}\left(2.24\right)$$

in the case of an RBD. An estimator for $D$ can be derived from the experimental design, conditionally on the measurement error model and the sampling design. Therefore the covariance matrix (2.22) is conveniently stated implicitly as the expectation over the measurement error model and the sampling design. A design-based estimator for this covariance matrix is given by

$$c\widehat{o}v(C{\widehat{\overline{Y}}}_{GREG})={\text{E}}_m{\text{E}}_{s}C\widehat{D}{C}^t\text{,}\left(2.25\right)$$

with $\widehat{D}$ a $KL\times KL$ diagonal matrix with elements

$${\widehat{d}}_{kl}=\frac{1}{n_{kl}(n_{kl}-1)}{\displaystyle \sum _{i=1}^{n_{kl}}{\left(\frac{n_{++}(y_{ikl}-{\widehat{b}}_{kl}^t{x}_i)}{N{\pi }_i}-\frac{1}{n_{kl}}{\displaystyle \sum _{i\text{'}=1}^{n_{kl}}\frac{n_{++}(y_{i\text{'}kl}-{\widehat{b}}_{kl}^t{x}_{i\text{'}})}{N{\pi }_{i\text{'}}}}\right)}^2},\left(2.26\right)$$

in the case of a CRD and

$${\widehat{d}}_{kl}={\displaystyle \sum _{b=1}^B\frac{1}{n_{bkl}(n_{bkl}-1)}{\displaystyle \sum _{i=1}^{n_{bkl}}{\left(\frac{n_{b++}(y_{ikl}-{\widehat{b}}_{kl}^t{x}_i)}{N{\pi }_i}-\frac{1}{n_{bkl}}{\displaystyle \sum _{i\text{'}=1}^{n_{bkl}}\frac{n_{b++}(y_{i\text{'}kl}-{\widehat{b}}_{kl}^t{x}_{i\text{'}})}{N{\pi }_{i\text{'}}}}\right)}^2}},\left(2.27\right)$$

in the case of an RBD. Proofs for (2.22) and (2.25) are given by van den Brakel (2010) and resemble the derivation of the covariance matrix for single factor experiments, given by van den Brakel and Renssen (2005) and van den Brakel(2008).

The results for (2.22) and(2.25) are obtained under the condition that a constant H-vector $a$ exists such that $a^t{x}_i=1$ for all $i\in U$. This is a rather weak condition, since it implies that a weighting model is used that at least uses the size of the finite population as a priori information. See van den Brakel and Renssen (2005) or van den Brakel (2008) for a more detailed discussion.

Since the $KL$ subsamples are drawn without replacement from a finite population, there is a nonzero design covariance between elements of ${\widehat{\overline{Y}}}_{GREG}$. From that point of view, it is remarkable that (2.25) has a structure as if the subsamples are drawn independently through sampling with replacement using unequal selection probabilities. This gives rise to an attractive variance estimation procedure for embedded experiments, since no design covariances between the subsample estimates appear in (2.25) and no second order inclusion probabilities are required in the variance estimators (2.26) and (2.27). This result is obtained since the covariance matrix of the contrasts between ${\widehat{\overline{Y}}}_{GREG}$ is derived instead of the covariance matrix of ${\widehat{\overline{Y}}}_{GREG}$  itself. A detailed interpretation of this result is given by van den Brakel and Renssen (2005) or van den Brakel (2008). See van den Brakel and Binder (2000) and Hidiroglou and Lavallée (2005) for approximations of the covariance matrix of ${\widehat{\overline{Y}}}_{GREG}$

The design-based estimators ${\widehat{\overline{Y}}}_{GREG}$ and $c\widehat{o}v(C{\widehat{\overline{Y}}}_{GREG})$ can be used to construct a design-based Wald statistic to test the hypotheses described in section 2.2:

$$W={\widehat{\overline{Y}}}_{GREG}^t{C}^t{\left(C\widehat{D}{C}^t\right)}^{-1}C{\widehat{\overline{Y}}}_{GREG}.\left(2.28\right)$$

Design-based inferences are generally based on normal large-sample approximations to construct confidence intervals for point estimates or p-values and critical regions for test statistics. Under this approach it follows under the null hypothesis that the Wald statistic is asymptotically distributed as a central chi-squared random variable, where the number of degrees of freedom equals the number of contrasts specified in the hypothesis.

The Wald statistic for the hypotheses about the main effects and interactions are given by (2.28) using the contrast matrix $C_A$, $C_B$, or $C_{AB}$. Under the null hypothesis, it follows that $W\to {\chi }_{[K-1]}^2$ for the test about the main effects of factor $A$, $W\to {\chi }_{[L-1]}^2$ for the test about the main effects of factor $B$ and $W\to {\chi }_{[(K-1)(L-1)]}^2$ for the test about interactions, where ${\chi }_{[p]}^2$ denotes a central chi-squared distributed random variable with p degrees of freedom.

The Wald test for the main effects can be further simplified. Expressions are developed for the Wald test for the main effects for factor $A$. Similar expressions can be derived for the main effects of factor $B$. Denote

$${\widehat{\overline{Y}}}_{A;GREG}={({\widehat{\overline{Y}}}_{1.;greg},\mathrm{...},{\widehat{\overline{Y}}}_{K.;greg})}^t,with{\widehat{\overline{Y}}}_{k.;greg}=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\widehat{\overline{Y}}}_{kl;greg}}$$

$${\widehat{D}}_A=\text{Diag}\left({\widehat{d}}_{1.}\text{,}\dots ,{\widehat{d}}_{K.}\right)\text{, with }{\widehat{d}}_{k.}=\frac{1}{L^2}{\displaystyle \sum _{l=1}^L{\widehat{d}}_{kl}}.\left(2.29\right)$$

It follows that $C_A{\widehat{\overline{Y}}}_{GREG}={\tilde{C}}_A{\widehat{\overline{Y}}}_{A;GREG}$ and $C_A\widehat{D}{C}_A^t={\tilde{C}}_A{\widehat{D}}_A{\tilde{C}}_A^t$. With the matrix inversion lemma, the Wald statistic for the main effects of factor $A$ can be simplified to:

$$\begin{array}{cc}\hfill W& ={\widehat{\overline{Y}}}_{A;GREG}^t{\tilde{C}}_A^t{\left({\tilde{C}}_A{\widehat{D}}_A{\tilde{C}}_A^t\right)}^{-1}{\tilde{C}}_A{\widehat{\overline{Y}}}_{A;GREG}^{}\hfill \\ & ={\widehat{\overline{Y}}}_{A;GREG}^t\left({\widehat{D}}_A^{-1}-\frac{1}{\text{Trace}({\widehat{D}}_A^{-1})}{\widehat{D}}_A^{-1}{j}_{(K-1)}{j}_{(K-1)}^t{\widehat{D}}_A^{-1}\right){\widehat{\overline{Y}}}_{A;GREG}^{}\left(2.30\right)\hfill \\ & ={\displaystyle \sum _{k=1}^K\frac{{\widehat{\overline{Y}}}_{k.;greg}^2}{{\widehat{d}}_{k.}}}-{\left({\displaystyle \sum _{k=1}^K\frac{1}{{\widehat{d}}_{k.}}}\right)}^{-1}{\left({\displaystyle \sum _{k=1}^K\frac{{\widehat{\overline{Y}}}_{k.;greg}^2}{{\widehat{d}}_{k.}}}\right)}^2.\hfill \end{array}$$

Finally note that the HT estimator (2.19) does not meet the condition that a constant H-vector $a$ exists such that $a^t{x}_i=1$ for all $i\in U$. The minimum use of auxiliary information used in the GREG estimator is obtained with a weighting scheme that only uses the size of the finite population as a priori knowledge, i.e. $(x_i)=1$ and ${\omega }_i^2={\omega }^2$  (Särndal et al. 1992, section 7.4). Under this weighting scheme it follows that

$${\widehat{\overline{Y}}}_{kl;greg}={\left({\displaystyle \sum _{i=1}^{n_{kl}}\frac{1}{{\pi }_i^{*}}}\right)}^{-1}\left({\displaystyle \sum _{i=1}^{n_{kl}}\frac{y_{ikl}}{{\pi }_i^{*}}}\right)\equiv {\tilde{y}}_{kl},\left(2.31\right)$$

and $({\widehat{b}}_{kl})={\tilde{y}}_{kl}$ . Expression (2.31) can be recognized as Hájek's ratio estimator for a population mean, (Hájek 1971). This weighting scheme satisfies the condition that a constant H-vector $a$ exists such that $a^t{x}_i=1$ for all $i\in U$. Therefore an approximately design-unbiased estimator for the covariance matrix of the contrasts between subsample estimates is given by (2.26) and (2.27) for a CRD and an RBD respectively, where ${\widehat{b}}_{kl}^t{x}_i={\tilde{y}}_{kl}$. Estimator (2.31) is preferable above the HT estimator (2.19), since (2.31) is more stable and the covariance matrix of the contrasts between (2.31) always has the relatively simple form of (2.25).

2.4 Special cases

It will be shown for two special cases that the design-based Wald statistic is equal to the F-test of a standard analysis of variance. Therefore, an ANOVA-type  pooled variance estimator for the diagonal elements of $\widehat{D}$ should be considered as an alternative for (2.26) or (2.27). Such a pooled variance estimator for a CRD is given by

$${\widehat{d}}_{kl}^p=\frac{1}{n_{kl}(n_{++}-KL)}{\displaystyle \sum _{k\text{'}=1}^K{\displaystyle \sum _{l\text{'}=1}^L{\displaystyle \sum _{i=1}^{n_{k\text{'}l\text{'}}}{\left(\frac{n_{++}(y_{ik\text{'}l\text{'}}-{\widehat{b}}_{k\text{'}l\text{'}}^t{x}_i)}{N{\pi }_i}-\frac{1}{n_{k\text{'}l\text{'}}}{\displaystyle \sum _{i\text{'}=1}^{n_{k\text{'}l\text{'}}}\frac{n_{++}(y_{i\text{'}k\text{'}l\text{'}}-{\widehat{b}}_{k\text{'}l\text{'}}^t{x}_{i\text{'}})}{N{\pi }_{i\text{'}}}}\right)}^2}}},\left(2.32\right)$$

and for an RBD by

$${\widehat{d}}_{kl}^p={\displaystyle \sum _{b=1}^B\frac{1}{n_{bkl}(n_{b++}-KL)}{\displaystyle \sum _{k\text{'}=1}^K{\displaystyle \sum _{l\text{'}=1}^L{\displaystyle \sum _{i=1}^{n_{bk\text{'}l\text{'}}}{\left(\frac{n_{b++}(y_{ik\text{'}l\text{'}}-{\widehat{b}}_{k\text{'}l\text{'}}^t{x}_i)}{N{\pi }_i}-\frac{1}{n_{bk\text{'}l\text{'}}}{\displaystyle \sum _{i\text{'}=1}^{n_{bk\text{'}l\text{'}}}\frac{n_{b++}(y_{i\text{'}k\text{'}l\text{'}}-{\widehat{b}}_{k\text{'}l\text{'}}^t{x}_{i\text{'}})}{N{\pi }_{i\text{'}}}}\right)}^2}}}}.\left(2.33\right)$$

Now consider a CRD that is embedded in a self-weighted sample, i.e. ${\pi }_i=n_{++}/N$, with equally sized subsamples, i.e. $n_{kl}=n_{k\text{'}l\text{'}}=n_s$. The inclusion probabilities for all units in the $KL$ subsamples are given by ${\pi }_i^{*}=n_s/N$. Let $\overline{y}=(1/n_s){\displaystyle {\sum }_{i=1}^{n_s}{y}_{ikl}}$. Under Hájek's ratio estimator (2.31) and the pooled variance estimator (2.32) it follows that ${\widehat{\overline{Y}}}_{kl;greg}={\overline{y}}_{kl}$, ${\widehat{b}}_{kl}={\overline{y}}_{kl}$, and

$${\widehat{d}}_{kl}^p=\frac{1}{n_s\left(n_{++}-KL\right)}{\displaystyle \sum _{k^{\prime }=1}^K{\displaystyle \sum _{l^{\prime }=1}^L{\displaystyle \sum _{i=1}^{n_s}{\left(y_{i{k}^{\prime }{l}^{\prime }}-{\overline{y}}_{k^{\prime }{l}^{\prime }}\right)}^2}}}\equiv \frac{{\widehat{S}}_{p;\text{ CRD}}^2}{n_s}.$$

The parameter estimates of the $K$ levels of factor $A$ averaged over the $L$ levels of factor $B$ are denoted as

$${\overline{y}}_{k.}=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\overline{y}}_{kl}}=\frac{1}{n_{k+}}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{i=1}^{n_s}{y}_{ikl}}},k=1,\dots ,K,\left(2.34\right)$$

with $n_{k+}={\displaystyle {\sum }_{l=1}^L{n}_{kl}}$. The diagonal elements of ${\widehat{D}}_A$ are now given by

$${\widehat{d}}_{k.}^p=\frac{1}{L^2}{\displaystyle \sum _{l=1}^L{\widehat{d}}_{kl}^p}=\frac{1}{L^2}{\displaystyle \sum _{l=1}^L\frac{{\widehat{S}}_{p;\text{ CDR}}^2}{n_s}}=\frac{{\widehat{S}}_{p;\text{ CDR}}^2}{n_{k+}},k=1,\dots ,K.\left(2.35\right)$$

Let ${\overline{y}}_{\mathrm{..}}=(1/n_{++}){\displaystyle {\sum }_{k=1}^K{\displaystyle {\sum }_{l=1}^L{\displaystyle {\sum }_{i=1}^{n_s}{y}_{ikl}}}}$. Inserting (2.34) and (2.35) into (2.30), gives rise to the following expression for the Wald statistic of the main effects of factor $A$

$$W=\frac{1}{{\widehat{S}}_{p;CRD}^2}\left({\displaystyle \sum _{k=1}^K{n}_{k+}{\overline{y}}_{k.}^2-n_{++}{\overline{y}}_{\mathrm{..}}^2}\right).\left(2.36\right)$$

Note that $W/(K-1)$  in (2.36) corresponds with the F-statistic for the main effects of an analysis of variance for the two-way layout with interactions, (Scheffé 1959, chapter 4). Under the null hypothesis and the assumption of normally and independently distributed errors, the F-statistic in the two-way layout follows an F-distribution with $(K-1)$ and $(n_{++}-KL)$ degrees of freedom, which is denoted as $F_{[n_{++}-KL]}^{[K-1]}$. If $n_{++}\to \infty$, then $F_{[n_{++}-KL]}^{[K-1]}\to {\chi }_{[K-1]}^2/(K-1)$. Consequently the F-statistic and the Wald statistic have the same limit distribution.

Now consider an RBD that is embedded in a self-weighted sampling design with equal subsample sizes, thus ${\pi }_i=n_{+++}/N$ and $n_{kl}=n_{k\text{'}l\text{'}}=n_s$, with $n_{+++}={\displaystyle {\sum }_{b=1}^B{n}_{b++}}$. Let ${\overline{y}}_{bkl}=(1/n_{bkl}){\displaystyle {\sum }_{i=1}^{n_{bkl}}{y}_{ikl}}$. Furthermore, it is assumed that the fraction of sampling units assigned to each treatment combination within each block is equal, i.e. $n_{bkl}/n_{b++}=n_s/n_{+++}$, and that the block sizes are sufficiently large to assume that $n_{b++}/(n_{b++}-KL)\approx 1$. Under Hájek's ratio estimator (2.31) and the pooled variance estimator (2.33) it follows that ${\widehat{\overline{Y}}}_{kl;greg}={\overline{y}}_{kl}$, ${\widehat{b}}_{kl}={\overline{y}}_{kl}$, and

$$\begin{array}{c}{\widehat{d}}_{kl}^p={\displaystyle \sum _{b=1}^B\frac{1}{n_{bkl}(n_{b++}-KL)}{\left(\frac{n_{b++}}{n_{+++}}\right)}^2{\displaystyle \sum _{k\text{'}=1}^K{\displaystyle \sum _{l\text{'}=1}^L{\displaystyle \sum _{i=1}^{n_{bk\text{'}l\text{'}}}{\left(y_{ik\text{'}l\text{'}}-{\overline{y}}_{bk\text{'}l\text{'}}\right)}^2}}}}\\ \approx \frac{1}{n_s{n}_{+++}}{\displaystyle \sum _{b=1}^B{\displaystyle \sum _{k\text{'}=1}^K{\displaystyle \sum _{l\text{'}=1}^L{\displaystyle \sum _{i=1}^{n_{bk\text{'}l\text{'}}}{\left(y_{ik\text{'}l\text{'}}-{\overline{y}}_{bk\text{'}l\text{'}}\right)}^2}}}}\equiv \frac{{\widehat{S}}_{p;RBD}^2}{n_s}.\end{array}$$

The parameter estimates of the $K$ levels of factor $A$ averaged over the $L$ levels of factor $B$ and the blocks are denoted as

$${\overline{y}}_{.k.}=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\overline{y}}_{kl}}=\frac{1}{n_{+k+}}{\displaystyle \sum _{b=1}^B{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{i=1}^{n_{bkl}}{y}_{ikl}}}},k=1,\dots ,K,\left(2.37\right)$$

where $n_{+k+}={\displaystyle {\sum }_{b=1}^B{\displaystyle {\sum }_{l=1}^L{n}_{bkl}}}$. The diagonal elements of ${\widehat{D}}_A$ are given by

$${\widehat{d}}_{k.}^p=\frac{1}{L^2}{\displaystyle \sum _{l=1}^L{\widehat{d}}_{kl}^p}=\frac{{\widehat{S}}_{p;\text{ RBD}}^2}{n_{+k+}},k=1,\dots ,K.\left(2.38\right)$$

Let ${\overline{y}}_{\mathrm{...}}=(1/n_{+++}){\displaystyle {\sum }_{b=1}^B{\displaystyle {\sum }_{k=1}^K{\displaystyle {\sum }_{l=1}^L{\displaystyle {\sum }_{i=1}^{n_{bkl}}{y}_{ikl}}}}}$. If these results are inserted into (2.30), then the expression for the Wald statistic of the main effects of factor $A$ can be simplified to

$$W=\frac{1}{{\widehat{S}}_{p;\text{ RBD}}^2}\left({\displaystyle \sum _{k=1}^K{n}_{+k+}{\overline{y}}_{.k.}^2-n_{+++}{\overline{y}}_{\dots }^2}\right).\left(2.39\right)$$

It can be recognized that $W/(K-1)$  in (2.39) corresponds with the F-statistic for the main effects of an analysis of variance for the three-way layout with interactions, (Scheffé 1959, chapter 4). As in the case of a CRD, this Wald and F-statistic have the same limit distribution.

Previous | Next