[EMAIL PROTECTED] 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
~~~~~~~~~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
