tusk <- read.csv("http://www.stat.ufl.edu/~winner/data/elephant_ivory_isotope.csv") attach(tusk); names(tusk) X <- cbind(delta13C,delta15N,delta18O,delta2H,delta34S) ## X <- X+100 ## X <- (X/100)^3 Region.Af <- Region[1:487] Region.Af.id <- ifelse(Region.Af == "Central Africa", 1, ifelse(Region.Af == "East Africa", 2, ifelse(Region.Af == "Southern Africa", 3, ifelse(Region.Af == "West Africa", 4, 0)))) 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 g <- 4; p <- 5 S1 <- SSCP1/(n1-1); S2 <- SSCP2/(n2-1); S3 <- SSCP3/(n3-1); S4 <- SSCP4/(n4-1) Sp <- W / (N-g) M <- (N-g)*log(det(Sp)) - ((n1-1)*log(det(S1)) + (n2-1)*log(det(S2)) + (n3-1)*log(det(S3)) + (n4-1)*log(det(S4))) u1 <- (1/(n1-1)) + (1/(n2-1)) + (1/(n3-1)) + (1/(n4-1)) - (1/(N-g)) u2 <- (2*p^2+3*p-1) / (6*(p+1)*(g-1)) u <- u1*u2 Box.C <- (1-u) * M C.df <- p*(p+1)*(g-1)/2 Box.C.CV <- qchisq(.95,C.df) Box.C.PV <- 1-pchisq(Box.C,C.df) Box.out <- cbind(Box.C, C.df, Box.C.CV, Box.C.PV) colnames <- c("Box's M", "df", "X2(.05)", "P-value") round(Box.out, 3) X.reg <- data.frame(X1=delta13C[1:487],X2=delta15N[1:487],X3=delta18O[1:487], X4=delta2H[1:487], X5=delta34S[1:487] , Region.Af.id) # install.packages("heplots") library(heplots) boxM(cbind(X1,X2,X3,X4,X5) ~ factor(Region.Af.id), data=X.reg) ## Normal Probability Plots by Variable/Region par(mfrow=c(2,2)) qqnorm(X1[,1]); qqline(X1[,1]) qqnorm(X2[,1]); qqline(X2[,1]) qqnorm(X3[,1]); qqline(X3[,1]) qqnorm(X4[,1]); qqline(X4[,1]) qqnorm(X1[,2]); qqline(X1[,2]) qqnorm(X2[,2]); qqline(X2[,2]) qqnorm(X3[,2]); qqline(X3[,2]) qqnorm(X4[,2]); qqline(X4[,2]) qqnorm(X1[,3]); qqline(X1[,3]) qqnorm(X2[,3]); qqline(X2[,3]) qqnorm(X3[,3]); qqline(X3[,3]) qqnorm(X4[,3]); qqline(X4[,3]) qqnorm(X1[,4]); qqline(X1[,4]) qqnorm(X2[,4]); qqline(X2[,4]) qqnorm(X3[,4]); qqline(X3[,4]) qqnorm(X4[,4]); qqline(X4[,4])