Adjusted Design Effect Model for Longitudinal Survey Data

Abstract

In this research, we propose a longitudinal adjusted design effect model that integrates both the design effect and time effect, offering reliable estimators of variance. We investigate the asymptotic properties of certain types of estimators. Through simulation studies, we show the exceptional performance of the proposed method, demonstrating its superiority in producing accurate standard error estimators compared to existing approaches.

1. Introduction

To calculate variance estimates in survey samples, commonly used methods include Taylor Series Linearization (TSL), jackknife, and Balanced Repeated Replication (BRR; interested readers can refer to Cohen (1979), Burt and Cohen (1984), Rao and Wu (1988), Rao (1988) and Wolter (2007), etc.). However, in many large-scale sample surveys such as the Current Population Survey (CPS) and the Canadian Labour Force Survey (CLFS), calculating thousands of variance estimates using these standard methods requires a significant amount of labor.

The Generalized Variance Functions (GVFs) have been used by CPS since 1947 (U.S. Census Bureau 2006). GVF develops regression estimates using a set of variance estimates obtained through BRR or TSL methods, along with their associated survey statistics. It allows for the direct calculation of standard errors for a large volume of potential survey statistics through a fitted regression model. Valliant (1987) demonstrated that the GVF model produces consistent estimates of the variance for a certain class of superpopulation models. However, while many current surveys follow the same households at regular time intervals, GVF does not take advantage of the longitudinal feature of the data.

Zhang et al. (2019) proposed a longitudinal generalized variance function (LGVF) by incorporating the time effect into the modeling of longitudinal survey data. The proposed LGVF maintains the property of GVF, where the estimated variances are often more stable than the direct estimates as they smooth out some of the variability from variable to variable when the design effects (determined by comparing the variance of a variable under a particular survey design with the variance under a simple random sampling [SRS] design with the same size) for the variables are similar to each other. However, a challenge faced by both GVF and LGVF is selecting a group of variables with similar design effects (deffs) to build a regression line, which is often practically difficult. Johnson and King (1986) studied GVF estimators using a national survey of reading ability among young adults and found that one way to markedly improve upon the GVF model is to use prior information about the design effect (deff) of individual estimators.

Building upon the groundwork laid by Zhang et al. (2019) and Johnson and King (1986), this research introduces the Longitudinal Adjusted Design Effect model (LADE). The LADE integrates both design and time effects using a model-based approach. The structure of the paper is as follows: Section 2 provides a comprehensive review of Generalized Variance Function (GVF) and Longitudinal Generalized Variance Function (LGVF) models; Section 3 establishes the theoretical framework, introduces the LADEs, and outlines the asymptotic properties of these models; Section 4 presents results from simulation studies conducted with CPS data. The research concludes in Section 5.

2. Generalized Variance Functions (GVFs) and Longitudinal Generalized Variance Functions (LGVFs)

In this section, we briefly review GVF and LGVF models. More detailed description of GVF can be found from textbooks by Wolter (2007), and Lohr (2021). Detailed description of LGVF can be found from Zhang et al. (2019).

2.1. Generalized Variance Functions (GVF)

Let ${\hat{T}}_{\upsilon }$ represent a survey statistic for variable $\upsilon$ , such as the estimated number of persons employed. Define ${\hat{p}}_{\upsilon }$ as the estimated proportion of employment, calculated as ${\hat{p}}_{\upsilon }={\hat{T}}_{\upsilon }/M$ , where $M$ is the population size. Consider $d_{\upsilon }$ as the design effect of ${\hat{p}}_{\upsilon }$ and $m$ as the sample size. We denote the sampling variance of ${\hat{p}}_{\upsilon }$ based on a given design as $\mathrm{var}({\hat{p}}_{\upsilon })=d_{\upsilon }\times p_{\upsilon }(1-p_{\upsilon })/m.$ Introduce the relative variance (relvar) of ${\hat{p}}_{\upsilon }$ as

$$\mathrm{relvar}({\hat{p}}_{\upsilon })=\frac{\mathrm{var}({\hat{p}}_{\upsilon })}{{[E({\hat{p}}_{\upsilon })]}^2}=a_{\upsilon }+\frac{b_{\upsilon }}{E({\hat{T}}_{\upsilon })},$$ (1)

where $a_{\upsilon }=‒d_{\upsilon }/m$ , and $b_{\upsilon }=M{d}_{\upsilon }/m$ . The relative variance is often recognized as the square of the coefficient of variation (CV). Let $V$ denote the number of variables under consideration. We propose a regression model linking a set of relative variances $\mathrm{relvar}({\hat{p}}_{\upsilon })$ to $E({\hat{T}}_{\upsilon })$ as follows:

$$\mathrm{relvar}({\hat{p}}_{\upsilon })=a+b/E({\hat{T}}_{\upsilon })\mathrm{for}\upsilon =1,2,\cdots ,V$$

assuming that $d_{\upsilon }$ values across these variables are the same, that is, $a_{\upsilon }=a$ and $b_{\upsilon }=b$ . Let $\hat{a}$ and $\hat{b}$ be the regression estimators of $a$ and $b$ . The GVF relative variance is predicted by $\hat{a}+\hat{b}/{\hat{T}}_{\upsilon }$ . A GVF estimate for $\mathrm{var}(\hat{T})$ is given by the following function,

