Univariate Normality

stiff <- read.table("resources/data/T4-3.DAT", header = FALSE)[,1:4]
head(stiff)

colnames(stiff) <- c("x1","x2","x3","x4")
attach(stiff)

###########
## Univariate normality
###########
par(mfrow = c(2,2))
hist(x1)
hist(x4)
qqnorm(x1);qqline(x1)  # Q-Q plot of x1
qqnorm(x4);qqline(x4)  # Q-Q plot of x4

# Shapiro Test for Univariate Normality
shapiro.test(x1)

## 
##  Shapiro-Wilk normality test
## 
## data:  x1
## W = 0.93068, p-value = 0.05118

shapiro.test(x2)

## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.91274, p-value = 0.01746

shapiro.test(x3)

## 
##  Shapiro-Wilk normality test
## 
## data:  x3
## W = 0.93258, p-value = 0.05751

shapiro.test(x4)

## 
##  Shapiro-Wilk normality test
## 
## data:  x4
## W = 0.96127, p-value = 0.3337

Bivariate Normality

#Scatter Plot Matrix
pairs(stiff)   

#Bivariate Percentage
par(mfrow = c(1,1))  
plot(x1, x4, pch = 20)
quan_chisq = ellipse(cor(stiff[,c(1,4)]), 
                     scale = c(sd(x1),sd(x4)), 
                     center = c(0,0), 
                     level = 0.5, 
                     which = c(1, 2), 
                     npoints = 100)
lines(quan_chisq[,1] + mean(x1),
      quan_chisq[,2] + mean(x4),
      col = "blue",
      lwd = 2)

# Mahalanobis distance between x1 & x4
m_dist14 <- mahalanobis(x = stiff[,c(1,4)],
                        center = colMeans(stiff[,c(1,4)]),
                        cov = cov(stiff[,c(1,4)]))

# 60% of Mahalanobis distance <= 50th quantile of chi-square
sum(m_dist14<=qchisq(.5,2))/30  

## [1] 0.6

# Calculate Mahalanobis distance
m_dist <- mahalanobis(x = stiff, 
                      center = colMeans(stiff),
                      cov = cov(stiff))

sum (m_dist <= qchisq(.05,4))

## [1] 3

par(mar=c(5,5.5,4,2)+0.1)  # by hand calculation
plot(sort(m_dist),qchisq(((1:30)-0.5)/30,4),xlab=expression(d[j]^2),
     ylab=expression(paste(chi[4]^2,"( ",frac(j-.5,30)," )")),pch=20)
abline(a=0,b=1,lwd=2)

# Directly by function mvn()
mvn(stiff, 
    mvnTest = "royston", 
    multivariatePlot="qq")  

## $multivariateNormality
##      Test        H     p value MVN
## 1 Royston 9.823858 0.009534665  NO
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    x1        0.9307    0.0512    YES   
## 2 Shapiro-Wilk    x2        0.9127    0.0175    NO    
## 3 Shapiro-Wilk    x3        0.9326    0.0575    YES   
## 4 Shapiro-Wilk    x4        0.9613    0.3337    YES   
## 
## $Descriptives
##     n     Mean  Std.Dev Median  Min  Max    25th    75th      Skew
## x1 30 1906.100 324.9866 1863.0 1325 2983 1715.25 2057.25 1.0380842
## x2 30 1749.533 318.6065 1680.0 1170 2794 1595.50 1888.75 1.1435912
## x3 30 1509.133 303.1783 1466.0 1002 2412 1295.75 1623.75 0.9800274
## x4 30 1724.967 322.8436 1674.5 1176 2581 1520.25 1880.75 0.5978431
##       Kurtosis
## x1  2.03586397
## x2  1.94986381
## x3  0.99683699
## x4 -0.04626509

Detecting Outliers

# Standardized Value
print(data.frame(scale(x1),scale(x2),scale(x3),scale(x4)))

