[R] Manova and contrasts
Yves Rosseel
Yves.Rosseel at UGent.be
Wed Jun 2 17:56:00 CEST 2004
Dear Eduardo,
> Hi R-users
> I'm trying to do multivariate analysis of variance of a experiment with
> 3 treatments, 2 variables and 5 replicates.
> The procedure adopted in SAS is as follow, but I'm having difficulty in
> to implement the contrasts for comparison of all treatments in R.
[See the original mail for the SAS input and output...]
What you probably want is a general function to test multivariate linear
hypotheses of the type 'LBM=K' with B the parameter matrix.
Unfortunately, as far as I know, no such a function is available in R or
any package on CRAN. (For the univariate case, there are some functions
for testing linear hypotheses in the packages car and gregmisc) It is,
however, fairly easy to write such a function yourself. All you need to
do is to compute the corresponding SSCP matrix for your particular set
of contrasts, and 'borrow' some code from the 'summary.manova' function
(only Wilks lambda is computed here):
mlh <- function(fit, L, M) {
if(!inherits(fit, "maov"))
stop("object must be of class \"manova\" or \"maov\"")
if( is.null(dim(L)) )
L <- t(L)
rss.qr <- qr( crossprod( fit$residuals %*% M ) )
X <- as.matrix( model.matrix(fit) )
B <- as.matrix( fit$coef )
H <- t(M) %*% t(L%*%B) %*%
as.matrix(solve(L%*%solve(t(X)%*%X)%*%t(L))) %*%
(L%*%B)%*%M
eig <- Re(eigen(qr.coef(rss.qr, H), symmetric = FALSE)$values)
q <- nrow(L); df.res <- fit$df.residual
test <- prod(1/(1 + eig))
p <- length(eig)
tmp1 <- df.res - 0.5 * (p - q + 1)
tmp2 <- (p * q - 2)/4
tmp3 <- p^2 + q^2 - 5
tmp3 <- if(tmp3 > 0) sqrt(((p*q)^2 - 4)/tmp3) else 1
wilks <- test
F <- ((test^(-1/tmp3) - 1) * (tmp1 * tmp3 - 2 * tmp2))/p/q
df1 <- p * q
df2 <- tmp1 * tmp3 - 2 * tmp2
Prob <- pf(F, df1, df2, lower.tail = FALSE)
out <- list(wilks=wilks, F=F, df1=df1, df2=df2, Prob=Prob)
out
}
To test the first contrast in your example (contrast 'TESTE vs TURFE'
TRA 1 -1 0), you can proceed as follows:
# your data (copied from your original mail)
### R
X1 =
c(4.63,4.38,4.94,4.96,4.48,6.03,5.96,6.16,6.33,6.08,4.71,4.81,4.49,4.43,
4.56)
X2 =
c(0.95,0.89,1.01,1.23,0.94,1.08,1.05,1.08,1.19,1.08,0.96,0.93,0.87,0.82,
0.91)
Trat = as.factor(c(rep("TESTE",5),rep("TURFE",5), rep("TURNA",5)))
Y = cbind(X1,X2)
fit = manova(Y ~ Trat)
# construct a 'L-matrix' (in this case only one row)
L <- rbind( c(0, -1, 0) )
# note: your contrast is '1 -1 0'
# adding a zero for the intercept gives
# L = (0 1 -1 0)
# removing 'redundant' coefficients gives
# L = (0 -1 0)
# The coefficient '1' should not be
# specified because the corresponding parameters were fixed
# to zero, and (unlike SAS or SPSS!) these
# redundant parameters are not included in the parameter matrix
# M-matrix is 2x2 identity matrix
M <- diag(2)
# test multivariate linear hypothesis
mlh(fit, L, M)
# the output is the same as in your SAS output:
$wilks
[1] 0.03437322
$F
[1] 154.5083
$df1
[1] 2
$df2
[1] 11
$Prob
[1] 8.896343e-09
Since multivariate linear hypotheses are widely used in psychology and
related disciplines, it would be nice if a more robust and cleaned-up
version of the 'mlh' function would eventually be added to R...
Yves Rosseel.
--
Dr. Yves Rosseel -- http://allserv.ugent.be/~yrosseel/
Department of Data Analysis, Ghent University
Henri Dunantlaan 1, B-9000 Gent, Belgium
More information about the R-help
mailing list