Linearization Variance Estimators for Longitudinal Survey Data 

A. Demnati and J. N. K. Rao 

A. Demnati, Social Survey Methods Division, 

Statistics Canada, Ottawa, Ontario, Canada, K1A 0T6 

Abdellatif.Demnati@statcan.ca 

J. N. K. Rao, School of Mathematics and Statistics, 

Carleton University, Ottawa, Ontario, Canada, K1S 5B6 

JRao@math.carleton.ca 

Summary 

In survey sampling, Taylor linearization is often used to obtain variance estimators for nonlinear finite population 

parameters, such as ratios, regression and correlation coefficients, which can be expressed as smooth functions of totals. 

Taylor linearization is generally applicable to any sampling design, but it can lead to multiple variance estimators that are 

asymptotically design unbiased under repeated sampling. The choice among the variance estimators requires other 

considerations such as (i) approximate unbiasedness for the model variance of the estimator under an assumed model, (ii) 

validity under a conditional repeated sampling framework. Demnati and Rao (2001) proposed a new approach to deriving 

Taylor linearization variance estimators that leads directly to a unique variance estimator that satisfies the above 

considerations. In this paper, we extended the work of Demnati and Rao (2001) to deal with longitudinal surveys which 

lead to dependent observations and to multiple weights on the same unit. We consider a variety of longitudinal sampling 

designs, covering panel surveys, household panel surveys as well as rotating surveys. 

Key Words: Multiple weights; Repeated surveys; Taylor linearization. 

1. Introduction 

Taylor linearization is a popular method of variance estimation for complex statistics such as ratio and regression 

estimators and logistic regression coefficient estimators. It is generally applicable to any sampling design that permits 

unbiased variance estimation for linear estimators, and it is computationally simpler than a resampling method such as the 

jackknife. However, it can lead to multiple variance estimators that are asymptotically design unbiased under repeated 

sampling. The choice among the variance estimators, therefore, requires other considerations such as (i) approximate 

unbiasedness for the model variance of the estimator under an assumed model, (ii) validity under a conditional repeated 

sampling framework. For example, in the context of simple random sampling and the ratio estimator, Y

ˆ = ( y / x)X , of the R 

population total Y , Royall and Cumberland (1981) showed that a commonly used linearization variance estimator 

2 θ = N 2

(n

−1 

− N −1)s does not track the conditional variance of Y

ˆ

 given x , unlike the jackknife variance estimator L z R

2 θ J . Here y and x are the sample means, X is the known population total of an auxiliary variable x , s is the sample z

variance of the residuals zi = yi − ( y / x)xi and (n, N) denote the sample and population sizes. By linearizing the 

2

jackknife variance estimator,θ J , we obtain a different linearization variance estimator, θ = (X / x) θ L JL , which also

tracks the conditional variance as well as the unconditional variance, where X = X / N is the mean of x . As a result, θ JL 

or θ J

 may be preferred over θL . Valliant (1993) obtained θ JL for the post-stratified estimator and conducted a simulation 

study to demonstrate that both θ J and θ JL possess good conditional properties given the estimated post-strata counts. 

Särndal, Swensson and Wretman (1989) showed that θ JL is both asymptotically design unbiased and asymptotically model 

unbiased in the sense of Em (θ JL ) = Varm (Y

ˆ

R ) , where Em denotes model expectation and Varm (Y

ˆ

R ) is the model variance 

of Y

ˆ

 under a Αratio model≅: Em ( yi ) = β xi ; i =1,..., N and the yi ’s are independent with model variance R

2 Varm ( yi ) =σ xi , σ

2 > 0 . Thus, θ JL is a good choice from either the design-based or the model-based perspective. 

Demnati and Rao (2001) proposed a new approach to variance estimation that is theoretically justifiable and at the same 

time leads directly to a θ JL -type variance estimator for general designs. Afterward, Demnati and Rao (2002) extended their 

1 

Page 138 
 

 

 

 

 

 

method to the case of missing responses when adjustment for complete nonresponse and imputation based on smooth 

