[R] anova.mlm for single model (one-way repeated measured anova)

Prof Brian Ripley ripley at stats.ox.ac.uk
Sat Aug 12 08:59:22 CEST 2006


On Sat, 12 Aug 2006, takahashi kohske wrote:

> Dear list members:
> 
> I'd like to one-way repeated measured anova by using mlm.
> I'm using R-2.3.1 and my code is:
> 
> dat<-matrix( c(9,7,8,8,12,11,8,13,     6,5,6,3,6,7,10,9,
>                10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14),
>             ncol=4, dimname=list(s=1:8, c=1:4))
> mlmfit<-lm(dat~1)
> anova(mlmfit, X=~1)
> Error: ceiling(length.out) : Non-numeric argument to mathematical function
> 
> this error occurs in anova.mlm
> 
>         if (rk > 0) {
>             p1 <- 1:rk
>             comp <- object$effects[p1, ]
>             asgn <- object$assign[object$qr$pivot][p1]
>             nmeffects <- c("(Intercept)", attr(object$terms,
>                 "term.labels"))
>             tlabels <- nmeffects[1 + unique(asgn)]
>              ix <- split(seq(length = nrow(comp)), asgn)  #HERE
>             ss <- lapply(ix, function(i) crossprod(comp[i, ,
>                 drop = FALSE]))
>             df <- sapply(split(asgn, asgn), length)
>         }
> 
> because nrow(comp) returns NULL.
> 
> in my memory, R-2.2.* ( or may be R-2.3.0) can correctly handle this code.
> so, I think this is a kind of side-effect of fixing PR#8679.
> 
> currently, i can workaround as follows:
> 
> anova(mlmfit, update(mlmfit, ~0), X=~1)
> 
> this code returns correct answer.
> 
> 
> I don't know whether this behavior is correct or bug.


Yes, it is a bug.  The line

            comp <- object$effects[p1, ]

should be

            comp <- object$effects[p1, , drop=FALSE]

I am changing this in 2.3.1 patched and R-devel.


-- 
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