From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>

Date: Fri 18 Feb 2005 - 21:57:08 EST

Date: Fri 18 Feb 2005 - 21:57:08 EST

Dear R-devel'ers

Below is an outline for a set of routines to improve support for multivariate linear models and "classical" repeated measurements analysis. Nothing has been coded yet, so everything is subject to change as loose ideas get confronted by the harsh realities of programming.

Comments are welcome. They might even influence the implementation...

-pd

General considerations:

- S3 class based to fit existing code
- similar to lm/glm code

fit <- lm(Y~...) creates basic mlm object (already does)

SSD(obj) returns object of class "SSD":

$SSD matrix of sums of squares & products

$df degrees of freedom.

Methods for (lm and) mlm.

estVar(obj) obj$SSD/obj$df (could have methods for lm/mlm too)

Convenience functions:

Tr is the trace operator sum(diag(M))

proj is the projection operator possibly generalized to matrices.
Rg: matrix rank (not sure we really need it, but see below)

mauchley.test(obj, M=diag(ncol=p), T = proj(X, orth=TRUE),

X = matrix(rep(1,p)))

(p = ncol(obj$SSD))

Test of sphericity, i.e. that the obj represents a empirical covariance matrix S with true value proportional to M or that TST' is proportional to TMT'. Alternatively, give X with the property that TX == 0. (One sticky bit is that you can't really just use proj() because T must have maximal rank. What is the current best practice for dealing with that? qr() pivoting?)

summary.mlm

print.summary.mlm

vcov.mlm

summary.mlm could be a little smarter than just coordinatewise summary.lm. It could at least provide the estimated residual covariance matrix (or SSD structure).

vcov currently inherits from "lm" leading to a completely arbitrarily scaled matrix. The correct matrix is a Kronecker product of the unscaled covariance matrix and estVar.

anova.mlm, anova.mlmlist, drop1.mlm, add1.mlm

These can (seemingly...) be obtained by relatively small modifications of their lm counterparts. The actual test calculations need to be excised from summary.manova (generalize? e.g. mvlin.test(SSD1, SSD2, method="Pillai")). It should be possible to wedge in tests under sphericity assumptions (with Greenhouse-Geisser and Huynh-Feldt corrections), as well as transformation/conditioning matrices (see below).

A more radical idea is to say that these are all different kinds of MANOVA tables and extend the *.manova functions to understand them. This would break the symmetry with lm, though, since you then need to use summary() to get a meaningful listing.

Notes on conditional tests (lower priority):

Consider Y = (Y1,Y2) and a linear hypothesis matrix H.

The test of zero intercept in the (multivariate) regression of H %*% Y2 on H %*% Y1 is of some interest, possibly after a linear transformation of Y. This is the test for additional information, but it is also the correct (ML) way of utilizing known-zero effects, e.g. pretreatment measurements. There is even a neat trick to fitting a linear structure across the responses by regressing on null-space contrasts.

To some extent, conditional tests can be handled just by moving variables to the r.h.s. of the linear model specification, but there might be a point in having a more evocative interface, especially where transformed Y's are involved. This could be formula-based or matrix-based: contrasts=ginv(contr.sdif(4)) or formula based: ystruct=~index.

-- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 ______________________________________________ R-devel@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-develReceived on Fri Feb 18 21:31:25 2005

*
This archive was generated by hypermail 2.1.8
: Fri 18 Feb 2005 - 23:30:13 EST
*