tusk <- read.csv("http://www.stat.ufl.edu/~winner/data/elephant_ivory_isotope.csv") attach(tusk); names(tusk) Region <- factor(Region) X <- cbind(delta13C,delta15N,delta18O,delta2H,delta34S) X <- X[1:487,] X1 <- X[1:120,] X2 <- X[121:157,] X3 <- X[158:418,] X4 <- X[419:487,] n1 <- nrow(X1); I_n1 <- diag(n1); J_n1 <- matrix(rep(1,n1^2),n1,n1) n2 <- nrow(X2); I_n2 <- diag(n2); J_n2 <- matrix(rep(1,n2^2),n2,n2) n3 <- nrow(X3); I_n3 <- diag(n3); J_n3 <- matrix(rep(1,n3^2),n3,n3) n4 <- nrow(X4); I_n4 <- diag(n4); J_n4 <- matrix(rep(1,n4^2),n4,n4) n <- nrow(X); I_n <- diag(n); J_n <- matrix(rep(1,n^2),n,n) (xbar1 <- (1/n1) * (t(X1) %*% rep(1,n1))) (xbar2 <- (1/n2) * (t(X2) %*% rep(1,n2))) (xbar3 <- (1/n3) * (t(X3) %*% rep(1,n3))) (xbar4 <- (1/n4) * (t(X4) %*% rep(1,n4))) (xbar <- (1/n) * (t(X) %*% rep(1,n))) SSCP1 <- t(X1) %*% (I_n1 - (1/n1)*J_n1) %*% X1 SSCP2 <- t(X2) %*% (I_n2 - (1/n2)*J_n2) %*% X2 SSCP3 <- t(X3) %*% (I_n3 - (1/n3)*J_n3) %*% X3 SSCP4 <- t(X4) %*% (I_n4 - (1/n4)*J_n4) %*% X4 SSCP <- t(X) %*% (I_n - (1/n)*J_n) %*% X B <- n1*((xbar1-xbar)%*%t(xbar1-xbar)) + n2*((xbar2-xbar)%*%t(xbar2-xbar)) + n3*((xbar3-xbar)%*%t(xbar3-xbar)) + n4*((xbar4-xbar)%*%t(xbar4-xbar)) W <- SSCP1 + SSCP2 + SSCP3 + SSCP4 B W SSCP Lambda.star <- det(W) / det(B+W) p <- ncol(X); g <- 4 (Bartlett.TS <- -(n-1-(p+g)/2)*log(Lambda.star)) (Bartlett.CV <- qchisq(.95,p*(g-1))) (Bartlett.PV <- 1-pchisq(Bartlett.TS, p*(g-1))) ### Approximate F-test df.trt <- g-1 df.err <- n-g f <- df.err - (p-df.trt+1)/2 d2.1 <- (df.trt*p)^2 - 4 d2.2 <- p^2 + df.trt^2 - 5 if (d2.2 <= 0) d2.2 <- d2.1 d2 <- d2.1 / d2.2 d <- sqrt(d2) lambda <- (p*df.trt-2)/4 (F.df1 <- p*df.trt); (F.df2 <- f*d-2*lambda) (F.TS <- ((1-Lambda.star^(1/d))/Lambda.star^(1/d))*(F.df2/F.df1)) (F.CV <- qf(.95,F.df1,F.df2)) (F.PV <- 1-pf(F.TS,F.df1,F.df2)) ### Bonferroni CI's Sp.vec <- diag(W/(n-g)) num.comp <- p*g*(g-1)/2 t.bon <- qt(1-.05/(2*num.comp),n-g) x1.x2.lo <- (xbar1-xbar2) - t.bon * sqrt(Sp.vec*((1/n1)+(1/n2))) x1.x2.hi <- (xbar1-xbar2) + t.bon * sqrt(Sp.vec*((1/n1)+(1/n2))) x1.x2.out <- cbind(xbar1-xbar2,x1.x2.lo,x1.x2.hi) colnames(x1.x2.out) <- c("Mean Diff", "Lower Bound", "Upper Bound") round((x1.x2.out),3) x1.x3.lo <- (xbar1-xbar3) - t.bon * sqrt(Sp.vec*((1/n1)+(1/n3))) x1.x3.hi <- (xbar1-xbar3) + t.bon * sqrt(Sp.vec*((1/n1)+(1/n3))) x1.x3.out <- cbind(xbar1-xbar3,x1.x3.lo,x1.x3.hi) colnames(x1.x3.out) <- c("Mean Diff", "Lower Bound", "Upper Bound") round((x1.x3.out),3) x1.x4.lo <- (xbar1-xbar4) - t.bon * sqrt(Sp.vec*((1/n1)+(1/n4))) x1.x4.hi <- (xbar1-xbar4) + t.bon * sqrt(Sp.vec*((1/n1)+(1/n4))) x1.x4.out <- cbind(xbar1-xbar4,x1.x4.lo,x1.x4.hi) colnames(x1.x4.out) <- c("Mean Diff", "Lower Bound", "Upper Bound") round((x1.x4.out),3) x2.x3.lo <- (xbar2-xbar3) - t.bon * sqrt(Sp.vec*((1/n2)+(1/n3))) x2.x3.hi <- (xbar2-xbar3) + t.bon * sqrt(Sp.vec*((1/n2)+(1/n3))) x2.x3.out <- cbind(xbar2-xbar3,x2.x3.lo,x2.x3.hi) colnames(x2.x3.out) <- c("Mean Diff", "Lower Bound", "Upper Bound") round((x2.x3.out),3) x2.x4.lo <- (xbar2-xbar4) - t.bon * sqrt(Sp.vec*((1/n2)+(1/n4))) x2.x4.hi <- (xbar2-xbar4) + t.bon * sqrt(Sp.vec*((1/n2)+(1/n4))) x2.x4.out <- cbind(xbar2-xbar4,x2.x4.lo,x2.x4.hi) colnames(x2.x4.out) <- c("Mean Diff", "Lower Bound", "Upper Bound") round((x2.x4.out),3) x3.x4.lo <- (xbar3-xbar4) - t.bon * sqrt(Sp.vec*((1/n3)+(1/n4))) x3.x4.hi <- (xbar3-xbar4) + t.bon * sqrt(Sp.vec*((1/n3)+(1/n4))) x3.x4.out <- cbind(xbar3-xbar4,x3.x4.lo,x3.x4.hi) colnames(x3.x4.out) <- c("Mean Diff", "Lower Bound", "Upper Bound") round((x3.x4.out),3) ### MANOVA function in R tusk.fit <- manova(cbind(delta13C,delta15N,delta18O,delta2H,delta34S) ~ Region, subset=(Region != "Asia")) summary(tusk.fit, test="Wilks")