[Rd] wishlist: function mlh.mlm to test multivariate linear hypotheses of the form: LBT'=0 (PR#8680)

[Yves.Rosseel@UGent.be]

Full_Name: Yves Rosseel
Version: 2.2.1
OS: 
Submission from: (NULL) (157.193.116.152)


The code below sketches a possible implementation of a function 'mlh.mlm' which
I think would be a good complement to the 'anova.mlm' function in the stats
package. It tests a single linear hypothesis of the form H_0: LBT'= 0 where B is
the matrix of regression coefficients; L is a matrix with rows giving linear
combinations of the regressions coefficients; the transformation matrix T has
the same meaning as in the anova.mlm function. An example and some bare-bones
code is listed below (code depends on the Pillai, Wilks etc. functions defined
in src/library/stats/R/mlm.R).

Example model: 3 dependents, 1 between-subjects factor with 4 levels

set.seed(123)
Y <- cbind(rnorm(100), rnorm(100), rnorm(100)) 
A <- factor(rep(c(1,2,3,4), each=25))
fit <- lm(Y ~ A)

Example 1: simple contrast: compare level 3 versus level 4 (multivariate)
(note: first zero in l1 corresponds to the intercept, not the first level of A)
l1 <- c(0, 0, 1, -1)
L <- rbind(l1)
mlh.mlm(fit, L=L, test="Wilks")

     Wilks   approx F     num Df     den Df     Pr(>F)
 0.9874192  0.3992218  3.0000000 94.0000000  0.7538689

Example 2: suppose the three dependents are three timepoints (time); is there a
contrast*time interaction (using the contrast above: level 3 versus level 4)

t1 <- c(1,-1,0); t2 <- c(0,1,-1)
T <- rbind(t1,t2)
mlh.mlm(fit, L=L, T=T, test="Wilks")

     Wilks   approx F     num Df     den Df     Pr(>F)
 0.9889929  0.5286555  2.0000000 95.0000000  0.5911205



Code:
------------------------------------------------------------------------------------
mlh.mlm <-
    function(object, L = null, T = diag(nrow = p),
             test = c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"))
{
    # test the multivariate linear hypothesis LBT'=0
    # B = matrix of regression coefficients
    # L = matrix, each row is a linear combination of the parameters
    # T = transformation matrix

    if(!inherits(object, "mlm"))
        stop("object must be of class \"mlm\"")

    if( is.null(L) )
        stop("L matrix is not specified.")

    # L must be a matrix
    if( is.null(dim(L)) )
        L <- t(L)

    if( nrow(object$coef) != ncol(L) )
        stop("nrow(object$coef) != ncol(L)")

    p <- ncol(SSD(object)$SSD)
    ssd <- SSD(object)
    df.res <- ssd$df
    rss.qr <- qr(T %*% ssd$SSD %*% t(T))

    X <- as.matrix( model.matrix(object) )
    B <- as.matrix( object$coef )

    df <- nrow(L)
    ss <- t(L %*% B) %*%
          as.matrix(solve(L %*% solve(t(X) %*% X) %*% t(L))) %*%
          (L %*% B)

    eigs <- Re(eigen(qr.coef(rss.qr,
                             T %*% ss %*% t(T)),
                             symmetric = FALSE)$values)

    test <- match.arg(test)
    stats <- switch(test,
                           "Pillai" = Pillai(eigs, df, df.res),
                           "Wilks"  =  Wilks(eigs, df, df.res),
                           "Hotelling-Lawley" = HL(eigs, df, df.res),
                           "Roy"    = Roy(eigs, df, df.res)
                    )
    stats[5] <- pf(stats[2], stats[3], stats[4], lower.tail = FALSE)
    names(stats) <- c(test, "approx F", "num Df", "den Df", "Pr(>F)")

    stats
}