functions of observed values, in particular ratio imputation, are used. 

While several methods have been proposed to correctly estimate the variance of an estimator under crosssectional data, methods of variance estimation from longitudinal data are limited, despite an increase in recent years in the 

number of longitudinal surveys. The dynamic nature of longitudinal populations such as births and deaths, the limitation of 

longitudinal samples to take into account such features of longitudinal populations, and multiple weights are one of the 

difficulties inherent to longitudinal data in addition to the customary cross-sectional difficulties such as complex designs, 

missing data and so on. This paper is a first attempt to extend the Demnati and Rao (2001,2002) method to the case of 

longitudinal survey data. Section 2 gives a brief account of the method for the case of cross-sectional data, and section 3 

presents the extension to longitudinal survey data. 

2. Cross-Sectional Data 

We give a brief account of the Demnati and Rao (2001) method under full response. Suppose an estimator θ

ˆ

 of a 

parameter θ can be expressed as a differentiable function g(Yˆ

) of estimated totals Yˆ = (Y

ˆ

1

,...,Y

ˆ 

m )

T , where 

Y

ˆ

j = ∑i∈U di (s) yij is an estimator of the population Y j , j =1,.., m , where di (s) = 0 if the unit i is not in the sample s , 

U is the set of population units, and θ = g(Y ) with Y = (Y ,...,Y )

T . We may write θ

ˆ

 as θ

ˆ = f (d(s), Ay 1 m ) and

j

th T θ = f (1, Ay ) , where Ay

 is an m × N matrix with column y j = ( y1 j ,..., ymj ) , j =1,..., N , d (s) = (d1(s),..., d N (s))T 

and 1 is the N -vector of 1's. For example, if θ

ˆ

 denotes the ratio estimator Y

ˆ = X (∑i∈U di (s) yi ) /(∑i∈U di (s)xi R ) , then

m = 2 , y1i = yi , y2i = xi and f (1, Ay ) reduces to the total Y , noting that X (Y / X ) = Y . Note that Y

ˆ

R

 is a function of 

d(s) , y , x and the known total X , but we dropped X for simplicity and write Y

ˆ = f (d(s), y , x) . If the Horvitz- R 

Thompson weights are used, then di (s) =1/π i for i ∈s , where π i

 is the probability of selecting unit i in the sample s . 

Let f (b, Ay ) = f (b) for arbitrary real numbers b = (b ,...,bN )

T . Demnati and Rao (2001) showed that the Taylor 1 

linearization of θ

ˆ −θ , namely 

T (Yˆ θ −Y ),

ˆ −θ = g(Yˆ

) − g(Y ) ≈ (∂g(a) / ∂a)

a=Y 

is equivalent to 

θ

ˆ −θ ≈ ∑N 

=1 (∂f (b) / ∂bk ) (dk (s) −1) k b=1 (2.1) 

~ 

= z 

T (d(s) − 1), 

)

