Multivariate Statistical Modeling with Survey

Data

Tihomir Asparouhov

Muthen & Muthen

Bengt Muthen

UCLA

1This research was supported by SBIR grant R43 AA014564-01 from NIAAA to Muthen

& Muthen.

1
Abstract

We describe an extension of the pseudo maximum likelihood (PML) estimation method developed by Skinner (1989) to multistage strati ̄ed cluster sampling

designs, including  ̄nite population and unequal probability sampling. We conduct simulation studies to evaluate the performance of the proposed estimator.

The estimator is also compared to the general estimating equation (GEE) method

for linear regression implemented in SUDAAN. We investigate the distribution of

the likelihood ratio test (LRT) statistic based on the pseudo log-likelihood value

and describe an adjustment that gives correct chi-square distribution. The performance of the adjusted LRT is evaluated with a simulation study based on the

Behrens-Fisher problem in a strati ̄ed cluster sampling design.

2
1 Introduction

Estimation of simple statistical models such as linear and logistic regressions with survey data is well established and widely used. These models

are however inadequate for analyzing large multivariate data sets that are

being made available by governmental agencies and other research institutions. Increasingly analysts are turning to advanced multivariate models

to better penetrate these complex data structures. Simultaneous regression

equations, structural equation models, time series models, log-linear models,

mixture models, mixed models, latent class models, latent variable models

and combinations of these are frequently the analysts' choice. Methods for

estimating such models however with data obtained by multistage survey

designs are not well established. Frequently analysts use methods designed

for simple random sampling followed by an ad-hoc adjustment for variance

estimation, see Stapleton (2005) for a review of such methods. These methods however are somewhat arbitrary and their theoretical properties are not

well known. Until recently many statistical packages implemented such ad

hoc methods as well, see Asparouhov (2005a).

Skinner (1989) introduced the pseudo maximum likelihood (PML) method

which can be used to estimate any general multivariate parametric model

with data from a complex survey design which includes strati ̄cation, cluster

sampling, and unequal probability sampling with replacement. This method

is in fact applicable to a more general sampling design which includes strati-

 ̄ed multistage sampling with unequal probability sampling with replacement

at the primary sampling stage while allowing for with and without replacement unequal probability sampling on subsequent stages. This sampling

design is known as the WR sampling design and is pioneered and implemented in the software package SUDAAN (RTI, 2002). SUDAAN however

is based on the general estimating equations (GEE) methodology and is only

capable of estimating simple univariate models such as linear and logistic

regression. Mplus, Version 3 (Muthen & Muthen, 1998-2004), implements

the PML method for the WR design and many multivariate models with

normal, discrete and other parametric distributions for observed and latent

variables. More information on Mplus modeling capabilities can be obtained

at www.statmodel.com. Other multivariate modeling packages, such as LISREL (SSI, 2005) have recently adopted the PML method as well.

The three fundamental sampling designs, WR, WOR and WORUNEQ,

pioneered in SUDAAN, are widely used in practice and are being adopted

3
in other software packages. The WOR design is a strati ̄ed multistage sampling design with equal probabilities without replacement sampling at the

PSU level and equal probabilities with or without replacement sampling at

the subsequent stages. The WORUNEQ is a strati ̄ed multistage design

with unequal probabilities without replacement sampling at the PSU level

and with or without replacement equal probabilities sampling at subsequent

stages.

The contributions of this paper are as follows. In Section 2, we expand

Skinner's (1989) PML method to the WOR and WORUNEQ designs. We

also discuss the asymptotic properties of the PML estimator. It is surprising

that this °exible estimation method has not been developed yet for these

common sampling designs. In Section 3 we conduct a simulation study on a

factor analysis model estimated from a two-stage WOR design. In Section

4 we evaluate the performance of four di®erent estimators for samples with

small number of PSUs. The four estimators are PML, implemented in Mplus,

the GEE method implemented in SUDAAN, the GEE with exchangeable

correlation method implemented in SUDAAN and the bias corrected PML

method we propose in this article. In Section 5 we investigate the distribution of the likelihood ratio test statistic based on the pseudo log-likelihood

value and describe an adjustment that gives a correct chi-square distribution.

The e®ects of various complex sampling features on the distribution of the

LRT statistic are illustrated with a simulation study based on the BehrensFisher problem in a strati ̄ed cluster sampling design. All computations are

performed with Mplus, Version 3 (Muthen & Muthen, 1998-2004), unless

explicitly noted.

2 Pseudo Maximum Likelihood Estimation in

Multistage Sampling

In this section we describe the pseudo maximum likelihood estimation for a

general parametric model and the three sampling designs WR, WOR, and