##      scale.x1.  scale.x2.   scale.x3.   scale.x4.
## 1  -0.05261755 -0.3092634  0.17107644  0.16426945
## 2   1.52898605  0.9367877  1.90602908  1.46211166
## 3   0.65510390 -0.1554687  1.00886726  1.53954854
## 4  -0.80341770 -0.3845914 -1.31649701 -0.59461204
## 5   0.21508578  0.5224835  0.34589106  0.48950437
## 6  -0.59725537 -0.1178047 -0.23132702 -0.55434486
## 7   0.11354314 -0.2025487 -0.78545638 -0.16716042
## 8   0.60894816  0.2211714  0.68562513  0.46162709
## 9   3.31367493  3.2782337  2.97800552  2.65154223
## 10 -0.49571272 -0.4693354 -0.41273841 -0.67204892
## 11 -0.60340947 -0.4975834  0.02924572 -0.17955033
## 12  0.43047927  0.4942355  0.38877012  0.53596650
## 13 -0.20339299  0.2870835  0.28322167  0.04966286
## 14 -0.12031265 -0.2025487 -0.05321401 -0.14547810
## 15 -0.14492905 -0.3155407 -0.39624647 -0.03396898
## 16  0.14739069  1.2537931 -1.08560978 -1.37517585
## 17 -1.78807364 -1.8189625 -1.67272303 -1.70041077
## 18 -1.49883096 -1.1880903 -0.84812577 -1.29154401
## 19 -0.24031759 -0.3626207  0.30631040  0.09302751
## 20 -0.55725372 -0.4881674 -0.64692404 -0.24459731
## 21  1.13820072  1.3793398  0.12489900  1.19572877
## 22 -0.02184705 -0.4253941 -0.28739963 -0.76807066
## 23 -0.84034230 -0.7423995 -0.72278698 -0.64726912
## 24  0.47663501  0.3686888  0.45143951  0.96651559
## 25 -0.15416020 -0.8051729 -0.50509331 -0.59461204
## 26 -0.55109962 -1.0594050 -0.89430321 -0.79285046
## 27  0.80587934  0.4597102  0.63285091  0.33772807
## 28 -0.77264721 -0.2339354 -0.31378674 -0.39637361
## 29  1.29205321  1.7308706  1.83346452  1.57671825
## 30 -1.28036042 -1.1535650 -0.97346455 -1.36588342

# Generalized Squared Distance--Mahalanobis distance
rank(m_dist)

##  [1]  3 25 27 24  8 12 22 10 29  4 11  2 16  1  6 30 19 20  7  9 28 23  5
## [24] 14 21 18 13 17 26 15

cbind(1:30,m_dist)[order(m_dist),]

##              m_dist
##  [1,] 14  0.1295714
##  [2,] 12  0.4635159
##  [3,]  1  0.6000129
##  [4,] 10  0.7665400
##  [5,] 23  0.7962096
##  [6,] 15  1.0792484
##  [7,] 19  1.3632124
##  [8,]  5  1.3980776
##  [9,] 20  1.4649908
## [10,]  8  1.4876570
## [11,] 11  1.9307771
## [12,]  6  2.2191409
## [13,] 27  2.3816050
## [14,] 24  2.5385575
## [15,] 30  2.5838186
## [16,] 13  2.6959024
## [17,] 28  2.9951752
## [18,] 26  3.3979804
## [19,] 17  3.5018290
## [20,] 18  3.9900603
## [21,] 25  4.5767867
## [22,]  7  4.9883498
## [23,] 22  5.0557446
## [24,]  4  5.2076098
## [25,]  2  5.4770196
## [26,] 29  6.2837628
## [27,]  3  7.6166439
## [28,] 21  9.8980384
## [29,]  9 12.2647550
## [30,] 16 16.8474070

# By MVN Function
mvn(stiff[-c(16,9),], mvnTest = "royston", multivariatePlot="qq")  

## $multivariateNormality
##      Test        H   p value MVN
## 1 Royston 1.098338 0.6271166 YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    x1        0.9847    0.9439    YES   
## 2 Shapiro-Wilk    x2        0.9664    0.4871    YES   
## 3 Shapiro-Wilk    x3        0.9672    0.5070    YES   
## 4 Shapiro-Wilk    x4        0.9660    0.4779    YES   
## 
## $Descriptives
##     n     Mean  Std.Dev Median  Min  Max    25th    75th       Skew
## x1 28 1865.929 262.1619 1857.5 1325 2403 1711.50 2049.75 0.08994538
## x2 28 1697.964 244.8618 1663.0 1170 2301 1593.25 1847.50 0.39091767
## x3 28 1488.643 253.1536 1466.0 1002 2087 1307.25 1617.25 0.49661284
## x4 28 1710.250 277.9986 1674.5 1176 2234 1528.75 1876.25 0.25921958
##      Kurtosis
## x1 -0.5084972
## x2  0.1961808
## x3  0.0516768
## x4 -0.6484407

Transformations