T ~ ~ ~ ~ where ∂g(a) / ∂a = (∂g(a) / ∂a ,..., ∂g(a) / ∂a and z = (z1

,..., zN )

T

 with zk = ∂f (b) / ∂bk b 1 . It follows from (2.1) th at 1 m =

~ ˆ ~ a variance estimator of θ

ˆ is approximately given by the variance estimator of the estimated total ∑ di (s)zi = Y (z , d (s)) ; 

~ θ

ˆ ≈θ

~ 

that is, var( ) (z , d(s)) , where θ( y, d(s)) denotes the variance estimator of Y

ˆ = Yˆ

( y, d(s)) in operator notation usin g 

~ the vector of design weights d(s) . Now we replace zk

 by zk = ∂f (b) / ∂bk , since zk =s are unknown, to get a b=d(s) 

linearization variance estimator 

θ ˆ

L (θ ) =θ(z, d(s)) . (2.2) 

Note that θL (θ

ˆ

) given by (2.2) is simply obtained from the formula θ( y, d(s)) for Y

ˆ = Yˆ

( y, d(s)) by replacing yi by zi 

~ for i ∈s . Note that we do not first evaluate the partial derivatives ∂f (b) / ∂bk

 at b = 1 to get zk

 and then substitute 

~

estimates for the unknown components of z . Our method, therefore, is similar in spirit to Binder(1996)’s approach. The 

~ variance estimator θ (θ

ˆ

) is valid because zk

 is a consistent estimator of zk . L 

2 

Page 139 
 

 

 

Suppose θ

ˆ

 is the ratio estimator Y

ˆ = X (∑ di (s) yi ) /(∑ di (s)xi ) , where ∑ denotes summation over i ∈U . Then R 

f (b) = X (∑ bi yi ) /(∑ b x i ) = XY ˆ( y, b) / Y

ˆ

i (x, b) and

X 

ˆ

zk = ∂f (b) / ∂b =

Xˆ (yk − Rxk ).

k b d(s) =

For simple random sampling, θL (Y

ˆ

R ) =θ(z, d(s)) agrees with θ JL = (X / x)

2θ L . 

The above derivation is simple and natural. On the other hand, in the standard linearization method, θ

ˆ

 is first 

expressed in terms of elementary components Y

ˆ

1

,...,Y

ˆ as g(Yˆ

) and the partial derivatives ∂g(a) / ∂a j m are then evaluated 

at a =Y . It is interesting to note that all the components Y

ˆ

j

 use the same weights di (s) and our approach always takes 

partial derivatives of f (b) with respect to bk at b = d(s) . It is not necessary to first express θ

ˆ

 in terms of elementary 

components. Demnati and Rao (2001) applied the method to a variety of problems, covering regression calibration 

estimators of a total Y and other estimators defined either explicitly or implicitly as solutions of estimating equations. 

They obtained a new variance estimator for a general class of calibration estimators that includes generalized raking ratio 

and generalized regression estimators. They also extended the method to two-phase sampling and obtained a variance 

estimator that makes fuller use of the first phase sample data compared to traditional linearization variance estimators. 

3. Longitudinal Data 

As in the usual finite population situation, we identify a population of size N during a fixed period of time by a 

set of indices P = {1,..., N} . Let P

(t)

 denote the cross-sectional population at time t and let J (t) denote the cross-sectional i

J (t) (t) (t) (t)

population membership indicator for unit i , i = 1,..., N , i.e. i = 1 if i ∈ P and Ji = 0 if i ∉ P with 

J (t) ) P

(t) 

+ = ∑i Ji

(t) 

