[Rd] wishlist: function mlh.mlm to test multivariate linear hypotheses of the form: LBT'=0 (PR#8680)
Yves.Rosseel@UGent.be
Yves.Rosseel at UGent.be
Mon Mar 13 18:09:19 CET 2006
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
}
More information about the R-devel
mailing list