$$\widehat{\mathrm{var}}({\hat{T}}_{\upsilon })=\hat{a}{\hat{T}}_{\upsilon }^2+\hat{b}{\hat{T}}_{\upsilon }.$$ (2)

2.2. LGVF

Zhang et al. (2019) introduced the Longitudinal Generalized Variance Function (LGVF), incorporating longitudinal information and addressing for time effects. Let $V$ denote the number of variables in the grouping, and $\tau$ represent the number of considered time periods, resulting in $V·\tau$ observations for LGVF. The subsequent discussion pertains to $t=1,2,\dots ,\tau$ , and $\upsilon =1,2,\dots ,V$ . To capture the time effect, we define

$$e_t=M_t/\overline{M},\mathrm{where}\overline{M}=(M_1+M_2+\cdots +M_{\tau })/\tau ,$$ (3)

where $M_t$ denotes the population size at time $t$ .

Let ${\hat{T}}_{t\upsilon }$ and ${\hat{p}}_{t\upsilon }$ denote the interested totals and proportions for variable $\upsilon$ at time $t$ . Let $\mathrm{var}({\hat{p}}_{t\upsilon })$ and ${\vartheta }_{t\upsilon }$ denote the variance and direct estimator of relative variance of ${\hat{p}}_{t\upsilon }$ based on the sampling design. Let $d_{t\upsilon }$ be the deff of variable $\upsilon$ at time $t$ , $a_{t\upsilon }=‒d_{t\upsilon }/m,b_{t\upsilon }=\overline{M}{d}_{t\upsilon }/m$ , and assume $a_{t\upsilon }=a_t=a$ and $b_{t\upsilon }=b_t=b$ . For estimating the relative variance, Zhang et al. (2019) proposed the LGVF model,

$$E({\vartheta }_{t\upsilon })=a+b·\frac{e_t}{T_{t\upsilon }}.$$

The predicted relative variance by LGVF value is

$${\hat{\vartheta }}_{t\upsilon }=\hat{a}+\hat{b}·\frac{e_t}{{\hat{T}}_{t\upsilon }}.$$ (4)

Under certain conditions, the ratios of relative variances and predicted relative variances from LGVF converge in probability to 1. The study by Zhang et al. (2019) demonstrated that LGVF was more efficient in reducing mean squared prediction errors (MSPE) compared to GVF.

3. The Longitudinal Adjusted Design Effects Models

In this section, we begin by introducing notation that closely adheres to conventions set by Valliant (1987) and Zhang et al. (2019). Subsequently, we present the Longitudinal Adjusted Design Effect models (LADEs), seamlessly integrating both the design effect and time effect. Finally, we explore the properties exhibited by specific types of estimators within the framework of LADE.

3.1. Notation

In a stratified two-stage cluster sampling design, various indices are employed to represent distinct levels within the sampling process. Specifically, $h$ serves as the stratum index, $i$ as the primary sampling unit (PSU) index, and $j$ as the secondary sampling unit (SSU) index, representing elements within the PSU. For instance, consider neighborhoods within a county as strata. In the first stage, households (PSUs) within the neighborhood are randomly selected, and in the second stage, some of the individuals (SSUs) residing in the chosen households are randomly interviewed. This approach facilitates a more manageable data collection process and ensures an accurate representation of diverse areas within the broader population.

At the PSU level, we denote $N_t$ as the total number of psus in the population at time $t$ , and $N_{\mathrm{th}}$ as the number of psus in stratum $h$ at time $t$ . Hence, the total number of PSUs in the population at time $t$ is given by $N_t=\sum _{h=1}^H{N}_{\mathrm{th}}$ . At the SSU level, we use $M_{\mathrm{thi}}$ to represent the number of SSUs in psu $i$ within stratum $h$ at time $t$ . Therefore, the total number of elements in stratum $h$ at time $t$ is denoted by $M_{\mathrm{th}}=\sum _{i=1}^{N_h}{M}_{\mathrm{thi}}$ , and the total number of SSUs in the population at time $t$ is denoted by $M_t=\sum _{h=1}^H{M}_{\mathrm{th}}$ .

Moving to the sample level, we denote $n_{\mathrm{th}}$ as the number of PSUs in the sample within stratum $h$ , and $n_t$ as the number of PSUs in the sample at time $t$ . Hence, the total number of PSUs in the sample at time $t$ is given by $n_t=\sum _{h=1}^H{n}_{\mathrm{th}}$ . We assume that the sample size remains constant over time, that is, $n_t=n$ for $t=1,2,\dots ,\tau$ . The assumption of a constant sample size over time is a convenient simplification for theoretical considerations. In practice, researchers should be aware of the real-world complexities introduced by non-response and attrition. Furthermore, we use $m_{\mathrm{thi}}$ to represent the number of elements in the sample from the ith PSU within stratum $h$ . Therefore, the total number of elements in the sample within stratum $h$ at time $t$ is denoted by $m_{\mathrm{th}}=\sum _{i=1}^{n_{\mathrm{th}}}{m}_{\mathrm{thi}}$ . Consequently, the total number of elements in the sample over all strata at time $t$ is denoted by $m_t=\sum _{h=1}^H{m}_{\mathrm{th}}$ .

