[R] means over factors in mlm terms

[Michael Friendly]

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