WORUNEQ. Suppose that the log-likelihood for individual i is Li and the

model parameters are μ. Let Ti be the vector of  ̄rst derivative of Li with

respect to μ. Suppose that wi are the weights produced by the complex

sample design, i.e., wi = 1=pi, where piis the probability that individual i

is included in the sample. Let n be the size of the sample population and

4
N be the size of the whole target population. The true model parameters μ0

are de ̄ned as the unique values that maximize the likelihood of the target

population

X

N

L0 = Li:

i=1

The PML ^ estimates μ are de ̄ned as the parameters that maximize the

weighted sample log-likelihood

Xn

L = wiLi:

i=1

These estimates are obtained by solving the weighted score equations

Xn

wiTi = 0:

i=1

For large sample size the weighted sample score equation is an approximation

to the total score equation

X X

n N

wiTi 1⁄4 Ti = 0:

i=1 i=1

(1)

which is solved by ^ the true parameter μ0. Thus the PML estimate μ is a consistent ^ estimate of μ. The asymptotic variance of μ is given by the asymptotic

theory for maximization estimators (see Amemiya (1985), Chapter 4)

(L

00)¡1V ar(L0

)(L

00)¡1

; (2)

where 0 and 00 denote the  ̄rst and the second derivatives of the weighted samP ple log-likelihood. The middle term V ar(L

0

) = V ar(

n

i=1 wiTi) is computed

according to the formulas for the variance of the weighted estimate of the

total described in Cochran, Chapter 11 (1977) or RTI (2002) taking the appropriate design into account. To describe V ar(L

0

) we index the individual

observations by membership in each of the sampling stages. That is, individual i1; i2; i3::: is individual in strata i1, PSU i2, secondary sampling unit i3,

etc. Let ni1:::ilbe the number of sampling subunits in sampling unit i1:::il,

i.e., ni1is the number of PSUs in strata i1, ni1i2is the number of secondary

5
sampling units in PSU i2 in strata i1, etc. Let Zi1:::ir = wi1:::ir Ti1:::ir and let

r be the total number of sampling stages.

X X Zi1:::il = Zi1:::il

il+1 = Zi1:::ir

;

il+1 il+1;:::;ir

1

Z1

i1:::il = Zi1:::il

;

ni1:::il

s

X 1 1 i1:::il = (Zi1:::il+1 ¡ Zi1:::il

)

T

(Zi1:::il+1 ¡ Zi1:::il

):

il+1

Suppose that

1⁄2

f

¤ =

fi1:::il

i1:::il

if the sampling in the i1:::il unit is WOR

0 otherwise

For the WR design, regardless of the number of sampling stages, the variance

of the score is given by

X n

V ar(L

0

) = i1s :

n ¡ 1

i1

i

i1 1

For the WOR design, for compactness, we describe the variance of the score

for a strati ̄ed 3 stage sampling design

V ar(L

0

) = V1 + V2 + V3;

where

X ¤ ni1i2

V2 = (1 ¡ f )fi

¤

1 i1i2

si1i2

n ¡ 1 i1;i2

i1i2

X n

V1 = (1 ¡ f

¤

)

i1

i1

si1

ni1 ¡ 1 i1

¤ ¤ 1 2 3 V3 = (1 ¡ fi1i2i3

)fi1fi

¤

1i2

si1i2i3

n ¡ 1 i1;i2;i3

i1i2i3

X ni i i

For the WORUNEQ design we describe the variance of the score again for a

strati ̄ed 3 stage sampling design. The probability that PSU i2 in stratum

i1 is selected is denoted by pi2 i1. The probability that both PSUs i2 and i

0

j 2

in stratum i1 are selected in the sample is denoted by pi2i

0

2

ji1

V ar(L

0

) = V1 + V2 + V3;

6
where

X

XX X V1 =

pi2ji1pi

0

2

ji1 ¡ pi2i

0

2

ji1

(Zi1i2 ¡ Zi1i

0

2

)

2

pi2i

0

2 i1 i2 i

0 ji1

2>i2

X ni1i2 V2 = (1 ¡ f

¤

)pi2ji1 i1i2si1i2

ni1i2 ¡ 1 i1;i2

ni1i2i3

f

¤ V3 = (1 ¡ f¤

)pi2ji1si1i2i3:

i1i2i3 i1i2 ni1i2i3 ¡ 1 i1;i2;i3

The above estimation method hinges on the approximation (1) of the

total score, which can be achieved if the number of PSU units is large and the

residuals of the score estimation within each PSU units satisfy Lindeberg's

extension of the central limit theorem (see Feller, 1968). If the number of PSU

