SR Millis wrote:
I've been using SAS proc mixed to analyze AB/AB
cross-over designs (both mixed and fixed effects
models). Out of curiosity, I compared results with
Stata's pkcross--using the Chow & Liu data (used in
the Reference manual example).
I'm getting very different results in SAS vs Stata,
eg, SAS is finding a significant treatment effect
while Stata isn't. What could be accounting for the
differences?
--------------------------------------------------------------------------------
Perhaps if you posted your SAS code for this example, someone might be able to
assist. The dataset is balanced, so PROC MIXED (with the default REML method)
should be giving you the same result as PROC GLM and PROC ANOVA.
Just by inspection, there doesn't seem to be a treatment effect.
Joseph Coveney
clear
set more off
use http://www.stata-press.com/data/r8/chowliu.dta
* Pooling for inspection
generate delta = period2 - period1
replace delta = -delta if seq == 2
ttest delta = 0, level(90) // P = 0.54, when pooled
bysort seq: ttest delta = 0, level(90) // Nothing hiding here
drop delta
*
pkshape id seq period1 period2, order(ab ba)
* Inspection redux--compare treat 1 period 1 with treat 2 period 2
* and treat 2 period 1 with treat 1 period 2
table treat period, contents(mean outcome sd outcome n outcome) format(%3.0f)
*
forvalues parameterization = 1/4 {
pkcross outcome, param(`parameterization') sequence(sequence) ///
treatment(treat) carryover(carry) period(period) id(id)
}
* Using the third parameterization,
anova outcome sequence / id|sequence treat treat*sequence
* Or, equivalently,
anova outcome sequence / id|sequence treat period
lincom _b[treat[2]] - _b[treat[1]], level(90)
* Random effects, FGLS (REML, here)
xi: xtreg outcome i.sequence i.period i.treat, i(id) re sa theta
adjust _Isequence_2 = 0 _Iperiod_2 = 0, by( _Itreat_2) xb se
lincom _b[_Itreat_2], level(90) // Note z = t above; use as
// t test statistic (with df), and P will be the same
* Random effects, full maximum likelihood
xi: xtreg outcome i.sequence i.period i.treat, i(id) mle nolog
estimates store A
adjust _Isequence_2 = 0 _Iperiod_2 = 0, by( _Itreat_2) xb se
quietly xi: xtreg outcome i.sequence i.period, i(id) mle nolog
lrtest A ., stats // Slightly anticonservative, as expected, but . . .
*
estimates drop _all
keep id sequence outcome treat
reshape wide outcome, i(id) j(treat)
xi: mvreg outcome1 outcome2 = i.sequence
lincom ([outcome2]_cons + 0.5 * [outcome2]_Isequence_2) - ///
([outcome1]_cons + 0.5 * [outcome1]_Isequence_2), level(90) // Same
* Equivalently,
quietly manova outcome1 outcome2 = sequence
matrix Representative_sequence = (1, 0.5, 0.5)
matrix Treatment_difference = (1, -1)
manovatest , test(Representative_sequence) ytransform(Treatment_difference)
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/
