﻿## One-Way MANOVA ##---------------------------------------------------------------------------------------------------------------------------------------------- #Example 1 #Data on student scores in aptitude, mathematics, language and general knowledge for students in technical #disciplines #(group 1), architecture (group 2) and medical technology students (group 3) as indicated in column 1. morel <- read.table(file = "datafiles/morel.dat", col.names=c("studentgroup", "aptitude", "mathematics", "language", "generalknowledge")) morel$studentgroup <- as.factor(morel$studentgroup) fit.lm <- lm(cbind(aptitude, mathematics, language, generalknowledge)~studentgroup , data = morel) # Note that by default, R puts the first treatment tau1'=(tau11,tau12,tau13,tau14)=0 fit.lm summary(fit.lm) #y_{ij}=mu_i+e_{ij}=mu+tau_i+e_{ij}: mu1=mu, mu2=mu+tau2, mu3=mu+tau3 # If we want to do the same as SAS, we have to use the contrasts option: #SAS puts the last treastment tau3'=(tau31,tau32,tau33,tau34)=0 fit.lm1 <- lm(cbind(aptitude, mathematics, language, generalknowledge)~studentgroup , data = morel, contrasts = list(studentgroup = contr.SAS)) fit.lm1 summary(fit.lm1) #H0:tau1=tau2=tau3=0 library(car) #Manova by using car #MANOVA fit.manova <- Manova(fit.lm) summary(fit.manova) #manova by using stats fit.anova<-manova(fit.lm) summary(fit.anova,test="Wilks") summary(fit.anova,test="Pillai") # Bonferroni simultaneous confidence intervels for treatments' differences trteff.diff(cbind(morel[,-1],morel[,1])) ##Test other hypothesis ##H0:Cbeta=0 C <- matrix(c(0, 1, 0, 0, 1, -1), ncol = 3, by = TRUE) newfit <- linearHypothesis(model = fit.lm, hypothesis.matrix = C) print(newfit) #H0: C beta M=0 M <- matrix(c(1, 0, 0, -1, 1, 0, 0, -1, 1, 0, 0, -1), nrow = 4, by = TRUE) newfit <- linearHypothesis(model = fit.lm, hypothesis.matrix = C, P = M) print(newfit) #for comparison fitp2<-lm((as.matrix(morel[,-1])%*%M)~morel[,1],data=morel) newfit2<-linearHypothesis(fitp2,hypothesis.matrix=C) print(newfit2) # # How to get other contrasts? Suppose we wanted \tau_2 = 0. ##y_{ij}=mu_i+e_{ij}=mu+tau_i+e_{ij}: mu2=mu, mu1=mu+tau1, mu3=mu+tau3 # Use R function contr.treatment # # morel$studentgroup <- C(object = morel$studentgroup, contr = contr.treatment(n = 3, base = 2)) # # Note that C() ["C" is in capitals] and is not the combine c() function. # # Also, base = 2 specifies that we want the second treatment to be the baseline fit.lm2 <- lm(cbind(aptitude, mathematics, language, generalknowledge)~studentgroup , data = morel) summary(fit.lm2) # # alternative (preferred way) defining C # # C1 <- matrix(c(1, 1, 0, 0, 0, -1, 0, -1, 1), ncol = 3, by = TRUE) newfit1 <- linearHypothesis(model = fit.lm, hypothesis.matrix = C1) C1 %*% fit.lm$coef fit.lm2$coef # compare with fit.lm1$coef (exactly the same) ##--------------------------------------------------------------------------------------------------------------------------------- ##Two-way MANOVA ## Example 2 ## peanuts <- read.table(file = "datafiles/T6-17.dat", col.names=c("location", "variety", "yield", "weight", "seedsize")) peanuts$location <- as.factor(peanuts$location) peanuts$variety <- as.factor(peanuts\$variety) fit.lm <- lm(cbind(yield, weight, seedsize) ~ location * variety, data = peanuts) fit.manova <- Manova(fit.lm) summary(fit.manova) ####Appendix programs trteff.diff<-function (data, level=0.95) { # Bonferroni Based Simultaneous CI for Treatments’ Difference for one-way MANOVA #data: Nx(p+1) numeric matrix or dataframe of data, N is the total sample size, p is number of variables, # the last column is the factors #level: confidence level of interval,default=0.95 Y<-data[,-ncol(data)] X<-data[,ncol(data)] p <- ncol(Y) g <- length(levels(X)) X <- as.numeric(X) data <- as.matrix(cbind(Y,X)) N <- length(X) meanvec <- matrix(apply(Y, 2, mean), ncol = 1) n <- matrix(rep(0, g), ncol = 1) trtmean <- matrix(rep(0, p * g), ncol = g) trt.effect <- matrix(rep(0, p * g), ncol = g) trt.cov <- matrix(rep(0, (p * p) * g), nrow = p) W <- matrix(rep(0, (p * p)), nrow = p) df2 <- 0 for (k in 1:g) { n[k] <- length(subset(X, X == k)) trtmean[, k] <- as.matrix(colMeans(subset(Y, X == k))) trt.effect[, k] <- trtmean[, k] - meanvec for (i in 1:p) { for (j in 1:p) { trt.cov[j, i + (k - 1) * p] <- cov(subset(Y[,i], X == k), subset(Y[, j], X == k)) } } W = W + (n[k] - 1) * trt.cov[, (1 + (k - 1) * p):(k *p)] } prob <- 1 - ((1 - level)/(p * g * (g - 1))) t <- qt(prob, (N - g)) d <- (g * (g - 1)/2) tau <- matrix(rep(0, (p * d)), ncol = 1) lower <- matrix(rep(0, (p * d)), ncol = 1) upper <- matrix(rep(0, (p * d)), ncol = 1) C<-matrix(0,nrow=d,g) kk<-1 for(k in 1:(g-1)) for(l in (k+1):g){ C[kk,k]<-1 C[kk,l]<--1 kk<-kk+1 } for (i in 1:p) { for (k in 1:(g - 1)) { for (l in (k + 1):g) { tau[(1 + (i - 1) * d):(i * d), ] <- C%*%(as.matrix(trt.effect[i, ])) lower[(1 + (i - 1) * d):(i * d), ] <- tau[(1 +(i - 1) * d):(i * d), ] - t * sqrt((W[i, i]/(N - g)) * ((1/n[k]) + (1/n[l]))) upper[(1 + (i - 1) * d):(i * d), ] <- tau[(1 +(i - 1) * d):(i * d), ] + t * sqrt((W[i, i]/(N - g)) * ((1/n[k]) + (1/n[l]))) } } } cat(trt.effect,"\n") cat("\n Bonferroni Based Simultaneous CI for Treatments’ Difference \n") SCI = data.frame(Trt = C, Estimate = round(tau,3), LowerCI = round(lower, 3), UpperCI = round(upper,3)) labs<-NA for(i in 1:p) for(k in 1:(g-1)) for(l in (k+1):g) labs<-c(labs,paste("tau",k,"-",l,i,sep="")) rownames( SCI)<-labs[-1] print(SCI) }