# yjPower" family (Yeo and Johnson(2000)),
# another modifiation of the Box-Cox family that allows a few negative values.
trans <- powerTransform(stiff, 
                        family = "yjPower")  
summary(trans)

## yjPower Transformations to Multinormality 
##    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## x1    0.0939           1      -0.9911       1.1790
## x2   -0.2802           0      -1.5129       0.9525
## x3    0.1478           1      -0.9585       1.2542
## x4    0.7546           1      -0.5292       2.0385
## 
##  Likelihood ratio test that all transformation parameters are equal to 0
##                                  LRT df    pval
## LR test, lambda = (0 0 0 0) 1.980807  4 0.73929

logstiff = log(stiff)

# Histogram to assess normality
hist(log(x1))

# Scatterplot
plot(log(x1),log(x4),pch=20)
# Add Eclipse
quan_chisq1 <- ellipse(cor(logstiff[,c(1,4)]), 
                       scale = c(sd(log(x1)),
                                 sd(log(x4))), 
                       center=c(0,0), 
                       t = sqrt(qchisq(.5, 2)), 
                       which = c(1, 2), 
                       npoints = 100)
# Add Lines
lines(quan_chisq1[,1] + mean(log(x1)), 
      quan_chisq1[,2] + mean(log(x4)), 
      col="blue", 
      lwd=2)

# Compare to test on untransformed stiff
mvn(logstiff, mvnTest = "royston", multivariatePlot="qq")

## $multivariateNormality
##      Test        H   p value MVN
## 1 Royston 1.652091 0.5200565 YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    x1        0.9742    0.6577    YES   
## 2 Shapiro-Wilk    x2        0.9626    0.3612    YES   
## 3 Shapiro-Wilk    x3        0.9795    0.8118    YES   
## 4 Shapiro-Wilk    x4        0.9831    0.9005    YES   
## 
## $Descriptives
##     n     Mean   Std.Dev   Median      Min      Max     25th     75th
## x1 30 7.539650 0.1631929 7.529941 7.189168 8.000685 7.447309 7.629120
## x2 30 7.452334 0.1721702 7.426545 7.064759 7.935230 7.374941 7.543648
## x3 30 7.301165 0.1910199 7.290123 6.909753 7.788212 7.166816 7.392488
## x4 30 7.436553 0.1834145 7.423268 7.069874 7.855932 7.326618 7.539424
##         Skew    Kurtosis
## x1 0.3660786  0.69728520
## x2 0.4814651  0.73170914
## x3 0.4107699  0.07436358
## x4 0.1581748 -0.44895116

Sample Statistics Sweat

We want to see if perspiration in Latin American healthy women in the age group 50-65 is the same as their counterparts in the North America. Perspiration is quantified in terms of the following variables
* X1: Sweat Rate
* X2: Sodium in Sweat
* X3: Potassium in Sweat

rm(list = ls())

sweat <- read.table("resources/data/T5-1.DAT")
head(sweat)

n <- nrow(sweat)
p <- ncol(sweat)
alpha <- 0.05

# Null Hypothesis Mean Values
mu0 <- matrix(c(4,50,10),ncol=1)

# Sample Means
xbar <- as.matrix(colMeans(sweat), ncol = 1)

# Sample Covariance Matrix
S <- cov(sweat)

Hotelling Test

# Calculating Hotelling
## By hand
Tsq <- n*t(xbar-mu0)%*%solve(S)%*%(xbar-mu0)

## By Mahalanobis distance
Tsq <- n*mahalanobis(x = colMeans(sweat), 
                     center = mu0, 
                     cov = cov(sweat))

# Critical point at alpha = 0.05
cri <- (((n-1)*p)/(n-p))*qf((1-alpha),p,n-p) 

# Tsq < cri Fail to reject H0
Tsq > cri

## [1] FALSE

# Compute P value
Tsq.mod <- ((n-p)/p)/(n-1)*Tsq
p.val <- 1 - pf(Tsq.mod, p, (n-p))

test.frame <- data.frame(stat = Tsq,
                         mod.stat = Tsq.mod,
                         crit = cri,
                         reject = (Tsq > cri),
                         reject.mod = (Tsq.mod > cri),
                         p = p.val)
test.frame

# Hotellings T^2 with Function
HotellingsT2(sweat, 
             mu = mu0)

## 
##  Hotelling's one sample T2-test
## 
## data:  sweat
## T.2 = 2.9045, df1 = 3, df2 = 17, p-value = 0.06493
## alternative hypothesis: true location is not equal to c(4,50,10)

