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.
I have already read manuals and other materials about manova in R, but
nothing about specific contrasts were found in them, just regarding buit
contrasts, such as contr.helmert, contr.poly, contr.treatment and
contr.sum. How to get output similar to that from SAS, mainly for the
three first contrasts, once that fourth contrast is contr.helmert?
### SAS
INPUT
proc glm;
classes TRA;
model X1 X2 = TRA;
contrast 'TESTE vs TURFE' TRA 1 -1 0;
contrast 'TESTE vs TURNA' TRA 1 0 -1;
contrast 'TURFE vs TURNA' TRA 0 1 -1;
contrast 'Ho:TESTE=TURFE=TURNA' TRA 1 1 -2,
TRA 1 -1 0;
manova h=TRA / printe printh;
run;
OUTPUT (some parts were removed)
...
Dependent Variable: X1
Source DF Sum of Squares Mean
Square F Value Pr > F
Model 2 7.24764000
3.62382000 94.96 0.0001
Error 12 0.45796000
0.03816333
Corrected Total 14 7.70560000
Source DF Type I SS Mean
Square F Value Pr > F
TRA 2 7.24764000
3.62382000 94.96 0.0001
Source DF Type III SS Mean
Square F Value Pr > F
TRA 2 7.24764000
3.62382000 94.96 0.0001
Contrast DF Contrast SS Mean
Square F Value Pr > F
TESTEMUNHA vs TURFE 1 5.14089000
5.14089000 134.71 0.0001
TESTEMUNHA vs TURNA 1 0.01521000
0.01521000 0.40 0.5397
TURFE vs TURNA 1 5.71536000
5.71536000 149.76 0.0001
Ho:TESTE=TURFE=TURNA 2 7.24764000
3.62382000 94.96 0.0001
...
Dependent Variable: X2
Source DF Sum of Squares Mean
Square F Value Pr > F
Model 2 0.09817333
0.04908667 6.22 0.0140
Error 12 0.09472000
0.00789333
Corrected Total 14 0.19289333
Source DF Type I SS Mean
Square F Value Pr > F
TRA 2 0.09817333
0.04908667 6.22 0.0140
Source DF Type III SS Mean
Square F Value Pr > F
TRA 2 0.09817333
0.04908667 6.22 0.0140
Contrast DF Contrast SS Mean
Square F Value Pr > F
TESTEMUNHA vs TURFE 1 0.02116000
0.02116000 2.68 0.1275
TESTEMUNHA vs TURNA 1 0.02809000
0.02809000 3.56 0.0837
TURFE vs TURNA 1 0.09801000
0.09801000 12.42 0.0042
Ho:TESTE=TURFE=TURNA 2 0.09817333
0.04908667 6.22 0.0140
...
Manova Test Criteria and F Approximations for the Hypothesis of no
Overall TRA Effect
H = Type III SS&CP Matrix for TRA E = Error SS&CP
Matrix
S=2 M=-0.5 N=4.5
Statistic Value F Num DF
Den DF Pr > F
Wilks' Lambda 0.02042803 32.9813 4
22 0.0001
Pillai's Trace 1.24157026 9.8222 4
24 0.0001
Hotelling-Lawley Trace 35.12690920 87.8173 4
20 0.0001
Roy's Greatest Root 34.75791615 208.5475 2
12 0.0001
...
Manova Test Criteria and Exact F Statistics for
the Hypothesis of no Overall TESTEMUNHA vs TURFE
Effect
H = Contrast SS&CP Matrix for TESTEMUNHA vs TURFE E Error
SS&CP Matrix
S=1 M=0 N=4.5
Statistic Value F Num DF
Den DF Pr > F
Wilks' Lambda 0.03437322 154.5083 2
11 0.0001
Pillai's Trace 0.96562678 154.5083 2
11 0.0001
Hotelling-Lawley Trace 28.09241176 154.5083 2
11 0.0001
Roy's Greatest Root 28.09241176 154.5083 2
11 0.0001
...
Manova Test Criteria and Exact F Statistics for
the Hypothesis of no Overall TESTEMUNHA vs TURNA
Effect
H = Contrast SS&CP Matrix for TESTEMUNHA vs TURNA E Error
SS&CP Matrix
S=1 M=0 N=4.5
Statistic Value F Num DF
Den DF Pr > F
Wilks' Lambda 0.65512498 2.8953 2
11 0.0977
Pillai's Trace 0.34487502 2.8953 2
11 0.0977
Hotelling-Lawley Trace 0.52642630 2.8953 2
11 0.0977
Roy's Greatest Root 0.52642630 2.8953 2
11 0.0977
...
Manova Test Criteria and Exact F Statistics for
the Hypothesis of no Overall TURFE vs TURNA
Effect
H = Contrast SS&CP Matrix for TURFE vs TURNA E = Error
SS&CP Matrix
S=1 M=0 N=4.5
Statistic Value F Num DF
Den DF Pr > F
Wilks' Lambda 0.03988589 132.3934 2
11 0.0001
Pillai's Trace 0.96011411 132.3934 2
11 0.0001
Hotelling-Lawley Trace 24.07152574 132.3934 2
11 0.0001
Roy's Greatest Root 24.07152574 132.3934 2
11 0.0001
...
Manova Test Criteria and F Approximations for
the Hypothesis of no Overall Ho:TESTE=TURFE=TURNA
Effect
H = Contrast SS&CP Matrix for Ho:TESTE=TURFE=TURNA E Error
SS&CP Matrix
S=2 M=-0.5 N=4.5
Statistic Value F Num DF
Den DF Pr > F
Wilks' Lambda 0.02042803 32.9813 4
22 0.0001
Pillai's Trace 1.24157026 9.8222 4
24 0.0001
Hotelling-Lawley Trace 35.12690920 87.8173 4
20 0.0001
Roy's Greatest Root 34.75791615 208.5475 2
12 0.0001
### 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)
summary.aov(fit)
summary(fit, test= "Wilks")
ANOVA:
Response X1 :
Df Sum Sq Mean Sq F value Pr(>F)
Trat 2 7.2476 3.6238 94.956 4.407e-08 ***
Residuals 12 0.4580 0.0382
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
Response X2 :
Df Sum Sq Mean Sq F value Pr(>F)
Trat 2 0.098173 0.049087 6.2187 0.01402 *
Residuals 12 0.094720 0.007893
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
MANOVA:
Df Wilks approx F num Df den Df Pr(>F)
Trat 2 0.020 32.981 4 22 5.302e-09 ***
Residuals 12
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
Any suggestions are much appreciated.
__________________________________________________________
Eng. Agr., M.Sc. Eduardo Dutra de Armas
__________________________________________________________
Centro de Energia Nuclear na Agricultura (CENA/USP)
Laboratório de Ecotoxicologia
Av.Centenário 303, C.P. 96, CEP 13400-970, Piracicaba, SP, Brasil
__________________________________________________________
---
[[alternative HTML version deleted]]
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