[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