Confidence Intervals for \(\mu\)

# One at a time intervals of mu
One <- matrix(NA,3,2)
for (i in 1:3){
  One[i,1] <- xbar[i] - qt(.975,n-1)*sqrt(S[i,i]/n) #left
  One[i,2] <- xbar[i] + qt(.975,n-1)*sqrt(S[i,i]/n) #right
}

## T^2 interval mu
Tsq <- matrix(NA,3,2)
c <- sqrt(((p*(n-1))/(n-p))*qf(.95,p,n-p))
for (i in 1:3){
  Tsq[i,1] <- xbar[i] - c*sqrt(S[i,i]/n)
  Tsq[i,2] <- xbar[i] + c*sqrt(S[i,i]/n)
}

## Bonferroni interval mu when m = 3 (m is the number of CIs being estimated).
m <- 3
Bonf <- matrix(NA,3,2)
for (i in 1:m){
  Bonf[i,1] <- xbar[i] - qt((1-alpha)/(2*m),n-1)*sqrt(S[i,i]/n) #left
  Bonf[i,2] <- xbar[i] + qt((1-alpha)/(2*m),n-1)*sqrt(S[i,i]/n) #right
}

colnames(One) = colnames(Tsq) =colnames(Bonf) =c("lower","upper")
row.names(One)=row.names(Tsq)=row.names(Bonf)=c("SweatRate","Sodium","Potassium")
list(One, Tsq, Bonf)

## [[1]]
##               lower    upper
## SweatRate  3.845840  5.43416
## Sodium    38.784779 52.01522
## Potassium  9.073601 10.85640
## 
## [[2]]
##               lower     upper
## SweatRate  3.397768  5.882232
## Sodium    35.052408 55.747592
## Potassium  8.570664 11.359336
## 
## [[3]]
##               lower     upper
## SweatRate  5.030216  4.249784
## Sodium    48.650435 42.149565
## Potassium 10.402995  9.527005

# For large sample inference, only need to change the critical point or multiplier using Z (replace to t) and chi-square (replace to F).

Confidence Region/Interval for Delta

plot(D)
CR <- ellipse(x = cor(D), 
              scale = sqrt(diag(Sd)/n1), 
              centre = as.vector(dbar),
              t = sqrt(cri1))
#n1*mahalanobis(as.vector(dbar),as.vector(CR[3,]),cov(D))
lines(CR, 
      type = "l",
      xlab = expression({delta}[1]),
      ylab = expression({delta}[2]), 
      main = expression("95% confidence ellipse"))
points(0,0,pch=20)

## One at a time intervals of for delta1 and delta2
One1 <- matrix(NA,2,2)
for (i in 1:2){
  One1[i,1] <- dbar[i] - qt(.975,n1-1)*sqrt(Sd[i,i]/n1) 
  One1[i,2] <- dbar[i] + qt(.975,n1-1)*sqrt(Sd[i,i]/n1) 
}

## T^2 interval for delta1 and delta2
Tsq1 <- matrix(NA,2,2)
c1 <- sqrt(((p1*(n1-1))/(n1-p1))*qf(.95,p1,n1-p1))
for (i in 1:2){
  Tsq1[i,1] <- dbar[i] - c1*sqrt(Sd[i,i]/n1) 
  Tsq1[i,2] <- dbar[i] + c1*sqrt(Sd[i,i]/n1) 
}

## Bonferroni Interval for delta1 and delta2 when m=2 (m is the number of CIs).
m <- 2
Bonf1 <- matrix(NA,2,2)
for (i in 1:m){
  Bonf1[i,1] <- dbar[i] - qt((1-0.05/(2*m)),n1-1)*sqrt(Sd[i,i]/n1)
  Bonf1[i,2] <- dbar[i] + qt((1-0.05/(2*m)),n1-1)*sqrt(Sd[i,i]/n1)
}

cols <- c("lower","upper")
rows <- c("BOD","SS")

colnames(One1) = colnames(Tsq1) = colnames(Bonf1) = cols
row.names(One1) = row.names(Tsq1) = row.names(Bonf1) = rows
out <- list(One1, Tsq1, Bonf1)
names(out) <- c("One-at-a-time Intervals", "T^2 Intervals", "Bonferroni Simultaneous Intervals")
out

