[Rd] Suggestions for enhanced routines for "mlm" models.

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
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...


General considerations:

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 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
Received 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