[R] Manova and contrasts

[Yves Rosseel]

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