units is small however the PML parameter estimates can be substantially

biased.

3 Factor Analysis Simulation Study

In this section we will evaluate the performance of the PML estimator for a

two-stage WOR design for a factor analysis model. The model is described

as follows

Yij = 1j +  ̧j ́i + "ij (3)

where i = 1; :::; n, n is the sample size, and j = 1; :::; 5, i.e., the observed vector for each individual is of dimension 5. Here 1j

is the intercept parameter,

 ̧j is the loading parameter,  ́iis the factor variable, and "ij is the residual

variable. The variables  ́i and "ij are normally distributed zero-mean variables with variances Ã and μj respectfully. The parameters we use for this

simulation study are as follows

£=(11; :::; 15;  ̧1; :::;  ̧5; μ1; :::; μ5; Ã) = (4)

(2; 2:7; 3:3; 4:5; 5:5; 1; 0:7; 1:3; 1:5; 0:5; 1; 1; 1; 1; 1; 1:2):

First we describe the target and the sample populations. We generate a multivariate target population of 50000 individuals with 5 normally distributed

outcomes with mean and variance given by model (3) with parameter values

given by (4). We impose the following two-level population structure on the

target population. We group the observations into 140 PSUs, the  ̄rst 120

7
Table 1: Bias of PML Parameter Estimates for Factor Analysis Model

n m 13  ̧3 μ3 Ã

200 20 0.054 0.045 -0.028 -0.129

500 50 0.007 0.009 -0.011 -0.035

1000 100 0.001 -0.003 -0.004 0.003

1400 140 0.003 -0.001 -0.001 0.003

are of size 250 and the remaining 20 are of size 1000. The observations are

not placed at random in the PSUs. They are placed according to an ordering

based on a function f. That is, the  ̄rst 250 observations with the highest values of f are placed in PSU 1, the second highest 250 are placed in PSU 2, etc.

After all 120 PSUs of size 250 are formed we form the remaining 20 PSUs of

size 1000 again according to the order given by the f function. This method

of constructing target population was used in Smith and Holmes (1989). The

choice of f is to some extent critical to the type of sampling we get. Suppose

that f is instead a random function independent of Y . The multi-stage sampling will then be equivalent to simple random sampling (SRS). In a model

with dependent variable Y and independent variable X, a function f that

depends only on X but not on Y produces non-informative random sampling.

The only way to produce informative sampling is to choose f which depends

on Y in addition to perhaps other variables. In this target population we

P choose fi = j Yij , which clearly induces informative sampling.

The target population is sampled with a two-stage WOR design. Equal

probability sampling is used at each stage. We vary the number m of PSUs

included in the sample while the number of units sampled from the i¡th

PSU remains constant, ni = 10. The total sample size is thus n = 10m.

The ratio between the sampling weights in the large PSUs and the sampling

weights in the small PSUs is 4. We use 500 replications, i.e., we sample the

target population 500 times and calculate the PML estimates and their 95%

con ̄dence intervals.

Table 1 shows the bias of the PML parameter estimates and Table 2 shows

the coverage of the PML con ̄dence intervals. We see that the performance of

the PML method is very good, bias is almost non-existent and the coverage

for the con ̄dence intervals is in line with expectation. The only exception

is the estimation of the Ã parameter which has larger bias and consequently

8
Table 2: Coverage of PML 95% Con ̄dence Intervals for Factor Analysis

Model

n m 13  ̧3 μ3 Ã

200 20 0.882 0.908 0.912 0.746

500 50 0.940 0.924 0.928 0.850

1000 100 0.950 0.950 0.946 0.926

1400 140 0.954 0.948 0.968 0.952

lower con ̄dence interval coverage for samples with small number of PSUs.

In the next section we explore further the bias that arises in model estimation

from samples with small number of PSU.

4 Small Number of PSUs

In this section we explore di®erent estimation techniques for dealing with

bias that arises in small number of PSUs samples. We conduct a simulation

study similar to the simulation study conducted in the previous section. For

simplicity we use a two-stage WR design on a smaller target population. Here

again equal probability sampling is used at each stage. The target population

of size 10000 is generated as in the previous section and 14 PSUs are formed,

12 of size 500 and 2 of size 2000. The sample population again has a varying

number m of PSUs while the number of units sampled from the i¡th PSU

remains constant, ni = 50. The total sample size is n = 50m. We use 500

replications in this simulation study as well. To be able to compare various

estimating techniques we choose a basic regression model for Y1 and Y2

Y1 = ® +  ̄Y2 + ":

Using the whole target population we get the true values of the parameters

as ® = 0:56 and  ̄ = 0:54. The variance parameter of the residual " is not