= N (t the size of the cross-sectional population , where ∑i denotes summation over the population 

units i . In this section, we consider variance estimation of a differentiable function of estimated cross-sectional totals 

under a series of longitudinal sampling designs assuming full response and no post-stratification. It is shown that a variance 

estimator can be obtained through variance estimation of a design-weighted estimator of total of a synthetic population. 

The case of estimators obtained as solution of estimating equations can be obtained along the lines of Demnati and Rao 

(2001), and the case of estimators obtained from survey data with missing responses can be obtained along the lines of 

Demnati and Rao (2002). In this version, details of these extensions are omitted for space reason. 

To discuss multiple weights inherent in repeated surveys, suppose an estimator θ

ˆ

 of a parameter θ can be 

= (Y

ˆ (1) expressed as a differentiable function θ

ˆ = g(Yˆ

) of estimated cross-sectional totals Yˆ ,...,Y

ˆ(T )

)

T , where Y

ˆ(t)

 is an 

J (t) (t) = (Y (1) estimator of the population total Y (t) = ∑iy at time t , t = 1,...,T , and θ = g(Y ) with Y ,...,Y (T )

)

T . We may i i 

(t ) ~ 

Y

ˆ(t) ˆ(t) (t) (t) (t) T (t) (t) (t) write as Y = ∑G

g ∑i d gi (s g ) yi = 1

G

(t ) A~ y with specified constants G , where 

d 

~ ~ (t) ~ A

(

~ 

t) 

= (d (t)(s 

(t) ),..., d (t)

(s ))T , d (t)(s g 

(t) ) are vectors of cross-sectional weights at time t , 1 is the G (t) -vector of 1 1 (t ) (t ) g G

(t )

d G G 

(t) (t) (t)

1’s and y = ( y ,..., y )

T . In this section we show that the estimator θ

ˆ

 can be written as θ

ˆ = f (d(s)) under a wi de 1 N 

class of longitudinal survey designs. Therefore, it follows from Demnati and Rao (2001) that a linearization varian ce 

estimator of θ

ˆ

 is given by 

θL (θ

ˆ

) =θ(z, d(s)) , (3.1) 

where 

d (s) = vech(Ad

T ) , (3.2) 

= (A

(1)T T (t) (t) (t) (t) (t) ( Ad d ,..., Ad 

(T )T ) , Ad = (d1

(s

1 ),..., d (t )

(s (t ) ))T 

, d g

t)

(s) are vectors of design weights, vech(A) denotes the 

G G 

vector obtained by stacking distinct column of A underneath each other, in order from left to right, and 

3 

Page 140 
 

 

 

 

 

/ ) ∂ I z k = ∂f (b bk . (3.3) b=d (s) 

3.1 Independent Samples 

(t) (t) d (t)

(s 

(t) Under independent samples we have G (t) 

= 1 , Y

ˆ(t) 

= Y

ˆ

( y , d (t)(s 

(t) )) = ∑i

y ) , and i i 

d(s) = vech(d (1) (s 

(1) ),..., d (T )

(s

(T ) )) , (3.4) 

d (t)(s

(t)

) = (d (t) (t) where (s 

(t) ),..., d (t)

(s ))T is the vector design weights at time t . If the Horvitz-Thompson weights are 1 N 

(t)

) = J (t) (t) (t) (t) (t)

used then d (t)(s a (s) / π

(t)

, a = 1 if unit i belongs to the cross-sectional sample s and a = 0 if not, i i i i i i 

(t)

/ J (t)

and π

(t) 

= Pr(i ∈ s = 1) .

i i 

3.2 Panel Sample 

Under panel sample situation, where observations are collected over time on the same sampled units, we have 

G (t) (t), d (0) = 1 , Y (s

ˆ(t) 

= Y

ˆ

( y 

(0) )) and 

d(s) = d (0) (s) , (3.5) 

where d (0) (s

(0) ) , i =1,..., N , are the basic survey weights of the unique sample selected at the initial wave (say at time i 

t = 0 ). 

3.3 Household Panel Sample 

Under household panel sample situation, where observations are collected over time not only on the sampled units 

selected at the initial wave but also on non-sampled units who join households containing at least one sampled unit, we 

(t) ( ~ t) (t have G (t) 

= 1 , Y

ˆ(t) = ∑ y di (s 

)

) , and i 

d(s) = d (0) (s) , (3.6) 

~ J (t) I (t) d (0) (s

(0) )α

(t) α

(t) J (t) I (t) (0) (s

(0) ) , I (t) where di

(t)

(s

(t)

) = ∑ j j ji j ji , α

(

ji

t) 

= i = 1/ ∑k k ki ak ji is an indicator variable indicating 

(0) (0) ) = d (0) (0) )E(a 

(0) (s

(0) if the unit j lives in the same household as unit i at time t , and a )) .