## $`One-at-a-time Intervals`
##           lower     upper
## BOD -18.8467298  0.119457
## SS   -0.4725958 27.018050
## 
## $`T^2 Intervals`
##          lower    upper
## BOD -22.453272  3.72600
## SS   -5.700119 32.24557
## 
## $`Bonferroni Simultaneous Intervals`
##          lower     upper
## BOD -20.573107  1.845835
## SS   -2.974903 29.520358

2-Sample Independent Confidence Intervals

# T^2-intervals
Clower <- t(A)%*%(xbar1 - xbar2) - sqrt(crit) * sqrt(diag((t(A)%*%Sp%*%A)*(1/n1 + 1/n2)))
Cupper <- t(A)%*%(xbar1 - xbar2) + sqrt(crit) * sqrt(diag((t(A)%*%Sp%*%A)*(1/n1 + 1/n2)))
cbind(Clower,Cupper)

##            [,1]        [,2]
## [1,] -0.2926638 -0.05776762
## [2,] -0.2365537 -0.05411666
## [3,] -0.3451377 -0.12906223

# Bonferroni Intervals
crit.b <- abs(qt(alpha/(2*p), (n1 + n2-2)))
Clower.b <- t(A)%*%(xbar1 - xbar2) - crit.b*sqrt(diag((t(A)%*%Sp%*%A)*(1/n1+1/n2)))
Cupper.b <- t(A)%*%(xbar1 - xbar2) + crit.b*sqrt(diag((t(A)%*%Sp%*%A)*(1/n1+1/n2)))
cbind(Clower.b,Cupper.b)

##            [,1]        [,2]
## [1,] -0.2734025 -0.07702893
## [2,] -0.2215940 -0.06907636
## [3,] -0.3274197 -0.14678026

# Assessing the normality and equal covariance assumptions 
par(mfrow = c(1,3))

# 3-d Scatterplot: Ellipsoidal? are the covariances the same?
sp <- scatterplot3d(X1, 
                    xlim = c(4.4,5.3), 
                    ylim = c(4.2,5), 
                    zlim = c(3.4,4.4), 
                    main = "Red = FEMALE")
sp$points3d(X2, col = "red")

#Visually compare the sample variances and covariances between the two groups
cov(X1) 

##             Length       Width      Height
## Length 0.011072004 0.008019142 0.008159648
## Width  0.008019142 0.006416726 0.006005271
## Height 0.008159648 0.006005271 0.006772758

cov(X2)

##            Length      Width     Height
## Length 0.02640563 0.02011195 0.02491758
## Width  0.02011195 0.01619045 0.01942430
## Height 0.02491758 0.01942430 0.02493980

Chi-Square Plots

# Chi-square plot for males and females separately.
distances <-function(D){
  dsq <- c()
  M <- colMeans(D) 
  S <- cov(D)
  n <- dim(D)[1]
  p <- dim(D)[2]
  for(i in 1:n){
    dsq <- c(dsq, as.matrix(D[i,]-M)%*%solve(S)%*%as.matrix(t(D[i,]-M)))
  }
  as.numeric(dsq)
}

#Chi-square plot for males
qc <- qchisq( ((1:n1)-1/2)/n1, df = p)
par(mar = c(4,4,1,1))
plot(sort(distances(X1)),
     qc,
     xlab = expression(d[j]^2),
     main="Male")
abline(a = 0, b = 1)

#Chi-square plot for females
qc <- qchisq( ((1:n2)-1/2)/n2, df = p)
par(mar = c(4,4,1,1))
plot(sort(distances(X2)), 
     qc,
     xlab = expression(d[j]^2),
     main = "Female")
abline(a = 0, b = 1)

One-Way Manova Function

test1 <- manova(cbind(D[,2],D[,4]) ~ as.character(G))
#c("Pillai", "Wilks", "Hotelling-Lawley", "Roy")
summary(test1,'Wilks') 

##                  Df    Wilks approx F num Df den Df    Pr(>F)    
## as.character(G)   2 0.038316   299.94      4    292 < 2.2e-16 ***
## Residuals       147                                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(test1)

##  Response 1 :
##                  Df Sum Sq Mean Sq F value    Pr(>F)    
## as.character(G)   2 11.345  5.6725   49.16 < 2.2e-16 ***
## Residuals       147 16.962  0.1154                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 2 :
##                  Df Sum Sq Mean Sq F value    Pr(>F)    
## as.character(G)   2 80.413  40.207  960.01 < 2.2e-16 ***
## Residuals       147  6.157   0.042                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise Comparison With Hotelling’s T²