included in this investigation because the methods implemented in SUDAAN

do not provide an estimate for this parameter. We compare four di®erent

estimation methods. The  ̄rst method is the PML method implemented in

Mplus. The second method is the GEE method implemented in SUDAAN.

This method is based on general estimating equations which are identical to

9
the PML score equations and as our simulation study con ̄rms the results

produced by the two methods are identical. This observation is valid also

for other models such as logistic regression. The third method is the GEE

method with exchangeable correlation implemented in SUDAAN. We denote

this method by GEE-Ex. The RSTEPS parameter that this method depends

on did not a®ect the results in our simulation signi ̄cantly and thus we only

report the results we obtained with the default RSTEPS=1.

The forth estimation method we examine in this simulation study is a

bias corrected PML method (BC) that we describe here. The  ̄rst step of

the BC estimation method is to construct estimates for the mean and the

variance/covariance of Y1 and Y2, by estimating the bias of the PML mean

and variance/covariance estimates. We illustrate this for a single Y variable.

The PML estimate for the mean is

P

i wiYi

1^PML = P :

i wi

The BC estimate for the mean is then

i wiYi

1 ^ ^BC = P ¡ C0;

i wi

P

where C^

0 is an estimate for the bias C0 of the PML estimate, i.e., if 1 is the

mean of Y

μP ¶

i wiYi C0 = E

P ¡ 1:

i wi

The term C0 is of the form

Z^

1 E(Z^1)

C = ¡ :

Z^

2 E(Z^2)

Formula 6.33 in [C] provides a method for estimating such a quantity. An

asymptotic estimate for C is

Ã ! 1 Z^

1

V ar(Z^

2) ¡ Cov(Z^1; Z^2) :

Z^2 ^

2 Z2

Both Z1 and Z2 are estimates of the total quantity for the variables Y and

the constant variable 1. Thus the variance/covariance terms above can be

estimated just as the variance/covariance of the total score estimates given

10
in Section 2. These estimates take into account the sampling design. The

PML estimate for the variance is

P ÃP !2

i wiY

2

i i wiYi

v^PML = P ¡ P :

i wi i wi

The BC variance estimate is

P ÃP !2

i wiY

2

i wiYi ^ v^BC = P

i ¡ P ¡ C1;

i wi i wi

where C^

1 is an estimate of the second order moment bias of the  ̄rst term

P 2 P

i wiYi = i wi constructed just as the bias estimate for the mean. The covariance term is estimated from the multivariate version of the above formula.

Once the mean and the variance/covariance estimates for Y1 and Y2 are constructed, we estimate the parameters μ=(®, ̄) by minimizing the quasi ML

 ̄t function

F(μ) = tr(^vBC v(μ)

¡1

)¡log jv^BCv(μ)

¡1

j+(^1BC ¡1(μ))Tv(μ)

¡1

(^1BC ¡1(μ));

where 1(μ) and v(μ) are the vector mean and variance of the (Y1; Y2) vector

expressed in terms of the model parameters μ and the following auxiliary

parameters: the mean Y2, the variance of Y2, and the variance of the residual

in the above regression equation.

We study the properties of these four estimators for samples with small

number of PSUs. Tables 3 and 4 show the bias and the MSE of the four

estimators on samples with 5, 10, 15, and 20 PSUs. The PML method and

the GEE method, as expected, produce identical results not only on average

but in individual replications as well and are reported in the same column.

The bias of the PML/GEE estimator is present for both the intercept and the

slope even for m = 20 but as expected this bias decreases as the number of

PSUs increases. The bias of the BC estimator is almost non-existent except

for m = 5. The BC method outperforms the PML/GEE estimator in terms

of both MSE and bias in this simulation. The BC method, however, may

not outperform the PML method in all situations. Examples in Cochran

(1977) show that sometimes this method reduces the bias while increasing

the MSE of the estimates. The estimator GEE-Ex performs very poorly.

This method produces large bias for both parameters and large MSE. It

seems also that this bias does not disappear as the number of sampled PSUs

11
Table 3: Bias and Mean Squared Error for the Intercept in Linear Regression.

m

PML/GEE

Bias

PML/GEE

MSE

GEE-Ex

Bias

GEE-Ex

MSE

BC

Bias

BC

MSE

250 5 0.434 0.632 1.576 2.855 0.179 0.442

500 10 0.220 0.281 1.634 2.889 0.016 0.211

750 15 0.132 0.185 1.664 2.902 -0.028 0.157

1000 20 0.103 0.131 1.675 2.917 -0.026 0.111

Table 4: Bias and Mean Squared Error for the Slope in Linear Regression.

n m

PML/GEE

Bias