k (sk (sk The value for α

(

ji

t) 

is obtained through the multiplicity approach, Sirken (1970). 

3.4 Multiple Samples 

New entrants who have zero probability to join the initial population units are not covered, unless a supplementary 

sample is taken. That’s why multiple samples are used at the estimation stage at time t under household panel surveys. 

Under multiple samples situation each sample s

(t)

, g = 1,..., G (t) can be viewed as a sample from the same cross-sectional g 

population using a different sampling frame F (t) . The frames are probably incomplete, but the union of the G (t) frames g

covers the entire cross-sectional population. The G (t)

 frames overlap. Let J ~

gi

(t) be the conditional frame membership 

~(t) J (t) ~(t)

indicator variable, i.e., J gi = 1 if i ∈ Fg 

(t)

| i = 1 and J gi = 0 otherwise. When combining the G (t) samples, the 

weights have to be adjusted to account for the multiplicity of the sampling frames: 

~ d gi

(t)

(s) = d (t) s α gi

(t , (3.7) gi ( ) ~ ) 

~ 

α

~(t) J (t) ~ where gi are constants satisfying the constraints ∑ g gi α gi

(t) 

= 1 gi and d (s) are the cross-sectional weights resulting (t) 

α

~(t) ~(t) ~ ~(t)

from frame Fg 

(t) . A first choice for gi is α gi = 1/ J +

(t

i

) , where J +i

 is the multiplicity of unit i at time t , i.e., the total 

α

~(t)

number of sampling frames reporting unit i at time t , (Sirken, 1970). A second choice for gi is 

4 

Page 141 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

α

~(t) ~(t) (t)

/ ~(t) (t) f (t)

gi = J gi f g ∑ g J gi f g

 in which case one wants to take into account the sampling fraction g of each sampling 

frame (Kalton and Anderson, 1986). Under multiple frames, the estimated total Y

ˆ(t)

 can be expressed as before as 

(t) (t)

1

G

T 

(t ) A~ y . Hence, estimation of the variance can be done along the line of sub-section 3.1-3.3. Details are omitted for 

d 

simplicity. 

3.5 Rotating Samples 

Rotating samples are used as a compromise between independent samples and panel surveys. Rotating samples 

control overlap between successive time periods and provide unbiased cross-sectional estimates while taking advantage of 

the correlation between two time periods. For example, in labor-force surveys which are conducted monthly in many 

countries, overlap is controlled by re-interviewing a high fraction of previous selected households (say G −1 rotation 

groups) while new households are selected for a first interview(say one rotation group) for a total of G rotation groups or 

panel samples. Each household remains in the same rotation group for a predetermined number of months. Unlike multiple 

samples, rotating groups do not overlap. A simple example of a rotating group is given by the following two steps: 

(a) Split the population at random into G groups U g (g = 1,..., G) each of size N / G . 

(b) Split each group at random into N / n samples each of size n / G . 

The following table illustrates a rotation pattern with G = 2 , N = 20 , n = 4 and 50% overlap. 

Table 1 : Example of Rotating Scheme 

Group Sample Time 

1 2 3 ... 

U1 s11 √

s12 √ √

... 

s15 

U 2 s21 √ √

s22 √

... 

s25 

Cross-sectional sample (1) 

s = s11 ∪ s21 

(2) 

s = s12 ∪ s21 

(3) 

s = s12 ∪ s22 

Each cross-sectional sample is composed of 2 samples. The sample s11 appears in the first cross-sectional sample 

(time t = 1) , the samples s12 and s21 appear in two time periods(time 2 and 3), while the sample s22 appears only in the 

third cross-sectional sample (time t = 3 ). 

The first step consists on selecting G non-overlap samples from the population with a complete coverage of the union of 

the G samples. Let ag1i (U g ) be the group U g , g = 1,..., G , membership indicator variables with ∑i ag1i (U g ) = Ng , 

∑ g ∑i ag1i (U g ) = N , and ag1i (U g )ah1i (U h ) = 1(g = h) . The second step is a repetition of the first step within each 

group. Let ags2i (s gs ) be the conditional sample s gs membership indicator variables with ∑i ags2i (s gs ) = ngs , 

(t) ˆ(t) ˆ (t) (t) (t) ∑s ∑i ags2i (s gs ) = Ng , and ags2i (s gs )agq2i (s gq ) = 1(s = q) . We have G = G , Y = ∑ g Y ( y , d g (s g )) / G , and 

(T ) T ) (T )

d(s) = vech(d (1) (s 

(1) ),..., d (1) (s (1) ),..., d (T )

(s ),..., d ((sG )) , (3.8) 1 1 G G 1 1 G 

(t) (1) with d (t)

(s

(t)

) = (d (t)(s 

(t) ),..., d (t)

(s 

(t) ))T 

, d (t)

(s

(t)

) = d (U )d (t)(s ) , d ) / E(ag1i (U g )) , g g g1 g gN g gi g g1i g g 2i g g1i (U g ) = J i ag1i (U g 

5 

Page 142 
 

 

 

 

d (t)(s 

(t)

) = a

(t)

(s

(t)

(s 

(t) (t) 

g 2i g gi g 

(t)

) / E(agi g 

(t) )) , and agi (s g ) = ags2i (s gs ) if s gs is the sample used at time t in group g . 

Under rotating samples, the variance of Y

ˆ(t)

 is given by Var(Y

ˆ(t)

) = EVar(Y

ˆ(t)

) +VarE(Y

ˆ(t)

) . Now 

VarE(Y

ˆ(t)

) = Var(Y (t)) = 0 , so an estimator of Var(θ

ˆ

) is given by 

θ L (θ

ˆ

) = ∑ g θ2

(z g , d g 2(s g )) , (3.9) 

where 

(1) ( ) ( ) T (T )

d g 2(s g ) = vech(d g 2(s g 

1 ),..., d g 2

(s g )) , (3.10) 

and θ2 (.) denotes a conditional variance estimator of a total. 

Composite estimation 

Under rotating samples, it is advantageous to use a composite estimator. Suppose the following simple composite 

estimator 

~(t) 

= Y

ˆ(t) 

−αˆ

(t)

(Y

ˆ(t) 

−Y

ˆ(t)

) − β

ˆ (t)

(Y

ˆ(t −1) 

−Y

ˆ(t−1) Y ) . (3.11) u m m 

where Y

ˆ(t) is an estimator of the population total using the unmatched sample, Yˆ(t) and Yˆ(t−1) are estimators of the u m m

β

ˆ (t)

population totals using matched samples, and αˆ

(t)

 and are assumed to be for simplicity regression parameters 

estimators as in the GREG approach. The case of optimal regression is covered in Demnati and Rao (2001). A variance

~ ~ 

estimator of the composite estimator, θ = g(Y ) , can be obtained using (3.9) with 

(T )

( (T )d g 2

(s g ) = vech(d g 

(0) 

2 (s g 

(0) ),..., d g 2 s g )) . (3.12) 

The composite estimator (3.11) has three components as in the AK composite estimator, but it is not recursive and 

estimates for different domains add up unlike in the case of AK composite estimator. 

Concluding Remarks 

We have presented a new approach to variance estimation under longitudinal survey data. A valid variance 

estimator is given under a variety of longitudinal sampling designs. Extension to estimators defined as solutions to survey 

weighted estimating equations is currently under investigation. 

References 

Binder, D. (1996), ΑLinearization Methods for Single Phase and Two-Phase Samples: A Cookbook Approach≅, Survey 

Methodology, 22, 17-22. 

Demnati, A. and Rao, J. N. K. (2001),Α Linearization Variance Estimators for Survey Data≅, Methodology Branch 

Working Paper, SSMD-2001-010E. Statistics Canada. 

Demnati, A. and Rao, J. N. K. (2002), ΑLinearization Variance Estimators for Survey Data With Missing Responses≅, in 

Proceeding of the Section Survey Research Methods, American Statistical Association. 

Kalton, G. and Anderson, D. (1986), “Sampling Rare Populations”, Journal of the Royal Statistical Society, Series A, 149, 

Part 1, pp.65-82. 

Royall, R. M., and Cumberland, W. G. (1981), ΑAn Empirical Study of the Ratio Estimator and Estimators of its 

Variance≅, Journal of the American Statistical Association, 76, 66-77. 

Särndal, C.-E., Swensson, B., and Wretman, J.H. (1989),ΑThe Weighted Residual Technique for Estimating the Variance 

of the General Regression Estimator of the Finite Population Total≅, Biometrika, 76, 527-537. 

Sirken, M. G. (1970),”Household Surveys with Multiplicity”, Journal of the American Statistical Association, 65, 257-266. 

Valliant, R. (1993),ΑPostsratification and Conditional Variance Estimation≅, Journal of the American Statistical 

Association, 88, 89-96. 

6 

Page 143