# Testing H0:mui-muj=delta0
delta0 <- rep(0,p)
xbar1 <- colMeans(X1) 
xbar2 <- colMeans(X2) 
xbar3 <- colMeans(X3)
Sp <- W/(n-g)

# Number of pairs of mu you are comparing
m <- g*(g-1)/2  
crit <- ((n-g)*p/(n-g-p+1))*qf(1- alpha/m,p,n-g-p+1)

# Between Species 1 and 2
T2 <- (N[1]*N[2]/(N[1]+N[2]))*t(xbar1-xbar2-delta0)%*%solve(Sp)%*%(xbar1-xbar2-delta0)
# Between Species 1 and 3
T2 <- (N[1]*N[3]/(N[1]+N[3]))*t(xbar1-xbar3-delta0)%*%solve(Sp)%*%(xbar1-xbar3-delta0)
# Between Species 2 and 3
T2 <- (N[2]*N[3]/(N[2]+N[3]))*t(xbar2-xbar3-delta0)%*%solve(Sp)%*%(xbar2-xbar3-delta0)

# Bonferoni Intervals

m <- p*g*(g-1)/2
crit.b <- qt(1- alpha/(2*m), n-g)

# Comparing Species 1 and 2
xbar1-xbar2 - crit.b*sqrt((1/N[1]+1/N[2])*diag(Sp))

##         V2         V4 
##  0.4763056 -1.1894645

xbar1-xbar2 + crit.b*sqrt((1/N[1]+1/N[2])*diag(Sp))

##         V2         V4 
##  0.8396944 -0.9705355

#Comparing Species 1 and 3
xbar1-xbar3 - crit.b*sqrt((1/N[1]+1/N[3])*diag(Sp))

##         V2         V4 
##  0.2723056 -1.8894645

xbar1-xbar3 + crit.b*sqrt((1/N[1]+1/N[3])*diag(Sp))

##         V2         V4 
##  0.6356944 -1.6705355

#Comparing Species 2 and 3
xbar2-xbar3 - crit.b*sqrt((1/N[2]+1/N[3])*diag(Sp))

##         V2         V4 
## -0.3856944 -0.8094645

xbar2-xbar3 + crit.b*sqrt((1/N[2]+1/N[3])*diag(Sp))

##          V2          V4 
## -0.02230564 -0.59053548

Box’s M Test for Equality of Covariance

boxM(cbind(plastic$V3,plastic$V4,plastic$V5),plastic$V1*plastic$V2)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  cbind(plastic$V3, plastic$V4, plastic$V5)
## Chi-Sq (approx.) = 7.1269, df = 6, p-value = 0.3093

boxM <-function(data, grouping)  {
    if (!inherits(data, c("data.frame", "matrix")))
      stop("'data' must be a numeric data.frame or matrix!")
    if (length(grouping) != nrow(data))
      stop("incompatible dimensions!")
    dname <- deparse(substitute(data))
    data <- as.matrix(data)
    grouping <- as.factor(as.character(grouping))
    p <- ncol(data)
    nlev <- nlevels(grouping)
    lev <- levels(grouping)
    dfs <- tapply(grouping, grouping, length) - 1
    if (any(dfs < p)) 
      warning("there are one or more levels with less observations than variables!")
    mats <- aux <- list()
    for(i in 1:nlev) {
      mats[[i]] <- cov(data[grouping == lev[i], ])
      aux[[i]] <- mats[[i]] * dfs[i]
    }
    names(mats) <- lev
    pooled <- Reduce("+", aux) / sum(dfs)
    logdet <- log(unlist(lapply(mats, det)))
    minus2logM <- sum(dfs) * log(det(pooled)) - sum(logdet * dfs)
    sum1 <- sum(1 / dfs) 
    Co <- (((2 * p^2) + (3 * p) - 1) / (6 * (p + 1) *
                                          (nlev - 1))) * (sum1 - (1 / sum(dfs)))
    X2 <- minus2logM * (1 - Co)
    dfchi <- (choose(p, 2) + p) * (nlev - 1)
    pval <- pchisq(X2, dfchi, lower.tail = FALSE)
    out <- structure(
      list(statistic = c("Chi-Sq (approx.)" = X2),
           parameter = c(df = dfchi),
           p.value = pval,
           cov = mats, pooled = pooled, logDet = logdet,
           data.name = dname,
           method = " Box's M-test for Homogeneity of Covariance Matrices"
      ),
      class = c("htest", "boxM")
    )
    return(out)
}