PML/GEE

MSE

GEE-Ex

Bias

GEE-Ex

MSE

BC

Bias

BC

MSE

250 5 -0.237 0.130 -0.672 0.466 -0.133 0.123

500 10 -0.125 0.063 -0.656 0.440 -0.036 0.059

750 15 -0.072 0.040 -0.642 0.419 -0.002 0.038

1000 20 -0.061 0.027 -0.649 0.427 -0.004 0.024

increases. A simulation study based on a logistic regression model produced

the same results. The PML/GEE method performed well as the number

of PSUs increases for logistic regression as well. In contrast the GEE-Ex

method produced large bias regardless of the number of PSUs in the sample.

5 Likelihood Ratio Test in Multistage Sampling

Hypotheses involving several parameters are frequently tested in multivariate

modeling. Wald's test can be used for such testing if the asymptotic variance/covariance of the parameter estimates is available. Wald's test however

requires additional calculations, which sometimes are quite complex. One

such example is the test of a factor analysis model against an unrestricted

covariance model. When maximum-likelihood estimation is performed however the likelihood ratio test (LRT) can be obtained without additional com12
putations and this test is frequently used for complex hypothesis testing.

In this section we show how the pseudo maximum likelihood can be used

to perform LRT for multistage sampling designs. The distribution of the

LRT statistic based on the maximized weighted log-likelihood value is not

a chi-square distribution. This distribution depends on the sampling design

just as the asymptotic covariance of the parameter estimates depends on the

sampling design. Here we describe an adjustment of the LRT statistic which

takes into account the sampling design and produces a test statistic with a

chi-square distribution. This adjustment is constructed similarly to the adjustments of the Yuan-Bentler (2000) and the Satorra-Bentler (1988) robust

chi-square tests for mean and variance structures. Similar  ̄rst and second

order adjustments are described also in Rao-Thomas (1989) for contingency

tables.

We assume a general hypothesis testing for two nested models M1 and

M ^

2. Let μi be the true parameter values and μi the parameters estimates for

model Mi that maximize the pseudo log-likelihood function Li. Let di be the

number of parameters in model Mi. The corrected LRT statistic is

T

¤ = c ¢ 2(L1 ¡ L2); (5)

where c is the correction factor

T r((L

00

1

)

¡1V ar(L1

0

)) ¡ T r((L

00

2

)

¡1V ar(L

0

2

)) c =

d1 ¡ d2

: (6)

The statistic T

¤ has approximately a chi-square distribution with d1 ¡ d2

degrees of freedom. The components T r((L

00)¡1

i V ar(L

0

i

)) are easily available

since they are part of the asymptotic covariance for the parameter estimates

given in (2). Justi ̄cation for this approximation is given in the Appendix.

We demonstrate the importance of the LRT adjustment with a simple

simulation study which incorporates both cluster and strati ̄ed sampling. For

simplicity we use a single outcome variable and compare the mean and the

variance of this outcome across two groups. Each of the two groups contains

three strata. Within each stratum we sample at random entire clusters. For

example the two groups can be private and public schools, the strata can

be di®erent regions in the country, the clusters can be the classrooms and

the students can be the individual observations. While in this example the

groups actually contain entire strata and clusters, this doesn't necessarily

have to be the case. For example the grouping variable could be gender

which is not nested above the strata and the cluster variables.

13
All six strata in our simulation study have equal size and we sample 200

observations from each by cluster sampling. Within each stratum the clusters

are of equal size. We denote the size of the clusters in stratum s in group g

by nsg. The cluster sizes in the six strata are as follows n11 = 5, n21 = 10,

n31 = 20, n12 = 10, n22 = 20, n32 = 40. The distribution of observations i in

cluster j in stratum s in group g is described by

Yijsg = 1sg +  ́jsg + "ijsg

where  ́jsg and "ijsg are zero mean normally distributed variables with variance 1, and the parameters 1sg are as follows 111 = 1, 121 = 2, 131 = 3,

112 = 0, 122 = 2, 132 = 4. Given our choice of parameters the total mean

in the two groups is 2. The total variance of y is, however, larger in the

second group. We test two hypotheses by LRT. The  ̄rst hypothesis T1 is

that the means in the two groups are equal and is also known as the BehrensFisher problem, see Sche®e (1970). The second hypothesis T2 is that both the

means and the variance parameters are equal in the two groups. The  ̄rst

test should not reject the hypothesis because the means are indeed equal.

The second test should, however, reject the hypothesis because the variances

are not equal. In addition the test statistic T1 should have a chi-square distribution with 1 degree of freedom because it tests just one constraint. Test

