3 Methodology

Iván A. Carrillo and Alan F. Karr

Previous | Next

3.1 Motivation

Assume that (in a non-survey context) interest lies in the $p\times 1$ vector parameter $\beta$ in the following model:

$\xi :(\begin{array}{ll}E\left[Y_{ij}|X_{ij}\right]={\mu }_{ij}=g^{-1}({X^{\prime }}_{ij}\beta )\mathrm{,}\hfill & j=\mathrm{1,2,}\dots \mathrm{,}J\mathrm{,}i=\mathrm{1,2,}\dots \hfill \\ \text{Var}\left[Y_{ij}|X_{ij}\right]=\varphi \nu ({\mu }_{ij})\mathrm{,}\hfill & j=\mathrm{1,2,}\dots \mathrm{,}J\mathrm{,}i=\mathrm{1,2,}\dots \hfill \\ \text{Cov}\left[Y_i|X_i\right]={\Sigma }_i\mathrm{,}\hfill & i=\mathrm{1,2,}\dots \hfill \\ Y_k\perp Y_l|X_k\mathrm{,}{X}_l\mathrm{,}\hfill & k\ne l=\mathrm{1,2,}\dots ;\hfill \end{array}$(3.1)

where $Y_{ij}$ is the response variable for subject $i$ at wave $j,X_{ij}$ is a $p\times 1$ vector of covariates, $Y_i={(Y_{i1}\mathrm{,}{Y}_{i2}\mathrm{,}\dots \mathrm{,}{Y}_{iJ})}^{\text{'}},$  $X_i=(X_{i1}\mathrm{,}{X}_{i2}\mathrm{,}\dots \mathrm{,}{X}_{iJ})$ is a $p\times J$ matrix; $g(\cdot )$ is a monotonic one-to-one differentiable "link function�; $\nu (\cdot )$ is the "variance function� with known form; and $\varphi >0$ is the "dispersion parameter.� Since, in general, the $J\times J$ covariance matrix ${\Sigma }_i$ is hard to specify, we model it as $\text{Cov}\left[Y_i|X_i\right]=V_i=A_i^{1/2}R(\alpha )A_i^{1/2},$ a "working� covariance matrix; where $A_i=\text{diag}\left[\varphi \nu ({\mu }_{i1}),\varphi \nu ({\mu }_{i2})\mathrm{,}\dots \mathrm{,}\varphi \nu ({\mu }_{iJ})\right]$ and $R(\alpha )$ is a "working� correlation matrix, both of dimension $J\times J,$ and $\alpha$ is a vector that fully characterizes $R(\alpha )$ (see Liang and Zeger 1986).

To estimate $\beta$ we select a (single-cohort) sample of $n$ elements from model $\xi$ and we (intend to) measure each of them at $J$ occasions. If all the elements in the sample respond at every single occasion $j,$ the task can be completed with the usual generalized estimating equation (GEE) methodology of Liang and Zeger (1986). However, in any study it is rarely the case that all subjects do respond at all waves. It is more common to have some elements in the sample who drop out of the study.

Under this situation, and assuming that the missing responses can be regarded as missing at random or MAR (see Rubin 1976), in particular that the dropout at a given wave does not depend on the current (unobserved) value, Robins, Rotnitzky and Zhao (1995) proposed to estimate $\beta$ by solving the estimating equations: ${\displaystyle {\sum }_{i=1}^n}(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}{\widehat{\Delta }}_i(y_i-{\mu }_i)=0,$ where ${\mu }_i={({\mu }_{i1}\mathrm{,}{\mu }_{i2}\mathrm{,}\dots \mathrm{,}{\mu }_{iJ})}^{\text{'}},$ ${\widehat{\Delta }}_i=\text{diag}\left[R_{i1}{\widehat{q}}_{i1}^{-1}\mathrm{,}{R}_{i2}{\widehat{q}}_{i2}^{-1}\mathrm{,}\dots \mathrm{,}{R}_{iJ}{\widehat{q}}_{iJ}^{-1}\right],R_{ij}$ is the response indicator for subject $i$ at wave $j,$ and ${\widehat{q}}_{ij}$ is an estimate of the probability that subject $i$ is observed through wave $j.$

For survey applications, one would use the estimating equation ${\displaystyle {\sum }_{i\in s}}\left[w_i(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}{\widehat{\Delta }}_i(y_i-{\mu }_i)\right]=0,$ where $w_i$ is the survey weight for subject $i.$ Another way of writing this equation is ${\displaystyle {\sum }_{i\in s}}(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}{\widehat{\Delta }}_{wi}(y_i-{\mu }_i)=0,$ with ${\widehat{\Delta }}_{wi}\text{=diag}\left[w_i{R}_{i1}{\widehat{q}}_{i1}^{-1}\mathrm{,}{w}_i{R}_{i2}{\widehat{q}}_{i2}^{-1}\mathrm{,}\dots \mathrm{,}{w}_i{R}_{iJ}{\widehat{q}}_{iJ}^{-1}\right]\cdot$

We notice that the diagonal elements of ${\widehat{\Delta }}_{wi}$ are simply wave-specific nonresponse-adjusted survey weights whenever the subject is observed, and are equal to zero whenever the subject is missing. This feature in and of itself suggests a solution to the multi-cohort problem, which will be presented in the next section.

3.2 A novel approach to combining cohorts in longitudinal surveys

Based on the discussion in the previous section, if we have a fixed-panel, fixed-panel-plus-'births', repeated-panel, rotating-panel, split-panel, or refreshment sample survey, we propose to estimate the superpopulation parameter $\beta$ in model $\xi$ by the solution to the estimating equations:

${\Psi }_s(\beta )={\displaystyle \sum _{i\in s}}\frac{\partial {{\mu }^{\prime }}_i}{\partial \beta }{V}_i^{-1}{W}_i(y_i-{\mu }_i)=0;$(3.2)

where the sum is over the sample $s,$ i.e., over all the elements selected (for the first time) in any of the samples $s_{1(1)},s_{2(2)},\dots ,s_{J(J)}.$ The diagonal matrix $W_i$ is $W_i=\text{diag}\left[I_i(U_1)w_{i1}\mathrm{,}{I}_i(U_2)w_{i2}\mathrm{,}\dots ,I_i(U_J)w_{iJ}\right],$   with $w_{ij}$ being the (nonresponse-adjusted) cross-sectional weight for subject $i$ at wave $j$ (as long as subject $i$ is part of sample $s_j$) and $I_i(U_j)$ is the indicator of whether subject $i$ belongs to finite population $U_j$ or not. In Section 3.2.1 we argue why this is a reasonable estimation procedure, and in Section 3.2.2 we discuss the missing value issue.

The cross-sectional weights $w_{ij},$ in $W_i,$ are such that the sample $s_j$ represents $U_j,$ when used in conjunction with said weights. This means that, for each observation $i$ in sample $s_j,$ there has to be a survey weight $w_{ij},$ which could be regarded as the number of units that such observation represents in $U_j.$ However, remember that the sample $s_j$ is composed of different sets of subjects, or different subsamples (the different cohorts), and the integration of these subsamples into a single cross-sectional weight variable $w_{ij}$ may not be a straightforward task.

For the SDR, the construction of the cross-sectional weight for wave $j$ is not too complicated as the different cohorts are selected independently, from non-overlapping populations. The base weight in that case is easy to compute, and all that remains is the adjustment for things like attrition and calibration to known totals in the population $U_j.$

On the other hand, in other situations, for example, when a frame of new members does not exist, the new cohort may need to be selected from the overall population at the given wave, or from a frame containing new members plus some old members, or from multiple frames. In such cases, the building of the cross-sectional weights may not be as straightforward, and the theory of multiple frames may need to be used. We refer the reader to the works of Lohr (2007) and Rao and Wu (2010), and references therein, for cases like that.

Expression (3.2) is a generalization of equation (2.25) in Vieira (2009). The latter is applicable only when all the subjects have the same number of observations or any missing responses can be regarded as missing completely at random or MCAR (see Rubin 1976). As discussed in Robins, et al. (1995), using such an equation when the missing responses are not MCAR produces inconsistent estimators; therefore, with a rotation scheme like that of the SDR, where not all subjects are dropped (or kept) with the same probabilities, its usage would not be appropriate. The adequacy of equation (3.2) in that case and when there are missing responses is addressed in sections 3.2.1 and 3.2.2, respectively. If all subjects have cross-sectional weights that do not vary over time (or have a single longitudinal weight) equation (3.2) reduces to equation (2.25) in Vieira (2009).

3.2.1 Unbiasedness

The unbiasedness property of the estimating function is important because, as Song (2007, Section 5.4) argues, it is the most crucial assumption in order to obtain a consistent estimator.

Let us define ${\beta }_N,$ the so-called "census estimator,� to be the solution to the following finite population estimating equation:

${\Psi }_U({\beta }_N)={\displaystyle \sum _{i\in U}}\frac{\partial {{\mu }^{\prime }}_i}{\partial {\beta }_N}{V}_i^{-1}{\text{I}}_i(\text{U})(y_i-{\mu }_i({\beta }_N))=0\mathrm{,}$(3.3)

where the sum is over $U,$ i.e., over all the elements who became members of the target population in any of $U_{1(1)},U_{2(2)},\dots ,U_{J(J)},$ and ${\text{I}}_i(\text{U})\text{=diag}\left[I_i(U_1)\mathrm{,}{I}_i(U_2)\mathrm{,}\dots I_i(U_J)\right].$ In order to show design-unbiasedness of the estimating function ${\Psi }_s(\beta ),$ we need to show that its design expectation is ${\Psi }_U(\beta )$ for any $\beta .$

The sampling design characteristics of a longitudinal survey can be thought of as those of a multiphase sample, as can be seen in Särndal, Swensson and Wretman (1992, Section 9.9). We therefore use the methodology of multiphase sampling for the derivations. We assume, without loss of generality, that there are only three waves; the derivations with just three waves show the patterns for general $J,$ with respect to unbiasedness and variance.

As we mentioned earlier, we assume that $w_{ij}$ is the cross-sectional weight for subject $i$ at wave $j,$ if that subject belongs to $s_j,$ and zero otherwise. From the theory of multiphase sampling we have that for $i\in s_{1(1)},w_{i1}={\pi }_{i1}^{-1},w_{i2}={\pi }_{i1}^{-1}{\pi }_{i2|s_{1(1)}}^{-1},$ and $w_{i3}={\pi }_{i1}^{-1}{\pi }_{i2|s_{1(1)}}^{-1}{\pi }_{i3|s_{2(1)}}^{-1};$ for $i\in s_{2(2)},w_{i2}={\pi }_{i2}^{-1}$ and $w_{i3}={\pi }_{i2}^{-1}{\pi }_{i3|s_{2(2)}}^{-1};$ and for $i\in s_{3(3)},w_{i3}={\pi }_{i3}^{-1};$ where ${\pi }_{ij}$ is the inclusion probability of subject $i$ in sample $s_{j(j)}$ and ${\pi }_{ij|s_{j-1(j^{\prime })}}$ is the conditional inclusion probability of subject $i$ in sample $s_{j(j^{\prime })}$ given $s_{j-1(j^{\prime })}.$

Using $E_p(\cdot )$ to denote the expectation with respect to the sampling design, we have:

$E_p\left[{\displaystyle \sum _{i\in s}}\frac{\partial {{\mu }^{\prime }}_i}{\partial \beta }{V}_i^{-1}{W}_i(y_i-{\mu }_i)\right]=E_p\left[{\displaystyle \sum _{j=1}^3}{\displaystyle \sum _{i\in s_{j(j)}}}{B}_i{W}_i{e}_i\right];$(3.4)

where $B_i=(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}$ and $e_i=y_i-{\mu }_i.$ For example, for ${\displaystyle {\sum }_{i\in s_{2(2)}}}{B}_i{W}_i{e}_i$ we obtain:

$\begin{array}{c}{E}_p\left[{\displaystyle \sum _{i\in s_{2(2)}}}{B}_i{W}_i{e}_i\right]=E\left\{E\left[{\displaystyle \sum _{i\in U_{2(2)}}}{B}_i{D}_i{e}_i|s_{2(2)}\right]\right\}=E\left\{{\displaystyle \sum _{i\in U_{2(2)}}}{B}_i{D}_i^{*}{e}_i\right\}\\ ={\displaystyle \sum _{i\in U_{2(2)}}}{B}_i{D}_i^{**}{e}_i\stackrel{\text{def}}{=}{\displaystyle \sum _{i\in U_{2(2)}}}{B}_i{\text{I}}_i(\text{U})e_i\mathrm{,}\end{array}$(3.5)

where $D_i=\text{diag}\left[\mathrm{0,}{I}_i(U_2)w_{i2}{I}_i(s_{2(2)})\mathrm{,}{I}_i(U_3)w_{i3}{I}_i(s_{3(2)})I_i(s_{2(2)})\right],$ 
$D_i^{*}=\text{diag}\left[\mathrm{0,}(I_i(U_2)w_{i2}\times I_i(s_{2(2)}))\mathrm{,}(I_i(U_3){\pi }_{i3|s_{2(2)}}{I}_i(s_{2(2)}))/({\pi }_{i2}{\pi }_{i3|s_{2(2)}})\right],$ and $D_i^{**}=\text{diag}\left[\mathrm{0,}(I_i(U_2){\pi }_{i2})/{\pi }_{i2}\mathrm{,}(I_i(U_3){\pi }_{i2})/{\pi }_{i2}\right];$   similarly we can show that $E_p\left[{\displaystyle {\sum }_{i\in s_{1(1)}}}{B}_i{W}_i{e}_i\right]={\displaystyle {\sum }_{i\in U_{1(1)}}}{B}_i{\text{I}}_i(\text{U})e_i$ and $E_p\left[{\displaystyle {\sum }_{i\in s_{3(3)}}}{B}_i{W}_i{e}_i\right]={\displaystyle {\sum }_{i\in U_{3(3)}}}{B}_i{\text{I}}_i(\text{U})e_i.$ From these expressions and equation (3.4) we conclude that $E_p\left[{\Psi }_s(\beta )\right]={\Psi }_U(\beta )$ for any $\beta ,$ which means that the estimating function ${\Psi }_s(\beta )$ is design-unbiased for the finite population estimating function.

Furthermore, as the target of inference is the superpopulation parameter, we need to guarantee that the model for ${\mu }_{ij}$ is such that $E_{\xi }(Y_{ij}-{\mu }_{ij})=0$ is satisfied, where $E_{\xi }(\cdot )$ represents the expectation with respect to model $\xi .$ For if this is the case, we have: 

$E_{\xi p}\left[{\Psi }_s(\beta )\right]\stackrel{\text{def}}{=}{E}_{\xi }{E}_p\left[{\Psi }_s(\beta )\right]=E_{\xi }\left[{\Psi }_U(\beta )\right]={\displaystyle \sum _{i\in U}}\frac{\partial {{\mu }^{\prime }}_i}{\partial \beta }{V}_i^{-1}{\text{I}}_i(\text{U})E_{\xi }(y_i-{\mu }_i)=0;$

so that the estimating function ${\Psi }_s(\beta )$ is model-design unbiased. The requirement $E_{\xi }(Y_{ij}-{\mu }_{ij})=0$ means that the mean model needs to be correctly specified; consequently, one needs to pay attention to residual diagnostics for the particular model being fitted.

3.2.2 A note on nonresponse

In the SDR, as in any other (longitudinal) survey, there is nonresponse. Some sampled individuals choose not to participate at all, whereas some subjects participate in some waves but not in others. The SDR remedies this situation by making a nonresponse adjustment to the cross-sectional survey weights.

Assume that the nonresponse adjustment at wave $j$ is a multiplication by the inverse of the estimated wave $j$ response probability ${\widehat{\pi }}_{rij}.$ For example, the nonresponse-adjusted weight for a person who did respond at wave 3 (and was first selected at wave 2), i.e., for $i\in r_{3(2)},$ would be $w_{ri3}={\pi }_{i2}^{-1}{\pi }_{i3|s_{2(2)}}^{-1}{\widehat{\pi }}_{ri3}^{-1}.$

We need to redefine the estimating equation, to include only the respondents, as ${\Psi }_r(\beta )={\displaystyle {\sum }_{i\in r}}(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}{W}_{ri}(y_i-{\mu }_i)=0,$ where the sum is over the respondent set $r,$ i.e., over all the elements who belonged for the first time in any of the respondent sets $r_{1(1)},r_{2(2)},\dots ,r_{J(J)},$ and the matrix $W_{ri}$ is $W_{ri}=\text{diag}\left[I_i(U_1)w_{ri1}\mathrm{,}{I}_i(U_2)w_{ri2}\mathrm{,}\dots ,I_i(U_J)w_{riJ}\right].$ Also, denote by $r_{j(j^{\prime })}$ the set of cohort $j^{\prime }$ respondents at wave $j.$ Obviously, $w_{rij}=0$ if $i\notin r_j={\displaystyle {\cup }_{j^{\prime }=1}^j{r}_{j(j^{\prime })}}.$

If additionally, the response mechanism $(R)$ can be assumed to be MAR, we then have, for example for ${\displaystyle {\sum }_{i\in r_{2(2)}}}{B}_i{W}_{ri}{e}_i:$

$E_R\left\{{\displaystyle \sum _{i\in r_{2(2)}}}{B}_i{W}_{ri}{e}_i\right\}=E_R\left\{{\displaystyle \sum _{i\in s_{2(2)}}}{B}_i{D}_i{e}_i\right\}={\displaystyle \sum _{i\in s_{2(2)}}}{B}_i{D}_i^{*}{e}_i={\displaystyle \sum _{i\in s_{2(2)}}}{B}_i{D}_i^{**}{e}_i\stackrel{\text{def}}{=}{\displaystyle \sum _{i\in s_{2(2)}}}{B}_i{W}_i{e}_i\mathrm{,}$(3.6)

where $D_i=\text{diag}\left[\mathrm{0,}{I}_i(U_2)w_{ri2}{I}_i(r_{2(2)})\mathrm{,}{I}_i(U_3)w_{ri3}{I}_i(r_{3(2)})\right],$
 $D_i^{*}\text{ = diag}\left[\mathrm{0,}(I_i(U_2){\pi }_{ri2})/({\pi }_{i2}\times {\widehat{\pi }}_{ri2})\mathrm{,}(I_i(U_3){\pi }_{ri3})/({\pi }_{i2}{\pi }_{i3|s_{2(2)}}{\widehat{\pi }}_{ri3})\right],$ and $D_i^{**}=\text{diag}\left[\mathrm{0,}{I}_i(U_2)w_{i2}\mathrm{,}{I}_i(U_3)w_{i3}\right].$ The third equality in (3.6) requires that the nonresponse model used for ${\widehat{\pi }}_{rij}$ satisfies $E_R\left[I_i(r_{j(j^{\prime })})\right]\stackrel{\text{def}}{=}{\pi }_{rij}={\widehat{\pi }}_{rij}.$ This means that in the model for ${\widehat{\pi }}_{rij}$ we have to include as much available information, thought to influence the nonresponse propensity, as possible, in order for this assumption (i.e., the MAR assumption) to be tenable. For example, if the nonresponse is thought to be independent across waves, one should include, in the model for ${\widehat{\pi }}_{rij},$ as many variables from the corresponding wave as possible. If, on the other hand, it is reasonable to assume that the response propensity at a given wave depends on previous responses (and possibly response history), then those responses should be included in the response model, and so on.

The design as well as the model-design unbiasedness follow immediately from (3.6) together with the previous section. Hereafter we therefore ignore the issue of nonresponse for notational simplicity.

3.3 Variance and variance estimation

We now develop a (Taylor Series) linearization for the variance of the proposed estimator. The basic technique is due to Binder (1983). For simplicity in the derivations and notation we divide through by $N;$ we redefine

${\Psi }_s(\beta )=N^{-1}{\displaystyle \sum _{i\in s}}\frac{\partial {{\mu }^{\prime }}_i}{\partial \beta }{V}_i^{-1}{W}_i(y_i-{\mu }_i)$ and ${\Psi }_U(\beta )=N^{-1}{\displaystyle \sum _{i\in U}}\frac{\partial {{\mu }^{\prime }}_i}{\partial \beta }{V}_i^{-1}{\text{I}}_i(\text{U})(y_i-{\mu }_i),$

where $N={\displaystyle {\sum }_{j=1}^J{N}_j}.$ Let $\widehat{\beta }$ be our estimator, which satisfies ${\Psi }_s(\widehat{\beta })=0,$ and let ${\beta }_N$ be the "census estimator,� which satisfies ${\Psi }_U({\beta }_N)=0.$ Assume ${\beta }_N-\beta =O_P(1/\sqrt{N_m})$ and $\widehat{\beta }-{\beta }_N=O_P(1/\sqrt{n_m}),$ with $N_m=\min \left\{N_1\mathrm{,}{N}_2\mathrm{,}\dots \mathrm{,}{N}_J\right\}$ and $n_m=\min \left\{n_1\mathrm{,}{n}_2\mathrm{,}\dots \mathrm{,}{n}_J\right\}.$ We can write the total error of $\widehat{\beta }$ as $\widehat{\beta }-\beta =(\widehat{\beta }-{\beta }_N)+({\beta }_N-\beta )=$ Sampling Error + Model Error. After some straightforward calculations, the total variance, or more precisely the total MSE, can be decomposed as:

$V_{\text{Tot}}=E_{\xi p}(\widehat{\beta }-\beta ){(\widehat{\beta }-\beta )}^{\text{'}}=V_{\text{Sam}}+2\otimes C_{\text{Sam}-\text{Mod}}+o(1/n_m)\mathrm{,}$(3.7)

where $2\otimes A=A+A^{\prime }$ for any matrix $A,V_{\text{Sam}}=E_{\xi }{V}_p$ is the "sampling variance� component, $2\otimes C_{\text{Sam}-\text{Mod}}$ is the cross "sampling-model variance� component, $V_p=$ $E_p\left[(\widehat{\beta }-{\beta }_N){(\widehat{\beta }-{\beta }_N)}^{\text{'}}\right],$ $C_{\text{Sam}-\text{Mod}}=E_p{C}_{\xi },$ and $C_{\xi }=E_{\xi }(\widehat{\beta }-\beta ){({\beta }_N-\beta )}^{\text{'}}.$ Furthermore, by Taylor series expansions we can obtain the following approximations: $\widehat{\beta }-{\beta }_N={\left[H({\beta }_N)\right]}^{-1}{\Psi }_s({\beta }_N)+o_P(1/\sqrt{n_m}),$ $\widehat{\beta }-\beta ={\left[\widehat{H}(\beta )\right]}^{-1}{\Psi }_s(\beta )+o_P(1/\sqrt{n_m}),$ and ${\beta }_N-\beta ={\left[H(\beta )\right]}^{-1}{\Psi }_U(\beta )+o_P(1/\sqrt{N_m}),$ where we define $H(\beta )=N^{-1}{\displaystyle {\sum }_{i\in U}}(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}{\text{I}}_i(\text{U})(\partial {\mu }_i/\partial \beta )$ and $\widehat{H}(\beta )=N^{-1}{\displaystyle {\sum }_{i\in s}}(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}{W}_i(\partial {\mu }_i/\partial \beta ).$

We then get, for $V_p$ and $C_{\xi }$ in (3.7),

$V_p={\left[H({\beta }_N)\right]}^{-1}{\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right]{\left[H({\beta }_N)\right]}^{-1}+o_P(1/n_m),$(3.8)

$\begin{array}{c}{C}_{\xi }={\left[\widehat{H}(\beta )\right]}^{-1}{E}_{\xi }\left[{\Psi }_s(\beta ){{\Psi }^{\prime }}_U(\beta )\right]{\left[H(\beta )\right]}^{-1}+o_P(1/n_m)\\ =N^{-1}{\left[\widehat{H}(\beta )\right]}^{-1}{\widehat{H}}_{\Sigma V}(\beta ){\left[H(\beta )\right]}^{-1}+o_P(1/n_m),\end{array}$(3.9)

where ${\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right]=E_p\left[{\Psi }_s({\beta }_N){{\Psi }^{\prime }}_s({\beta }_N)\right]$ and ${\widehat{H}}_{\Sigma V}(\beta )=N^{-1}{\displaystyle {\sum }_{i\in s}}\left[(\partial {{\mu }^{\prime }}_i/\partial \beta )V_i^{-1}{W}_i{\Sigma }_i\times V_i^{-1}(\partial {\mu }_i/\partial \beta )\right];$ the derivation of (3.9) can be found in the Appendix.

In conclusion, so far we have found that:

$\begin{array}{cc}\hfill V_{\text{Tot}}& =E_{\xi }{V}_p+2\otimes E_p{C}_{\xi }+o(1/n_m)\hfill \\ & =E_{\xi }\left\{{\left[H({\beta }_N)\right]}^{-1}{\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right]{\left[H({\beta }_N)\right]}^{-1}\right\}\hfill \\ & +2\otimes N^{-1}{E}_p\left\{{\left[\widehat{H}(\beta )\right]}^{-1}{\widehat{H}}_{\Sigma V}(\beta ){\left[H(\beta )\right]}^{-1}\right\}+o(1/n_m).\hfill \end{array}$(3.10)

In (3.10) all the terms can be estimated by "plugging in� the estimate $\widehat{\beta },$ except for the term ${\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right];$ this is the subject of the next section.

If the sampling fraction is small, i.e., $n\ll N,$ the first term in expression (3.10) is a good approximation for the total variance; i.e., the expression for $V_{\text{Tot}}$ is simply $E_{\xi }{V}_p$ (and lower order terms). If, on the other hand, the sampling fraction is large, both terms in (3.10) are required.

3.3.1 Design variance of the estimating function

In order to derive an expression for ${\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right],$ we assume $J=3,$ as before. The methodology is that of two-phase sampling (more precisely, multiphase sampling), as discussed in chapter 9 of Särndal, et al. (1992). After some derivations (see Appendix), and defining $B_i={(\partial {{\mu }^{\prime }}_i/\partial \beta )|}_{\beta ={\beta }_N}{V}_i^{-1},e_i=y_i-{\mu }_i({\beta }_N),e_{i(1\cdots 3)}=e_i,e_{i(2\cdots 3)}={(\mathrm{0,}{e}_{i2}\mathrm{,}{e}_{i3})}^{\prime },$ and $e_{i(3\cdots 3)}={(\mathrm{0,0,}{e}_{i3})}^{\prime },$ we obtain:

${\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right]={\displaystyle \sum _{j=1}^3{D}_{(j)}}={\displaystyle \sum _{j=1}^3}{\displaystyle \sum _{k=j}^3{D}_{(j)k}}\mathrm{,}$(3.11)

where $D_{(j)}\stackrel{\text{def}}{=}{N}^{-2}{\text{Var}}_p({\displaystyle {\sum }_{i\in s_{j(j)}}{B}_i{W}_i{e}_i})={\displaystyle {\sum }_{k=j}^3{D}_{(j)k}},$ for $j=\mathrm{1,2,3},$

$N^2{D}_{(j)j}\stackrel{\text{def}}{=}\text{Var}\left[{\displaystyle \sum _{i\in s_{j(j)}}}{w}_{ij}{B}_i{\text{I}}_i(\text{U})e_{i(j\dots 3)}\right]\mathrm{,}\text{for}j=\mathrm{1,2,3,}$

$N^2{D}_{(j-1)j}\stackrel{\text{def}}{=}E\left\{\text{Var}\left[{\displaystyle \sum _{i\in s_{j(j-1)}}}{w}_{ij}{B}_i{\text{I}}_i(\text{U})e_{i(j\dots 3)}|s_{j-1(j-1)}\right]\right\}\mathrm{,}\text{for}j=\mathrm{2,3,}$

$N^2{D}_{(1)3}\stackrel{\text{def}}{=}E\left\{E\left[\text{Var}({\displaystyle \sum _{i\in s_{3(1)}}}{w}_{i3}{B}_i{\text{I}}_i(\text{U})e_{i(3\dots 3)}|s_{2(1)}\mathrm{,}{s}_{1(1)})|s_{1(1)}\right]\right\}\mathrm{,}$

and in the Appendix we show that:

$N^2{D}_{(j)k}=\text{Var}\left[{\displaystyle \sum _{i\in s_{k(j)}}}{w}_{ik}{B}_i{\text{I}}_i(\text{U})e_{i(k\dots 3)}\right]-\text{Var}\left[{\displaystyle \sum _{i\in s_{k-1(j)}}}{w}_{i\mathrm{,}k-1}{B}_i{\text{I}}_i(\text{U})e_{i(k\dots 3)}\right]\mathrm{,}$

for $j=\mathrm{1,2,3},$ and $3\ge k>j.$ In general, we have proved the following

Property 3.1  The (design) variance of ${\Psi }_s({\beta }_N)$ can be decomposed as:

               $\begin{array}{l}{\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right]\\ =\frac{1}{N^2}{\displaystyle \sum _{j^{\prime }=1}^J}{\displaystyle \sum _{j=j^{\prime }}^J}\left\{{\text{Var}}_p\left[{\displaystyle \sum _{i\in s_{j(j^{\prime })}}}{w}_{ij}{B}_i{\text{I}}_i(\text{U})e_{i(j\dots J)}\right]-{\text{Var}}_p\left[{\displaystyle \sum _{i\in s_{j-1(j^{\prime })}}}{w}_{i\mathrm{,}j-1}{B}_i{\text{I}}_i(U)e_{i(j\dots J)}\right]\right\}\end{array}$(3.12)

               $=\frac{1}{N^2}{\displaystyle \sum _{j=1}^J}\left\{{\text{Var}}_p\left[{\displaystyle \sum _{i\in s_j}}{w}_{ij}{B}_i{\text{I}}_i(U)e_{i(j\dots J)}\right]-{\text{Var}}_p\left[{\displaystyle \sum _{i\in s_{j-1}}}{w}_{i\mathrm{,}j-1}{B}_i{\text{I}}_i(U)e_{i(j\dots J)}\right]\right\}\mathrm{,}$(3.13)

where we let $w_{i\mathrm{,}j-1}=0$ whenever $j=j^{\prime },w_{i0}=0,$ and to get (3.13) we have changed variables and used the independence among cohorts.

In (3.11), (3.12), and (3.13) we have assumed that the cohorts are design-independent. However, in some cases this assumption may not be tenable; an example of such a case is the multiple frame situation discussed in the first part of Section 3.2. Another instance in which it may not be appropriate to assume cohort independence is when weight adjustments cross cohorts, which is the case of the SDR; we discuss this issue in Section 5. Calculations for the case of three cohorts, in the Appendix, show that (3.13) holds for the variance terms even without independence. The Appendix also identifies conditions under which it is a good approximation for the covariance terms.

3.3.2 Estimation

The estimation of $V_{\text{Tot}}$ in (3.10) can be achieved as follows. $H({\beta }_N),\widehat{H}(\beta ),$ and $H(\beta )$ can be estimated by $\widehat{H}(\widehat{\beta }).$ ${\widehat{H}}_{\Sigma V}(\beta )$ can be estimated by ${\widehat{H}}_{\Sigma V}(\widehat{\beta }),$ where ${\Sigma }_i=\text{Cov}\left[Y_i|X_i\right]$ can be estimated by ${\widehat{e}}_i{{\widehat{e}}^{\prime }}_i.$

We use (3.13) in Property 3.1 to estimate ${\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right].$ As long as there is a method to estimate the variance of (cross-sectional) Horvitz-Thompson (H-T) estimators, expression (3.13) can be used. If we define $Z_{ij}=B_i{\text{I}}_i(\text{U})e_{i(j\dots J)},$ we notice that each of the terms involved in the computation of (3.13), terms like ${\text{Var}}_p\left[{\displaystyle {\sum }_{i\in s_j}}{w}_{ij}{Z}_{ij}\right],$ is simply the variance of a wave $­j$ H-T estimator. Obviously, the variance estimation method needs to account for the sampling design as well as for any nonresponse and calibration adjustments performed, but this does not present any additional complications beyond what is found in any cross-sectional problem, as everything is implemented cross-sectionally. The SDR uses replication to estimate variances of cross-sectional estimators, but any method of design variance estimation can be used.

We use the cross-sectional replicate weights that SDR provides, but we do not re-estimate the parameter of interest at each replicate. First, note that we require replication only for the estimation of the "meat� $({\text{Var}}_p\left[{\Psi }_s({\beta }_N)\right])$ of the design variance $(E_{\xi }{V}_p).$ Secondly, although $\widehat{\beta }$ does appear in the expression for the H-T estimator whose variance needs to be calculated (and re-calculated at each replicate), the work of Roberts, Binder, Kova�ević, Pantel and Phillips (2003), who apply the "estimating function bootstrap� (Hu and Kalbfleisch 2000) to survey data, show that in a setting like ours, it is not necessary to re-compute the estimator at each replicate, but that the full-sample estimator suffices. This simplification speeds up the computation of the replicate estimates.

As a way of illustration, say we currently are at wave $j,$ i.e., we are estimating the $j^{\text{th}}$ term in (3.13). The $r^{\text{th}}$ replicate of the first term is ${\displaystyle {\sum }_{i\in s_j}}{w}_{ij}^{(r)}{B}_i(\widehat{\beta }){\text{I}}_i(\text{U})e_{i(j\dots J)}(\widehat{\beta }),$ where $w_{ij}^{(r)}$ is the $r^{\text{th}}$ replicate weight for subject $i$ at wave $j,$ and the $r^{\text{th}}$ replicate of the second term is ${\displaystyle {\sum }_{i\in s_{j-1}}}{w}_{i,j-1}^{(r)}{B}_i(\widehat{\beta }){\text{I}}_i(\text{U})e_{i(j\dots J)}(\widehat{\beta }),$ where $w_{i\mathrm{,}j-1}^{(r)}$ is the $r^{\text{th}}$ replicate weight for subject $i$ at wave $j-1.$

Previous | Next