[R] means over factors in mlm terms
Michael Friendly
friendly at yorku.ca
Thu Nov 23 18:55:49 CET 2006
Thanks, Richard
Years ago, I recall reading a tech report on an elegant system
in APL for doing (M)ANOVA via a series of successive projections of
the response(s) on the series of (successively orthogonalized)
basis vectors of the design matrix. (Written by RMH) Perhaps this was
what you were thinking of.
The algorithm (in R-aplish) was quite lovely:
aov <- function ( Y, class, design ) {
err <- Y
fit <- NULL
X <- J(nrow(Y), 1) # intercept
k <- 1
loop: Xk <- xbasis( design[k,], class)
Xk <- Xk - proj(Xk, X) # orthogonalize to previous effects
X <- cbind(X, Xk) # append to what we've fitted
fitk <- proj(err, Xk) # project remaining Y on Xk
fit <- cbind(fit, fitk) # append to what we've fitted
err <- err - fit # and remove from Y
-> loop if nrow(design) < k <- k + 1
done:
SS <- t(fit) %*% fit
SSE <- t(err) %*% err
}
here, Y (n x p) and class (f x p) are the data and class/factor
variables, design (#terms x f) is a binary matrix whose rows
specify the factors that participate in each effect, e.g.,
design = ( 1 0 / 0 1 / 1 1 ) <--> y ~ A + B + A:B.
One nice thing was that projecting the current residuals from Y
on the basis vectors for a given X[k,] term gave the fitted values
for that effect, whose SS were the (type 1) SS for that effect. Another
was that it generalized automatically to the mlm, where Y'Y gives
SSCP matrices. It would be nice if someone (you?) developed this
for from aov() and model.tables() for R, where mlm methods are still
somewhat limited.
-Michael
Richard M. Heiberger wrote:
> ## Thank you, Michael, for the example. I believe you are looking
> ## for the natural extension to mlm of the proj function.
> ## Here is the lapply workaround from which that extension
> ## might be written
> ##
> ## Rich
>
> sapply(soils, class)
> soils$Group <- factor(soils$Group)
> soils$Block <- factor(soils$Block)
> sapply(soils, class)
>
> result <- list()
> for (i in names(soils)[6:14])
> result[[i]] <- aov(soils[[i]] ~ Block + Contour*Depth, data=soils)
>
> proj(result[[1]])
>
> lapply(result, proj)
>
> model.tables(result[[1]], type="means")
> lapply(result, model.tables, type="means")
>
> ## Exercise for the reader
> ## 1. defined proj.lm to be proj.aov when it makes sense
> ## 2. extend the above to the proj.mlm method
> ## 3. collapse the projections to tables
> ## (this is actually how model.tables is written).
> ## 4. extend model.tables to the lm and mlm when it makes sense.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT M3J 1P3 CANADA
More information about the R-help
mailing list