statistic T2 has two degrees of freedom because it tests two constraints. The

null hypothesis for the second test is not correct however and therefore the T2

test statistic is not expected to have a chi-square distribution with 2 degrees

of freedom. This test statistic is expected to be su±ciently large so that the

test is rejected.

To evaluate the e®ect of strati ̄cation and clustering on the test we compare  ̄ve di®erent methods for computing the LRT statistic. These methods

are as follows.

2 Method A. Adjusted robust LRT which takes both the clustering and

the strati ̄cation into account.

2 Method B. Adjusted robust LRT which takes only the clustering into

account and ignores the strati ̄cation.

2 Method C. Adjusted robust LRT which takes only the strati ̄cation into

account and ignores the clustering.

2 Method D. Adjusted robust LRT which ignores both the clustering and

the strati ̄cation.

14
2 Method E. Unadjusted LRT.

The results of the simulation study are presented in Table 5. We report

the average values of the T1 and T2 test statistics over 500 replications and the

rejection rates for the two tests based on the 5% rejection level. As expected

method A performs correctly producing a test statistic T1 with an average

value of approximately 1 and rejection rate of approximately 5%, while all

the other methods produced erroneous results. From the table we clearly

see that including the strati ̄cation information results in an increase of the

LRT statistic and the rejection rates, while including the cluster information

decreases the LRT statistic and the rejection rates. The result of not including the strati ̄cation information in the  ̄rst test is that there are virtually

no rejections, while the result of not including the cluster information is that

the test rejects the null hypothesis incorrectly an additional 47% of the time

above the nominal 5% level. Methods D and E both produce rejection rates

that are too high and in our simulation the results of the two methods are

quite close.

The most important e®ect of strati ̄cation is actually seen in the second

test. Methods C, D and E all have in°ated power largely because the clustering information is ignored. Method A rejects 76% of the time for this sample

size. As the sample size increases this rejection rate converges to 100%. Not

including the strati ̄cation information in method B results in a decrease of

power. As a result of that, method B does not reject the second hypothesis

as it should an additional 26% of the time.

It is clear from Table 5 that the sampling features in complex sampling

designs can a®ect dramatically the distribution of the LRT statistics and

erroneous conclusions can be reached if the sampling features are not accounted for. The adjusted LRT statistic provides an e®ective solution for

hypothesis testing with complex sampling data. The LRT adjusted statistic

is implemented in Mplus, Version 3 (Muthen & Muthen, 1998-2004) for a

wide variety of models and complex sampling designs.

6 Conclusion

In this article we demonstrated how the PML estimator can be used with

the three basic complex sampling designs WR, WOR and WORUNEQ. The

PML estimator can be utilized in advanced multivariate statistical modeling

to properly account for various features of complex sampling designs. The

15
Table 5: E®ect of Strati ̄cation and Clustering on the Chi-Square Test

Method A B C D E

T1 Average 1.042 0.349 9.141 5.052 4.984

T1 Rejection 0.054 0.002 0.524 0.380 0.380

T2 Average 12.827 8.057 75.884 61.236 53.856

T2 Rejection 0.760 0.500 0.990 0.982 0.980

PML parameter estimates are a®ected only by the sampling weights while

their standard errors are adjusted to re°ect the e®ects of strati ̄cation, cluster sampling, multistage sampling,  ̄nite population sampling and unequal

probability sampling. Our simulation studies showed that the PML method

performs very well as long as the number of PSUs is not small. When the

number of PSUs is small alternative estimator such as the bias corrected PML

method described here are preferable. Our comparison with the method implemented in SUDAAN showed that the GEE method is equivalent to the

PML method for linear and logistic regression. The GEE with exchangeable

correlation method performed poorly in our simulation study. The main advantage of the PML method however is its generality. This method can be

used to estimate any parametric model.

In this article we described an adjustment to the LRT statistic which

takes into account the complex sampling design. The unadjusted LRT can

lead to erroneous results when analyzing survey data, while the adjusted

LRT performs correctly. Because of its simplicity of use, the adjusted LRT

is a valuable alternative to other methods such as Wald's test.

The PML extension described in this article and the LRT adjustment

can also be used for multilevel models via the multilevel pseudo maximumlikelihood method described in Asparouhov (2005b).

16
7 Appendix

We follow the ideas of Yuan-Bentler (2000) to derive a general LRT correction

based on the PML method under complex sampling. The arguments below

also apply to any consistent estimator obtained by maximizing an objective

function l. Such estimators are called extremum estimators; see Amemiya

(1985), Chapter 4.

We assume a general hypothesis testing for two nested models M1 and

M ^

2. Let μi be the true parameter values and μi the parameters estimates