At time $t$ , we define ${\mathcal{S}}_{\mathrm{th}}$ as the set of sampled psus in stratum $h$ , ${\mathcal{R}}_{\mathrm{th}}$ as the set of nonsampled psus in stratum $h$ , ${\mathcal{S}}_{\mathrm{thi}}$ as the set of sampled elements within PSU $i$ in stratum $h$ , and ${\mathcal{R}}_{\mathrm{thi}}$ as the set of nonsampled elements within PSU $i$ in stratum $h$ . Assuming that a random variable $y_{t\upsilon \mathrm{hij}}$ is associated with the jth unit within PSU $i$ in stratum $h$ at time $t$ for variable $\upsilon$ , we can define a general type estimator of $T_{t\upsilon }$ as follows:

$${\hat{T}}_{t\upsilon }=\sum _h\sum _{i\in {\mathcal{S}}_{\mathrm{th}}}{\gamma }_{\mathrm{thi}}{\hat{T}}_{t\upsilon \mathrm{hi}}.$$ (5)

Here, ${\gamma }_{\mathrm{thi}}$ is a coefficient that depends on the specific survey design. Define ${\overline{y}}_{t\upsilon \mathrm{hi}}=\sum _{j\in {\mathcal{S}}_{\mathrm{thi}}}{y}_{t\upsilon \mathrm{hij}}/m_{\mathrm{thi}}$ , so that ${\hat{T}}_{t\upsilon \mathrm{hi}}=M_{\mathrm{thi}}{\overline{y}}_{t\upsilon \mathrm{hi}}$ .

For example, the Horvitz-Thompson estimator, which is commonly used when PSUs are selected with probabilities proportional to $M_{\mathrm{thi}}$ and an equal probability sample is selected within each sampled PSU at time $t$ , can be written as:

$${\hat{T}}_{\mathrm{tv},\mathrm{HT}}=\sum _h\sum _{i\in {\mathcal{S}}_{\mathrm{th}}}\left[M_{\mathrm{th}}{\hat{T}}_{t\upsilon \mathrm{hi}}{(n_{\mathrm{th}}{M}_{\mathrm{thi}})}^{-1}\right],$$ (6)

where ${\gamma }_{\mathrm{thi}}=M_{\mathrm{th}}(n_{\mathrm{th}}{M}_{\mathrm{thi}}{)}^{-1}$ .

For prediction purposes, the following model assumptions apply:

