[R] nlme vs aov with Error() for an ANCOVA
Prof Brian Ripley
ripley at stats.ox.ac.uk
Thu Jan 15 15:33:11 CET 2004
Well, I don't think this is ANCOVA as you seem to want to specify a random
slope for a covariate. aov() is not designed for that. It is also not
designed for assessing the size of fixed effects which seems the question
here.
As I understand it, you have only one observation for each value of `rt'
for each subject, and `rt' is an explanatory variable. For lme you have
specified a subject-dependent intercept and coefficient of `rt'. You
cannot do that in aov, where the argument of Error is supposed to be a
factor or a combination of factors. This is in the reference given by
help for aov or Error (on pp 157-9).
On Thu, 15 Jan 2004, Christoph Lehmann wrote:
> Hi
> I compouted a multiple linear regression with repeated measures on one
> explanatory variable:
> BOLD peak (blood oxygenation) as dependent variable,
>
> and as independent variables I have:
> -age.group (binaray:young(0)/old(1))
> -and task-difficulty measured by means of the reaction-time 'rt'. For
> 'rt' I have repeated measurements, since each subject did 12 different
> tasks.
> -> so it can be seen as an ANCOVA
>
> subject age.group bold rt
>
> subj1 0 0.08 0.234
> subj1 0 0.05 0.124
> ..
> subj1 0 0.07 0.743
>
> subj2 0 0.06 0.234
> subj2 0 0.02 0.183
> ..
> subj2 0 0.05 0.532
>
> subjn 1 0.09 0.234
> subjn 1 0.06 0.155
> ..
> subjn 1 0.07 0.632
>
> I decided to use the nlme library:
>
> patrizia.lme <- lme(bold ~ rt*age.group, data=patrizia.data1, random= ~
> rt |subject)
> > summary(patrizia.lme)
> Linear mixed-effects model fit by REML
> Data: patrizia.data1
> AIC BIC logLik
> 272.2949 308.3650 -128.1474
>
> Random effects:
> Formula: ~rt | subject
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> (Intercept) 0.2740019518 (Intr)
> rt 0.0004756026 -0.762
> Residual 0.2450787149
>
> Fixed effects: bold ~ rt + age.group + rt:age.group
> Value Std.Error DF t-value p-value
> (Intercept) 0.06109373 0.11725208 628 0.521046 0.6025
> rt 0.00110117 0.00015732 628 6.999501 0.0000
> age.group -0.03750787 0.13732793 43 -0.273126 0.7861
> rt:age.group -0.00031919 0.00018259 628 -1.748115 0.0809
> Correlation:
> (Intr) rt ag.grp
> rt -0.818
> age.group -0.854 0.698
> rt:age.group 0.705 -0.862 -0.805
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -3.6110596 -0.5982741 -0.0408144 0.5617381 4.8648242
>
> Number of Observations: 675
> Number of Groups: 45
>
> --end output
> #-> if the model assumptions hold this means, we don't have a
> significant age effect but a highly significant task-effect and the
> interaction is significant on the 0.1 niveau.
Nope. It means that you have two lines with a common non-zero intercept
and probably different slopes for the two age groups. However, as I
understand it rt=0 is an extraplolation to an physically impossible value,
so interpreting the intercept makes little sense.
> I am now interested, if one could do the analysis also using aov and the
> Error() option.
>
> e.g. may I do:
> > l <- aov(bold ~ rt*age.group + Error(subject/rt),data=patrizia.data1)
> > summary(l)
>
> Error: subject
> Df Sum Sq Mean Sq
> rt 1 0.0022087 0.0022087
>
> Error: subject:rt
> Df Sum Sq Mean Sq
> rt 1 40.706 40.706
>
> Error: Within
> Df Sum Sq Mean Sq F value Pr(>F)
> rt 1 2.422 2.422 10.0508 0.001592 **
> age.group 1 8.722 8.722 36.2022 2.929e-09 ***
> rt:age.group 1 0.277 0.277 1.1494 0.284060
> Residuals 669 161.187 0.241
> ---
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
>
>
> which looks weird
>
> or what would you recommend?
>
> thanks a lot
>
> Christoph
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-help
mailing list