for model Mi that maximize the pseudo log-likelihood function Li. We are

interested in the asymptotic distribution of the ^ test statistic T = 2(L2(μ2)

^

¡

L1(μ1)) that can be used to test the more restricted model M1 versus the less

restricted model M2. More speci ̄cally we are interested in the asymptotic

distribution of T when M1 is correct. Since M1 is correct μ2 is a function of

μ1 and L1(μ1) = L2(μ2). Let ¢ = @μ2=@μ1. Let Si = @Li(μ)=@μi and Hi =

¡n

¡1@2Li(μ)=(@μi)2

. Given some basic regularity conditions (see Amemiya,

Theorem 4.1.3) we have that

p

n(μ

^

i ¡ μi) = Op(1); (7)

where n is the number of observations. Using the Taylor expansion we get

that

1

Li(μ

^

i) = Li(μi) + Si(μi)(μ

^

i ¡ μi) ¡ n(μ

^

i ¡ μi)

THi(μi)(μ

^

i ¡ μi) + op(1)

2

(8)

and

p

0 = Si(μ

^

i) = Si(μi) ¡ nHi(μi)(μ

^

i ¡ μi) + op( n) (9)

Solving equation (9) for Si(μi) and substituting that in (8) gives us

1

Li(μ

^

i) = Li(μi) + n(μ

^

i ¡ μi)

THi(μi)(μ

^

i ¡ μi) + op(1)

2

(10)

Now

T = n(μ

^

2 ¡ μ2)

T H2(μi)(μ

^

2 ¡ μ2) ¡ n(μ

^

1 ¡ μ1)

TH1(μi)(μ

^

1 ¡ μ1) + op(1) (11)

The chain rule for di®erentiation gives us

S1 = ¢S2: (12)

17
Solving (9) for Si(μi) and substituting in (12) we get that

p p H1(μ1) n(μ

^

1 ¡ μ1)=¢H2(μ2) n(μ

^

2 ¡ μ2) + op(1) (13)

^ Solving now equation (13) for p

n(μ1 ¡ μ1) and substituting in (11) we get

μ ¶

T = n(μ

^

2 ¡μ2)

T H2(μ2)¡ H2(μ2)¢T H¡1

(μ1)¢H2(μ2) (μ

^

2 ¡μ2) +op(1) 1

(14)

From equation (9) we also see that the asymptotic distribution of

p

n(μ

^

i ¡ μi) ! N(0; Vi) (15)

where

Vi =

1

Hi(μi)

¡1V ar(Si(μi))Hi(μi)¡1

n

(16)

Elementary matrix algebra shows that the asymptotic distribution of T is

X ̧iÂ

2

1i

i

(17)

where Â

2

1i

are independent chi-square distributed random variables and  ̧i

are the eigenvalues of

μ ¶

E = V2 H2(μ2) ¡ H2(μ2)¢T H1

¡1

(μ1)¢H2(μ2) : (18)

The p-values of this distribution are easy to compute following a method

^ developed in Imhof (1961). Because μ1 and μ2 are not known we use μ1 and

^μ2 in equation (18).

p

By equation (9) we get that Si(μi) = Op( n): The chain rule for the

second derivative gives us

H1(μ1)=¢TH2(μ2)¢ + n

¡1S2@2

μ2=(@μ1)

2 = ¢T H2(μ2)¢ + op(1) (19)

This leads us to the following alternative computation of E

μ ¶

E2 = V2 H2(μ2) H2(μ2)¢T(¢H2(μ2)¢) ¡

¡1¢TH2(μ2) = E + op(1): (20)

While asymptotically equations (18) and (20) are equivalent, they will lead

to di®erent results for  ̄nite sample size. It is not clear which one of the two

should be preferred in speci ̄c applications.

18
Instead of computing the exact p-value of the weighted chi-square distribution (17) we can use the following adjusted test statistic. Let

d d

T

¤ = P

T = T:

i  ̧i T r(E)

(21)

where d is the number of parameter restrictions model M1 imposes, i.e., d

is the di®erence between the number of parameters in the two models. The

ratio

T r(E)

c =

d

is the correction factor. The distribution of T

¤

is approximated by a chisquare distribution with d degrees of freedom and thus its p-values are readily available. Again we can use E2 in formula (21) instead of E and get an

asymptotically equivalent statistic which in  ̄nite sample size may be substantially di®erent.

Now we derive one more formula for computing T

¤

. Using equations (12)

and (16) we get that

H1(μ1)V1H1(μ1)=¢H2(μ2)V2H2(μ2)¢T

: (22)

Now using formula (18) and (22) we get that