$$E(y_{tvhij})={\mu }_{tvh},\mathrm{cov}(y_{tvhij},y_{tv{h}^{\prime }{i}^{\prime }{j}^{\prime }})=(\begin{array}{cc}{\sigma }_{tvhi}^2& \mathrm{if}h=h^{\prime },i=i^{\prime },j=j^{\prime }\\ {\rho }_{tvhi}{\sigma }_{tvhi}^2& \mathrm{if}h=h^{\prime },i=i^{\prime },j\ne j^{\prime }\\ 0& \mathrm{otherwise}.\end{array}$$ (7)

where within a PSU $i$ in stratum $h$ at time $t$ , ${\sigma }_{t\upsilon \mathrm{hi}}$ signifies the variance of measurements and ${\rho }_{t\upsilon \mathrm{hi}}$ represents the intra-class correlation coefficient (ICC) for variable $\upsilon$ . The ICC measures the degree of correlation among observations within a cluster, indicating the level of homogeneity within clusters. A higher intra-class correlation implies greater similarity among observations within a cluster.

Similar formulations can be found in the works of Scott and Smith (1969), Royall (1976, 1986), and Burdick and Sielken (1979). A general direct variance estimator for $\mathrm{var}({\hat{T}}_{t\upsilon })$ used in this research is provided by Royall (1986):

$$s_{{\hat{T}}_{t\upsilon }}^2=\sum _h{n}_{\mathrm{th}}(n_{\mathrm{th}}-1{)}^{-1}\sum _{{\mathcal{S}}_{\mathrm{th}}}{\gamma }_{\mathrm{thi}}^2{r}_{t\upsilon \mathrm{hi}}^2,$$ (8)

where, $r_{t\upsilon \mathrm{hi}}={\hat{T}}_{t\upsilon \mathrm{hi}}-(\sum _{j\in {\mathcal{S}}_{\mathrm{th}}}{\gamma }_{\mathrm{thj}}{\hat{T}}_{t\upsilon \mathrm{hj}}/M_{\mathrm{ht}})M_{\mathrm{thi}}$ , and ${\gamma }_{\mathrm{thi}}$ is defined in Equation (6). The direct estimator of the relative variance is ${\vartheta }_{t\upsilon }=s_{{\hat{T}}_{t\upsilon }}^2/{\hat{T}}_{t\upsilon }^2$ .

Let $k_{t\upsilon \mathrm{hi}}=[1+(m_{\mathrm{thi}}-1){\rho }_{t\upsilon \mathrm{hi}}]/m_{\mathrm{thi}}$ , ${\alpha }_{1t\upsilon h}$ and ${\alpha }_{2t\upsilon h}$ are constants. Assume that

$${\sigma }_{t\upsilon \mathrm{hi}}^2={\sigma }_{t\upsilon h}^2={\alpha }_{1t\upsilon h}{\mu }_{\mathrm{th}}+{\alpha }_{2t\upsilon h}{\mu }_{t\upsilon h}^2.$$ (9)

For example, for a binary distribution, ${\alpha }_{1h}=1$ and ${\alpha }_{2h}=‒1$ since ${\sigma }_h^2={\mu }_h(1-{\mu }_h)={\mu }_h-{\mu }_h^2$ . According to the definition, $\mathrm{relvar}({\hat{T}}_{t\upsilon }-T_{t\upsilon })=\mathrm{var}({\hat{T}}_{t\upsilon }-T_{t\upsilon })/[E({\hat{T}}_{t\upsilon }){]}^2$ , where $\mathrm{var}({\hat{T}}_{t\upsilon })=\mathrm{var}\left(\sum _h\sum _{i\in {\mathcal{S}}_{\mathrm{th}}}{\gamma }_{\mathrm{thi}}{\hat{T}}_{t\upsilon \mathrm{hi}}\right)=\mathrm{var}(\sum _h\sum _{i\in {\mathcal{S}}_{\mathrm{th}}}{\gamma }_{\mathrm{thi}}{\hat{M}}_{\mathrm{thi}}{\overline{y}}_{t\upsilon \mathrm{hi}})$ . By expanding the variance and utilizing assumptions in Equation (7) and Equation (9), we can derive:

$$\begin{array}{c}\mathrm{relvar}({\hat{T}}_{t\upsilon }-T_{t\upsilon })\approx \sum _h{\pi }_{t\upsilon h}^2{\alpha }_{2t\upsilon h}{M}_{\mathrm{th}}^{-2}\sum _{{\mathcal{S}}_{\mathrm{th}}}{\gamma }_{\mathrm{thi}}^2{k}_{t\upsilon \mathrm{hi}}{M}_{\mathrm{thi}}^2\hfill \\ \multicolumn{1}{c}{+\left[\sum _h{\pi }_{t\upsilon h}{\alpha }_{1t\upsilon h}{M}_{\mathrm{th}}^{-2}\sum _{{\mathcal{S}}_{\mathrm{th}}}{\gamma }_{\mathrm{thi}}^2{k}_{t\upsilon \mathrm{hi}}{M}_{\mathrm{thi}}^2\right]/E({\hat{T}}_{t\upsilon })}\\ \multicolumn{1}{c}{=a_t+b_t/E({\hat{T}}_{t\upsilon }),}\end{array}$$ (10)

where ${\pi }_{t\upsilon h}=E({\hat{T}}_{t\upsilon h})/E({\hat{T}}_{t\upsilon })$ , and ${\hat{T}}_{t\upsilon h}=\sum _{i=1}^{N_h}{\hat{T}}_{t\upsilon \mathrm{hi}}$ .

3.2. The LADE Models

Recall that $\mathrm{relvar}({\hat{p}}_i)=‒d_i/m+M{d}_i/[\mathrm{mE}({\hat{T}}_i)]$ . If we divide both sides by $d_i$ , the ratio $\mathrm{relvar}({\hat{p}}_i)/d_i$ exhibits a perfect linear relationship with $1/E({\hat{T}}_i)$ . Figure 1 depicts the plots of the direct estimators of the relative variance ${\vartheta }_{2010\upsilon }$ (plot (a)) and the adjusted relative variance ${\vartheta }_{2010\upsilon }/{\hat{d}}_{2010\upsilon }$ (plot (b)) for eighteen variables against the reciprocal of the estimated total from a single simulation run using March 2010 CPS annual social and economic supplement (ASEC) data. For detailed information regarding the simulation setup, please consult Section 4. Notably, plot (b) clearly illustrates a perfect linear relationship between ${\vartheta }_{2010\upsilon }/{\hat{d}}_{2010\upsilon }$ and $1/{\hat{T}}_{2010\upsilon }$ , as expected.

Figure 1.

Direct estimates of relative variance (plot (a)) and adjusted direct estimates of relative variance (plot (b)) versus reciprocal of estimated total.

Motivated by this observation, we propose the Longitudinal Adjusted Design Effect (LADE) model, which adjusts the relative variance (relvar) by the design effect (deff) and incorporates time effect. To get an overview of the deffs of different variables, Figure 2 displays the design effects of the eighteen variables obtained from the mean of the five hundred simulation runs using the 2008 to 2010 ASEC data. It is observed that the design effects for the variables remain relatively stable over the three-year period. However, it is worth noting that the magnitudes of the design effects vary across the different variables.

Figure 2.

Design effects of the eighteen variables for 08 to 10 ASEC data from five hundred simulation runs. X-axis is the index of variables and Y-axis is the design effect for the variables from the five hundred simulation runs. For detailed information regarding the simulation setup, please consult Section 4.

To incorporate the design effects into the analysis, we adjust the relative variance by taking the mean of the design effects for each year across the eighteen variables. Let $\mathbf{\theta }=(a,b)\prime$ be the regression parameters. Let $a_{t\upsilon }=‒d_{t\upsilon }/m_t$ , $b_{t\upsilon }=\overline{M}{d}_{t\upsilon }/m_t$ , and $e_t$ be defined as in Equation (3). Using Equation (2), we can express the estimated variance of ${\hat{T}}_{t\upsilon }$ as:

$$\begin{array}{c}\multicolumn{1}{c}{\widehat{\mathrm{var}}({\hat{T}}_{t\upsilon })=\frac{-{\hat{d}}_{t\upsilon }}{m_t}{\hat{T}}_{t\upsilon }^2+\frac{M_t{\hat{d}}_{t\upsilon }}{m_t}{\hat{T}}_{t\upsilon }}\\ \multicolumn{1}{c}{=\frac{-{\hat{d}}_{t\upsilon }}{{\hat{\overline{d}}}_t}\frac{{\hat{\overline{d}}}_t}{m_t}{\hat{T}}_{t\upsilon }^2+\frac{M_t}{\overline{M}}\frac{{\hat{d}}_{t\upsilon }}{{\hat{\overline{d}}}_t}\frac{\overline{M}{\hat{\overline{d}}}_t}{m_t}{\hat{T}}_{t\upsilon }}\\ \multicolumn{1}{c}{={\hat{f}}_{t\upsilon }{\hat{a}}_{t\upsilon }^{*}{\hat{T}}_{t\upsilon }^2+e_t{\hat{f}}_{t\upsilon }{\hat{b}}_{t\upsilon }^{*}{\hat{T}}_{t\upsilon }}\end{array}$$

where $f_{t\upsilon }=d_{t\upsilon }/{\overline{d}}_t$ , ${\overline{d}}_t=\sum _{\upsilon =1}^V{d}_{t\upsilon }/V$ , $a_{t\upsilon }^{*}=‒{\overline{d}}_t/m_t$ , and $b_{t\upsilon }^{*}=\overline{M}{\overline{d}}_t/m_t$ ; ${\hat{f}}_{t\upsilon }$ , $\hat{\overline{d}}t,{\hat{a}}_{t\upsilon }^{*}$ , and ${\hat{b}}_{t\upsilon }^{*}$ are estimators of $f_{t\upsilon }$ , ${\overline{d}}_t$ , $a_{t\upsilon }^{*}$ , and $b_{t\upsilon }^{*}$ respectively. Therefore, $\widehat{\mathrm{var}}({\hat{T}}_{t\upsilon })/({\hat{T}}_{t\upsilon }^2*{\hat{f}}_{t\upsilon })={\hat{a}}_{t\upsilon }^{*}+{\hat{b}}_{t\upsilon }^{*}{e}_t/{\hat{T}}_{t\upsilon }$ . Define a set of adjusted estimators of relvars ${\vartheta }_{t\upsilon }^{*}$ as ${\vartheta }_{t\upsilon }^{*}={\vartheta }_{t\upsilon }/{\hat{f}}_{t\upsilon }$ for $t=1,2,\cdots ,\tau$ and $\upsilon =1,2,\cdots ,V,$ such that

$$E({\vartheta }_{t\upsilon }^{*})=a_{t\upsilon }^{*}+b_{t\upsilon }^{*}·\frac{e_t}{T_{t\upsilon }}.$$ (11)

Recall that ${\vartheta }_{t\upsilon }$ denote the variance and direct estimator of relative variance of ${\hat{p}}_{t\upsilon }$ based on the sampling design, and $E({\vartheta }_{t\upsilon })=a_0+b_0{e}_t/T_{t\upsilon }$

Let ${\mathbf{\vartheta }}_t^{*}=({\vartheta }_{t1}^{*},{\vartheta }_{t2}^{*},\cdots ,{\vartheta }_{\mathrm{tV}}^{*})\prime$ , ${\mathbf{\vartheta }}^{*}=({\mathbf{\vartheta }}_1^{*\prime },{\mathbf{\vartheta }}_2^{*\prime },\cdots ,{\mathbf{\vartheta }}_{\tau }^{*\prime })\prime$ , ${ϵ}_t=({ϵ}_{t1},{ϵ}_{t2},\cdots ,{ϵ}_{\mathrm{tV}})\prime$ , and $ϵ=({ϵ}_1\prime ,{ϵ}_2\prime ,\cdots ,{ϵ}_{\tau }\prime )\prime$ . Define ${\mathbf{X}}_t$ as the $V\times 2$ design matrix for time $t$ with the first column consisting of one second and second column being $(e_t/{\hat{T}}_{t1},e_t/{\hat{T}}_{t2},\cdots ,e_t/{\hat{T}}_{\mathrm{tV}})$ . Let $\mathbf{X}$ be the design matrix defined as $\mathbf{X}=(\mathbf{X}{\prime }_1,\cdots ,\mathbf{X}{\prime }_{\tau })\prime$ and let $\mathbf{\theta }=(a,b)\prime$ . Under the assumption that $a_{t\upsilon }^{*}=a_t^{*}=a$ and $b_{t\upsilon }^{*}=b_t^{*}=b$ , the LADE model in Equation (11) can be written in the matrix form as:

$${\mathbf{\vartheta }}^{*}=\mathbf{X}\mathbf{\theta }+ϵ.$$ (12)

The weighted least square estimators (WLS) of $\mathbf{\theta }$ are given by $\hat{\theta }=(\mathbf{X}\prime \mathbf{WX}{)}^{-1}\mathbf{X}\prime \mathbf{W}{\vartheta }^{*}$ , where $w_{t\upsilon }$ is a weight associated with variable $\upsilon$ at time $t$ , and $\mathbf{W}$ is a $V\tau \times V\tau$ matrix with the diagonal elements as $w_{t\upsilon }$ . Typically, $w_{t\upsilon }$ is chosen as the reciprocal of the variance of ${\vartheta }_{t\upsilon }^{*}$ when known. Alternatively, if the variances are unknown, we can approximate the weights by the reciprocal of the squared ${\vartheta }_{t\upsilon }^{*}$ . Following the work of Zhang et al. (2019), the estimators of $a$ and $b$ are given by,

$$\hat{b}=\frac{\sum _{t=1}^{\tau }\sum _{\upsilon =1}^V{\vartheta }_{t\upsilon }^{*}[e_t{\hat{T}}_{t\upsilon }^{-1}-{\overline{T}}_{-}]/w_{t\upsilon }}{\sum _{t=1}^{\tau }\sum _{\upsilon =1}^V{[e_t{\hat{T}}_{t\upsilon }^{-1}-{\overline{T}}_{-}]}^2/w_{t\upsilon }}={\hat{S}}_1/{\hat{S}}_2$$ (13)

and

$$\hat{a}={\overline{\vartheta }}^{*}-\hat{b}{\overline{T}}_{-}$$ (14)

where ${\overline{T}}_{-}=\sum _{t,\upsilon }(e_t^{-1}{\hat{T}}_{t\upsilon }{w}_{t\upsilon }{)}^{-1}/\sum _{t,\upsilon }{w}_{t\upsilon }^{-1}$ , ${\overline{\vartheta }}^{*}=\sum _{t,\upsilon }{\vartheta }_{t\upsilon }^{*}{w}_{t\upsilon }^{-1}/\sum _{t,\upsilon }{w}_{t\upsilon }^{-1}$ , ${\hat{S}}_1=\sum _{t=1}^{\tau }\sum _{\upsilon =1}^V{\vartheta }_{t\upsilon }^{*}[e_t{\hat{T}}_{t\upsilon }^{-1}-{\overline{T}}_{-}]/w_{t\upsilon }$ , and ${\hat{S}}_2=\sum _{t=1}^{\tau }\sum _{\upsilon =1}^V[e_t{\hat{T}}_{t\upsilon }^{-1}-{\overline{T}}_{-}{]}^2/w_{t\upsilon }$ . The predicted adjusted relvar of ${\hat{T}}_{t\upsilon }$ based on the estimated LADE model is

$${\hat{\vartheta }}_{t\upsilon }^{*}={\overline{\vartheta }}^{*}+\hat{b}[e_t{\hat{T}}_{t\upsilon }^{-1}-{\overline{T}}_{-}],\mathrm{and}{\hat{\vartheta }}_{t\upsilon }={\hat{\vartheta }}_{t\upsilon }^{*}{\hat{f}}_{t\upsilon }.$$ (15)

3.3. Properties of Proposed Estimators

In this section, we examine a specific class of estimators defined as ${\gamma }_{\mathrm{thi}}=g_{1\mathrm{th}}{g}_{2\mathrm{thi}}$ , representing the product of a nonrandom factor dependent solely on the stratum $\left(g_{1\mathrm{th}}\right)$ and another nonrandom factor dependent on both the stratum and cluster $\left(g_{2\mathrm{thi}}\right)$ . For instance, when PSUs are selected with probabilities proportional to stratum size $M_{\mathrm{thi}}$ , and an equal probability sample is chosen within each sampled PSU at time $t$ , we have $g_{1\mathrm{th}}=M_{\mathrm{th}}/n_{\mathrm{th}}$ and $g_{2\mathrm{thi}}=1/M_{\mathrm{thi}}$ .

Given the assumptions in Equation (7) and the structure of estimators ${\gamma }_{\mathrm{thi}}=g_{1\mathrm{th}}{g}_{2\mathrm{thi}}$ , we investigate the asymptotic properties of ${\hat{T}}_{t\upsilon }$ , $s_{{\hat{T}}_{t\upsilon }}^2$ , and ${\hat{\vartheta }}_{t\upsilon }$ when the number of PSUs in each stratum is large. Lemmas 1 to 3 (see Appendix) extend the work by Royall (1986). Under certain conditions, Theorem 1 demonstrates that the ratios of relative variances and predicted relative variances from the proposed LADEs converge in probability to 1. The proof of Theorem 1 is provided in the Appendix. The asymptotic normality then follows.

• Theorem 1. Under the model defined by Equation (7) with variance specification in Equation (9), the assumptions (i) to (xiv) listed in the Appendix, and ${\mu }_{4t\upsilon \mathrm{hi}}=E[{\hat{T}}_{t\upsilon \mathrm{hi}}-E({\hat{T}}_{t\upsilon \mathrm{hi}}){]}^4<\infty ,$ $a_{t\upsilon }^{*}=a_t^{*}=a$ and $b_{t\upsilon }^{*}=b_t^{*}=b$ for $t=1,2,\dots ,\tau$ , $\upsilon =1,2,\dots ,V$ , as $N_{\mathrm{th}},n_{\mathrm{th}}\to \infty$ ,

$$\frac{\mathrm{relvar}({\hat{T}}_{t\upsilon }-T_{t\upsilon })}{{\hat{\vartheta }}_{t\upsilon }}{\to }^{p}1.$$

• Proof. The proof is given in the Appendix. □

• Theorem 2. Under model in Equation (7) with variance specification in Equation (9), the assumptions (i) to (xiv) listed in the Appendix, and ${\mu }_{4\upsilon \mathrm{thi}}=E[{\hat{T}}_{t\upsilon \mathrm{hi}}-E({\hat{T}}_{t\upsilon \mathrm{hi}}){]}^4<\infty ,$ $a_{t\upsilon }^{*}=a_t^{*}=a$ and $b_{t\upsilon }^{*}=b_t^{*}=b$ for $t=1,2,\dots ,\tau$ , $\upsilon =1,2,\dots ,V$ , as $N_{\mathrm{th}},n_{\mathrm{th}}\to \infty$ ,

$$\frac{{\hat{T}}_{t\upsilon }-T_{t\upsilon }}{{\hat{T}}_{t\upsilon }{({\hat{\vartheta }}_{t\upsilon })}^{1/2}}{\to }^{d}N(0,1).$$

• Proof. The proof is a straightforward extension of work by Royall 1986. □

Theorem 2 enables the construction of confidence intervals using the normal distribution.

4. Simulation Studies

This section presents a simulation study that compares the performance of two types of estimators: LGVFs and LADEs. The study utilizes data from the CPS Annual Social and Economic Supplement (ASEC), specifically focusing on the data restricted to New Mexico for the years 2008 to 2011.

The analysis considers eighteen binary variables derived from the “Source of Income” section of the ASEC data. These variables capture characteristics such as self-employment, unemployment compensation, and more. In this context, a binary variable takes a value of 1 if a person possesses a specific characteristic and 0 if they do not possess that characteristic. The available ASEC data includes 2,059 observations for 2008, 2,188 observations for 2009, 2,108 observations for 2010, and 1,975 observations for 2011. Each household in the data is associated with an ultimate sampling unit (USU) defined by the CPS. However, the specific USU information is not publicly released.

To replicate the sampling design, the households within each year were sorted in ascending order based on their sequence numbers. Then, four consecutive households were combined to form a PSU while maintaining the original order. As a result, there were 205, 208, 193, and 184 PSUs created for the years 2008, 2009, 2010, and 2011, respectively. On average, each of these PSUs accommodates about ten individuals across the four households. Consequently, for example, the total number of individuals in 2008 can be calculated as the number of PSUs (205) multiplied by the average number of individuals residing in the four households, resulting in approximately 2,059 observations. The simulation study is conducted using data from March 2008 to March 2010 for bias square and empirical mean squared error (EMSE) comparisons, and data from March 2011 for mean squared prediction error (MSPE) comparisons. The simulation steps are as follows:

• (a) Within each year (2008, 2009, and 2010), a sample of $n=100$ PSUs is selected using probabilities proportional to size (PPS). Within each selected PSU, $m_i=4$ individuals are randomly selected without replacement.

• (b) Estimates ${\hat{T}}_{t\upsilon }$ for each sample year $t$ and variable $\upsilon$ are calculated using the Horvitz-Thompson estimator (Equation (6)), and the estimated variance $s_{{\hat{T}}_{t\upsilon }}^2$ is calculated using Equation (8). The direct estimator of relative variance is then calculated as ${\vartheta }_{t\upsilon }=s_{{\hat{T}}_{t\upsilon }}^2/({\hat{T}}_{t\upsilon }{)}^2$ .

• (c) The time adjustment, denoted as $e_t$ , is applied to the estimates ${\hat{T}}_{t\upsilon }$ by dividing each estimate by itself, that is, $e_t/{\hat{T}}_{t\upsilon }$ .

•   Here, $e_1=M_1/\overline{M}=1,978,390/1,967,487=1.0056(\mathrm{for}\mathrm{year}2009)$ ;

•    $e_2=M_2/\overline{M}=1,977,807/1,967,487=1.0052(\mathrm{for}\mathrm{year}2010)$ ; and $e_3=M_3/\overline{M}=1,946,264/1,967,487=0.9892(\mathrm{for}\mathrm{year}2008)$ .

• (d) Deff adjustment ${\hat{f}}_{t\upsilon }={\hat{d}}_{t\upsilon }/\widehat{\overline{d_t}}$ is applied to the relvar estimates ${\vartheta }_{t\upsilon }$ , so that ${\vartheta }_{t\upsilon }^{*}={\vartheta }_{t\upsilon }/{\hat{f}}_{t\upsilon }.$ Here, ${\hat{d}}_{t\upsilon }$ is the design effect of each variable and sample, and $\widehat{\overline{d_t}}$ is the average design effect calculated across all variables and simulation runs at time $t$ .

• (e) Two regression models are fitted. Model (16) is used to fit LGVFs, and Model (17) is used to fit LADEs.

$${\vartheta }_{t\upsilon }=a_0+b_0·\frac{e_t}{{\hat{T}}_{t\upsilon }},$$ (16)

and

$${\vartheta }_{t\upsilon }^{*}=a+b·\frac{e_t}{{\hat{T}}_{t\upsilon }}.$$ (17)

Fitting methods include ordinary linear regression (OLS) LGVF1 and LADE1 and WOLS with weights $w_{t\upsilon }=1/({\vartheta }_{t\upsilon }{)}^2$ for LGVF2 and $w_{t\upsilon }=1/({\vartheta }_{t\upsilon }^{*}{)}^2$ for LADE2.

• (f) The direct estimates of relvar calculated by Equations (6) and (8) is recorded, as well as the fitted relvar values by LGVF1, LGVF2, LADE1, and LADE2 multiplied by the deff adjustment $f_{t\upsilon }$ .

• (g) The simulation process is repeated $R=500$ times, following steps (a) to (f). During each iteration, the recorded relvars and fitted relvar values are collected. After the five hundred iterations, the average values of these recorded relvars and fitted relvar values are calculated. Subsequently, bias square, EMSE, and MSPE are computed using the provided formulas.

$${\mathrm{Bias}}^2=\sum _{t=1}^3\sum _{\upsilon =1}^{18}({\mathrm{Bias}}_{t\upsilon }{)}^2/54\mathrm{with}{\mathrm{Bias}}_{t\upsilon }=\frac{\sum _{r=1}^{500}{\hat{\vartheta }}_{t\upsilon }^{(r)}}{R}-{\vartheta }_{t\upsilon },$$ (18)

$$\mathrm{EMSE}=\sum _{t=1}^3\sum _{\upsilon =1}^{18}\mathrm{EMS}{E}_{t\upsilon }/54,\mathrm{with}\mathrm{EMS}{E}_{t\upsilon }=\frac{\sum _{r=1}^{500}{({\hat{\vartheta }}_{t\upsilon }^{(r)}-{\vartheta }_{t\upsilon })}^2}{R},$$ (19)

and

$$\mathrm{MSPE}=\sum _{r=1}^{500}{\mathrm{MSPE}}^{(r)}/R,\mathrm{with}\mathrm{MSP}{E}^{(r)}=\frac{\sum _{\upsilon =1}^{18}{({\hat{\vartheta }}_{(2011)\upsilon }^{(r)}-{\vartheta }_{(2011)\upsilon })}^2}{18},$$ (20)

where $(r)$ denotes the estimate from the rth simulation run.

The simulation results are summarized in Table 1 which compares the bias square, EMSE, and MSPE values for the estimators LGVF1, LGVF2, LADE1, and LADE2. The results indicate that LADEs consistently outperform LGVFs across all metrics, showcasing lower bias square, EMSE, and MSPE values. Notably, there is no discernible advantage in using the WLS estimators LGVF2 and LADE2 over the OLS estimators LGVF1 and LADE1. This absence of advantage can be attributed to the overall equal variance assumption of the regression model not being violated.

Table 1.

Simulation Results of a Comparison Between LGVFs and LADEs with $n_t=100$ . LGVF1 and LADE1 Are with Ordinary Linear Regression Fit, and LGVF2 and LADE2 Are with Weighted Least Square Fit.

	LGVF1	LGVF2	LADE1	LADE2
Bias	0.003250381	0.003314496	0.001336026	0.001165163
EMSE	0.012627135	0.017065547	0.007332843	0.008725467
MSPE	0.03385114	0.03887468	0.01813410	0.01679681

Figure 3 provides a graphical representation of the estimated relvars plotted against the logs of population totals. The solid line represents the direct estimates of relvar values calculated by Equations (6) and (8). The dashed line and dotted line represent the estimates from LGVFs 1 to 2, while the dot-dashed line and the thickest line represent the estimates from LADEs 1 to 2. It is observed that the WLS estimators outperform the OLS estimators in estimating the relvars, especially when the population totals of the variables are large. Figure 3 suggests that heteroscedasticity is more pronounced in the range of large values. WLS is then able to better handle the heteroscedasticity, leading to more efficient and less biased estimates of the relative variance compared to OLS.

Figure 3.

Logs of estimates of relvar plotted versus logs of population totals.

Although LADE2 exhibits a slight improvement over LGVF2, the distinction between LGVFs and LADEs is not evident as that in Table 1. This discrepancy can be attributed to the differing perspectives and scaling from which Table 1 and Figure 3 approach the data. Table 1 offers a comprehensive evaluation of bias in relative variance, whereas Figure 3 explores the intricate connection between the logarithm of relative variance and log totals.

5. Conclusions

In this study, we extended the Longitudinal Generalized Variance Functions (LGVFs) to a novel model known as Longitudinal Adjusted Design Effect models (LADEs), incorporating both design effects and time effects into the modeling process. We demonstrated that, under certain conditions, the ratios of relative variances and predicted relative variances obtained from LADEs converge to 1. A comparative analysis between LGVFs and LADEs revealed that the latter outperforms in reducing bias, empirical mean squared error, and mean squared prediction errors. While our simulation study supports the superiority of LADE over LGVF, we acknowledge its limitations and encourage users to carefully select the model based on their specific context.

Our research hasn’t specifically addressed survey nonresponse, but techniques like multiple imputation and weighted adjustments for regression methods can be seamlessly integrated into our framework. Notably, our analysis is currently confined to estimators derived from the LADE model without auxiliary variables in the regression model.

For future investigations, exploring mixed models and nonparametric smoothing methods for fitting regression models could be valuable. Additionally, a detailed examination of estimators in cases where variances within one PSU differ from another PSU could offer insightful perspectives.