Last month Phoenix Do posted to the list a question about whether -gllamm-
or another Stata command can handle cross-classified random effects.
Laurent Glance replied that he did not think that -gllamm- (which is set up
primarily for hierarchical models) could handle cross-classified random
effects models. I believe that -gllamm- *can* estimate crossed-random
effects models.
Harvey Goldstein describes how to set up the crossed-random effects model
for a hierarchical software package such as HLM (or MLwiN, or -gllamm-) in
Chapter 8 of his _Multilevel Statistical Models_ Second Edition. (London:
Edward Arnold, 1995), esp. pp. 116-17 and pp. 123-24. (This book is the
third in the series, Kendall's Library of Statistics, and is co-published in
the U.S.A. by Halsted Press, part of John Wiley & Sons.)
Harvey Goldstein's approach is to choose one of the crossed random effects
as the second level, and to take the other crossed effect and create dummy
(indicator) variables of it as random effects at a third level, and
to -constraint- the variance covariance matrix of this second random effect
so that the diagonal elements are equal and off-diagonal elements are zero.
I've illustrated the approach in the do-file below using the dataset from P.
E. Shrout & J. L. Fleiss, Intraclass Correlations: Uses in Assessing Rater
Reliability. _Psychological Bulletin_ 86:420-28, 1979.
Because there are multiple random effects at the third level, you'll need to
use the -eqs- option to specify equations for them. And because you need to
use -eqs-, you'll need to specify the equation for the random effect at the
second level of the hierarchy, as well. To do that, use the method
described in -gllamm-'s user's manual: create a variable equal to one so
that an intercept-only constant can be defined for the equation for the
random effect at the second level. Run -gllamm, noest- to get the names of
the elements of the variance-covariance matrix, and set up -constraints- to
make for the matrix structure of the random effects at the third level of
the hierarchy, which are really a dummy-variable rendition of the second of
the two crossed random effects. -gllamm- takes care of the
variance-covariance matrix at the second level, and you'll see the structure
of the block-diagonal Cholesky decomposition matrix with that random effect
represented as a single element and with zero covariance between it and
elements for the random effect at the third level (the crossed random
effect). The latter are set up with zero covariance, and this is seen in
the Cholesky decomposition matrix, as well. Use -matrix list e(chol)- after
convergence to see all of this.
Now, with the Shrout & Fleiss small dataset of four raters, -gllamm- will
need to estimate five random effects. -gllamm- takes its time estimating
numerous random effects, and if you have dozens of neighborhoods, it will
still be worthwhile to use an alternative software package, such as one of
those that Laurent mentioned. (Or hold off until Stata releases their
offering for mixed models, as Bobby Gutierrez mentioned on Monday.)
In addition, if you're working with normal (linear) models and are
interested in the three variance components in their own right, such as to
estimate a judges-as-random-effect intraclass correlation coefficient for an
unbalanced dataset, then -gllamm-'s maximum likelihood method will give more
biased estimates than will restricted maximum likelihood (REML) despite that
there's only a single fixed effect (the intercept, the two populations'
joint mean). So, for a project that I just finished to estimate ICC(2,1)
for an unbalanced dataset of physicians evaluating patients, I needed to
resort to S-Plus, anyway. (It's a little awkward to fit a simple unblocked
crossed-random effects model in R/S-Plus, too, by the way. See
https://stat.ethz.ch/pipermail/r-help/2003-October/039775.html , but also
http://maths.newcastle.edu.au/~rking/R/help/03b/6412.html .)
Joseph Coveney
----------------------------------------------------------------------------
clear
set more off
input byte target byte judge byte score
1 1 9
1 2 2
1 3 5
1 4 8
2 1 6
2 2 1
2 3 3
2 4 2
3 1 8
3 2 4
3 3 6
3 4 8
4 1 7
4 2 1
4 3 2
4 4 6
5 1 10
5 2 5
5 3 6
5 4 9
6 1 6
6 2 2
6 3 4
6 4 7
end
/* A constant is needed for the equation for the
intercept of the random effect at the second level. */
generate byte targets = 1
eq k: targets
/* Create dummy (indicator) variables for the crossed random effect
that will be assigned to the third level of the hierarchical model,
and assign each to an equation for intercept of a random effect. */
tabulate judge, generate(rater)
eq r1: rater1
eq r2: rater2
eq r3: rater3
eq r4: rater4
/* Define constraints to hold the diagonal elements of the Cholesky
decomposition matrix equal, so that a single variance will be
estimated for this crossed random effect, which is represented
by numerous dummy variables (numerous random effects) at the third
level. */
constraint define 1 rater1 = rater2
constraint define 2 rater2 = rater3
constraint define 3 rater3 = rater4
/* As an aside, you can use the constraints below in lieu of
-gllamm-'s -nocorrel- option; I've chosen the latter
in this illustration. */
constraint define 4 [con2_2_1]_cons = 0
constraint define 5 [con2_3_1]_cons = 0
constraint define 6 [con2_3_2]_cons = 0
constraint define 7 [con2_4_1]_cons = 0
constraint define 8 [con2_4_2]_cons = 0
/* Since the crossed random effect that's assigned to the
third level is represented by a bunch of dummy
variables, you'll need to create a constant here, too,
in order to define (identify) the third level's single
random effect (the crossed random effect that you're
assigning to the third level of the hierarchy).
You could go ahead and use the same constant that you've
created for the second level's equation. But the
printout would be confusing, so create another constant with
a more appropriate name for the random effect at this level. */
generate byte judges = 1
/* Estimate the model; I've chosen three integration points
in the interests of time--it still takes a while. */
gllamm score, i(target judges) nrf(1, 4) eqs(k r1 r2 r3 r4) nocorrel ///
constraints(1 2 3) nip(3) adapt trace allc
/* If your'e willing to wait, uncomment the next two lines.
matrix A = e(b)
gllamm score, i(target judges) nrf(1, 4) eqs(k r1 r2 r3 r4) nocorrel ///
constraints(1 2 3) adapt from(A) skip trace allc */
matrix list e(chol)
exit
----------------------------------------------------------------------------
From: "Glance, Laurent" <[email protected]>
Subject: st: RE: cross-classified random effects
Date: Tue, 23 Mar 2004 07:05:19 -0500
I assume that your outcome variable is binary.
Gllaam will not handle cross-classified data structures. MLWIN can handle
cross-classified models, but the syntax for cross-classified models is not
well documented. I have found PROC GLIMMIX to work well with
cross-classified data structures. The best reference for PROC GLIMMIX is
SAS System for Mixed Models by Ramon Littell - available through the SAS web
site.
Larry
-----Original Message-----
From: Do, Phoenix [mailto:[email protected]]
Sent: Monday, March 22, 2004 6:42 PM
To: [email protected]
Subject: st: cross-classified random effects
Hello,
I am trying to model risk behaviors for IDUs. Since my data is hierarchical
in nature, I want to apply a multilevel model using gllamm.
My dataset consists of an unbalanced panel in which some
people are observed only once while others are observed 2 to
10 times. We want to account for "neighborhood"
characteristics so we have neighborhood level variables in
our model. So level 1 would be the individual level time
series data. Level 2 would be person level. Level 3 would be
the neighborhood.
However, these people move and are not necessarily in the
same neighborhood throughout our study.
I believe what I need to do is use a crossed random effects
model.
Can you do this in gllaam? And if so, how? I haven't been able to find any
reference to this in the manual.
Thank you for any help you can offer,
Phoenix
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
