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@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 ______________________________________________ R-help@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.