I wrote: "Why wouldn't it just be -anova lntime case modality reader-?
. . ."
As David points out in his follow-up post about the differences in model
specification between StataCorp and UCLA, there can be an interest a
fixed-by-random interaction--a random-slope--in addition to the random
intercept. A classic example would be a clinic-by-treatment interaction in
a model that considers clinic as a random effect. (Is there variation in
treatment effect between clinics?) In Peter's case, it would be, "Is there
variation between physicians (patients) in magnitude of display-method
differences?" The interaction of the two random effects, physicians and
patients, is also estimable with ANOVA given that he has replications. But
this variance component is liable to be small, so Peter might have
convergence difficulties with -xtmixed-. Even those for display
method-by-patient and display method-by-physician might be small, or even
negative by chance, which will give rise to unpleasantries with -xtmixed-.
Because he has a balanced dataset, it will be better for him to set the
observed mean squares from -anova- equal to the expected mean squares and
solve for the sigma squares (sigma squareds?) to get these variance
components.
Note that for -xtmixed- (for -gllamm-, too), generating the interaction
indicator variable for a random slope that involves a categorical variable
in the fixed effect does not use the -xi- command. Instead, for nominal
categorical fixed effects that have more than two levels, use -egen
fixedXrandom, group(fixed random)-. (This egen-group interaction term
shouldn't then be used with -anova-; instead, use -anova-'s normal syntax
for specifying interaction terms.)
The do-file below illustrates this. It creates a fictional dataset for
study with the same design as Peter's (one fixed effect, two crossed random
effects). I call the fixed effect "A" and the two random effects, "B" and
"C". To translate the do-file to his study, he can think of these as A =
modality, B = cases, C = readers, and response = lntime. I put in a
positive variance component for the interaction of the two random effects,
but left the random slopes (A-by-B and A-by-C) with whatever they got by
chance. As you'd guess, variance components for both are miniscule. As it
turned out, one of them is actually negative and miniscule; the other is
just miniscule. This is actually the opposite of what Peter can expect
with his real data: modality-by-case and modality-by-reader undoubtedly are
positive and might even be discriminable from zero in his study, but it's
likely that case-by-reader would be more difficult to distinguish from zero.
(Exceptional diagnosticians unfazed by even the most difficult cases,
perhaps?)
As expected, -xtmixed- didn't converge and left us with the last
expectation-maximization iteration. Dropping the zero variance components
allowed convergence. Of course, -anova- didn't balk and allowed testing and
estimation of the variance components. There's good agreement in estimates
of variance components between -anova- and -xtmixed-, so it looks as if I
haven't riddled the formulas for expected mean squares with mistakes. (This
fictional example isn't really good to illustrate the likelihood ratio test
for a variance component, because we didn't get convergence with the full
model.)
Joseph Coveney
clear
set more off
set seed `=date("2005-08-02", "ymd")'
tempfile tmpfil0
set obs 20
generate byte C = _n
generate float Ceffect = invnorm(uniform())
save `tmpfil0', replace
clear
set obs 20
generate byte B = _n
generate float Beffect = invnorm(uniform())
cross using `tmpfil0'
erase `tmpfil0'
generate float BCeffect = invnorm(uniform())
forvalues i = 1/3 {
generate float response`i' = `i' + ///
2 * Beffect + 4 * Ceffect + BCeffect + ///
invnorm(uniform()) / 2
}
reshape long response, i(B C) j(A)
anova response A B C A*B A*C B*C, category(A B C)
local MSE = e(rss) / e(df_r)
local MSB = e(ss_2) / e(df_2)
local MSC = e(ss_3) / e(df_3)
local MSAB = e(ss_4) / e(df_4)
local MSAC = e(ss_5) / e(df_5)
local MSBC = e(ss_6) / e(df_6)
display "Variance Component B = " ///
(`MSB' - (e(df_1) + 1) * `MSBC' - `MSE') / ///
((e(df_1) + 1) * (e(df_3) + 1))
display "Variance Component C = " ///
(`MSC' - (e(df_1) + 1) * `MSBC' - `MSE') / ///
((e(df_1) + 1) * (e(df_2) + 1))
display "Variance Component AB = " ///
(`MSAB' - `MSE') / (e(df_3) + 1)
display "Variance Component AC = " ///
(`MSAC' - `MSE') / (e(df_2) + 1)
display "Variance Component BC = " ///
(`MSBC' - `MSE') / (e(df_1) + 1)
egen long AXB = group(A B)
egen long AXC = group(A C)
egen long BXC = group(B C)
compress
xi: xtmixed response i.A || ///
_all: R.B || _all: R.C || ///
_all: R.AXB || _all: R.AXC || ///
_all: R.BXC, nolrtest variance
estimates store A
xi: xtmixed response i.A || ///
_all: R.B || _all: R.C || ///
_all: R.BXC, nolrtest variance
lrtest A .
exit
*
* 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/