T r(E) = T r(V2H2(μ2)) ¡ T r(V2H2(μ2)¢TH

¡1

1

(μ1)¢H2(μ2)) =

T r(V2H2(μ2)) ¡ T r(¢H2(μ2)V2H2(μ2)¢T H

¡1

1

(μ1)) =

T r(V2H2(μ2)) ¡ T r(V1H1(μ1)):

Again ^ ^ since μ1 and μ2 are not know we approximate with μ1 and μ2

T r(E) = T r(V2H2(μ

^

2)) ¡ T r(V1H1(μ

^

1)): (23)

Formula (23) is the same as formula (5) and is also the formula implemented

in Mplus. This formula has several advantages. It is computationally more

e±cient then formulas (18) and (20) because it does not involve the computation of ¢. It can also be used to easily compute the proper LRT adjustment

when two nested hypothesis are involved as follows. Suppose that we have

three models M1, M2 and M3 and we have the test statistics T1

¤ and T2¤

for

testing M1 versus M3 and M2 versus M3. Suppose that the LRT statistics

19
have been computed according to formulas (21) and (23). Let the correction factors be c1 and c2 and the degrees of freedom d1 and d2. We want

to compute the LRT statistic T

¤

for testing M1 versus M2. Let the degrees

of freedom for that test be d and the correction factor be c. We have that

d = d1 ¡ d2 and

cd = T r(V2H2(μ

^

2)) ¡ T r(V1H1(μ

^

1)) = (T r(V3H3(μ

^

3)) ¡ T r(V1H1(μ

^

1)))¡

(T r(V3H3(μ

^

3)) ¡ T r(V2H2(μ

^

2))) = c1d1 ¡ c2d2:

Thus

c1d1 ¡ c2d2

c =

d

and

c1T

¤ ¡ c2T¤

T

¤ =

1 2

:

c

The exact same approach was outlined in Satorra-Bentler (1999) when applied to the Satorra-Bentler (1988) chi-square statistic.

20
8 References

Amemiya, T. (1985). Advanced Econometrics. Harvard University Press.

Asparouhov, T. (2005a). Sampling Weights in Latent Variable Modeling.

Structural Equation Modeling, 12, 411-434.

Asparouhov, T. (2005b). General Multilevel Modeling with Sampling

Weights. Accepted in Communications in Statistics{Theory and Methods.

Cochran, W. G.(1977). Sampling Techniques. John Wiley & Sons, third

edition.

Feller, W. (1968) Introduction to Probability Theory and its Applications, vol. 1. John Wiley & Sons, third edition.

Imhof, J.P. (1961) Computing the Distribution of Quadratic Forms in

Normal Variables. Biometrika 48, 419-429.

Muthen, L.K. and Muthen, B.O. (1998-2005). Mplus User's Guide. Third

Edition. Los Angeles, CA: Muthen & Muthen

Rao, J. N. K., & Thomas, D. R. (1989). Chi-Square Tests for Contingency Table. In Analysis of Complex Surveys (eds. C.J.Skinner, D.Holt and

T.M.F. Smith) 89-114, Wiley.

Research Triangle Institute (2002). SUDAAN User Manual Release 8.0,

Second Edition.

Satorra, A., & Bentler, P.M. (1988). Scaling Corrections for Chi-Square

Statistics in Covariance Structure Analysis. Proceedings of the Business and

Economic Statistics Section of the American Statistical Association, 308-313.

Satorra, A., & Bentler, P.M. (1999). A Scaled Di®erence Chi-square Test

Statistic for Moment Structure Analysis. UCLA Statistics Series # 260.

http://www.stat.ucla.edu/papers/preprints/260/

21
Sche®e, H. (1970). Practical Solutions of the Behrens-Fisher Problem.

Journal of the American Statistical Association, 65, 1501-1508.

Scienti ̄c Software International (2005). LISREL, Version 8.7.

Skinner, C. J. (1989). Domain Means, Regression and Multivariate Analysis. In Analysis of Complex Surveys (eds. C.J.Skinner, D.Holt and T.M.F.

Smith) 59-87, Wiley.

Smith, T., and Holmes, D. (1989) Multivariate Analysis. In Analysis of

Complex Surveys (eds. C.J.Skinner, D.Holt and T.M.F. Smith) 165-190, Wiley.

Stapleton, L. (2005) An Assessment of Practical Solutions for Structural

Equation Modeling with Complex Sample Data. Accepted in Structural

Equation Modeling.

Yuan, K., & Bentler, P. M. (2000) Three Likelihood-Based Methods for

Mean and Covariance Structure Analysis With Nonnormal Missing Data.

Sociological Methodology 30, 167-202.

22