R devel - Mar 2006 - anova.mlm (single-model case) does not handle factors? (PR#8679)

Full_Name: Yves Rosseel
Version: 2.2.1
OS: i686-pc-linux-gnu
Submission from: (NULL) (157.193.116.152)


Dear developers,

For the single-model case, the anova.mlm() function does not seem to handle
multi-parameter predictors (eg factors) correctly. A toy example illustrates the
problem:

Y <- cbind(rnorm(100),rnorm(100),rnorm(100))
A <- factor(rep(c(1,2,3,4), each=25))
fit <- lm(Y ~ A)
anova.mlm(fit)

gives

Error in T %*% ss[[i]] : non-conformable arguments

I think the problem lies in the computation of the 'ss' terms (line 237
in the
file mlm.R in the source code). Changing this (and the following line) by
something similar to what is done in summary.manova seems to resolve the
problem:

 comp <- as.matrix(fit$effects)[seq(along=asgn), ,drop=FALSE]
 uasgn <- unique(asgn)
 nterms <- length(uasgn)
 ss <- list(nterms)
 df <- numeric(nterms)
 for(i in seq(nterms)) {
    ai <- (asgn == uasgn[i])
    ss[[i]] <- crossprod(comp[ai, ,drop=FALSE])
    df[i] <- sum(ai)
 }

Yves.Rosseel at ugent.be writes:
> Full_Name: Yves Rosseel
> Version: 2.2.1
> OS: i686-pc-linux-gnu
> Submission from: (NULL) (157.193.116.152)
> 
> 
> Dear developers,
> 
> For the single-model case, the anova.mlm() function does not seem to handle
> multi-parameter predictors (eg factors) correctly. A toy example
illustrates the
> problem:
> 
> Y <- cbind(rnorm(100),rnorm(100),rnorm(100))
> A <- factor(rep(c(1,2,3,4), each=25))
> fit <- lm(Y ~ A)
> anova.mlm(fit)
> 
> gives
> 
> Error in T %*% ss[[i]] : non-conformable arguments
> 
> I think the problem lies in the computation of the 'ss' terms (line
237 in the
> file mlm.R in the source code). Changing this (and the following line) by
> something similar to what is done in summary.manova seems to resolve the
> problem:
> 
>  comp <- as.matrix(fit$effects)[seq(along=asgn), ,drop=FALSE]
>  uasgn <- unique(asgn)
>  nterms <- length(uasgn)
>  ss <- list(nterms)
>  df <- numeric(nterms)
>  for(i in seq(nterms)) {
>     ai <- (asgn == uasgn[i])
>     ss[[i]] <- crossprod(comp[ai, ,drop=FALSE])
>     df[i] <- sum(ai)
>  }
Thanks for pointing this out. Amazing that it hasn't cropped up before...

The culprit seems to be the line

ss <- lapply(split(comp, asgn), function(x) crossprod(t(x)))

in which "comp" (an effects by responses matrix) is effectively turned
into a vector by split(). And, adding insult to injury, had you
actually gotten a matrix result, you wouldn't want to transpose it
before calculating the cross products.

A simpler, but slightly dirty fix is to set 

ss <- lapply(split.data.frame(comp, asgn), crossprod)

(the dirtiness is that split.data.frame actually does the right thing
for matrices; possibly, an identically defined split.matrix() would be
cleaner, although I'm never quite happy with matrix functions that
aren't symmetric in rows/columns.)

-- 
   O__  ---- Peter Dalgaard             ?ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907