[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
Joseph Coveney <jcoveney@bigplanet.com> |

To |
Statalist <statalist@hsphsun2.harvard.edu> |

Subject |
Re: re: st: Specifying a repeated measures analysis in anova andxtmixed |

Date |
Wed, 03 Aug 2005 18:23:58 +0900 |

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/

- Prev by Date:
**Re: st: output format** - Next by Date:
**st: Test for cross-sectional correlation with small T** - Previous by thread:
**Re: re: st: Specifying a repeated measures analysis in anova andxtmixed** - Next by thread:
**st: nbreg with cluster gives NO LR test** - Index(es):

© Copyright 1996–2015 StataCorp LP | Terms of use | Privacy | Contact us | What's new | Site index |