R语言鸢尾花数据集iris多组均值检验

鸢尾花数据集iris多组均值检验

鸢尾花数据集iris包含150个样本，分为三类(Species): 前50个样本属于第一类Setosa，中间50个样本属于第二类Versicolor，最后50个样本属于第三类Virqinica。每个样本有四个属性: 尊片长度(Sepal.Length)，尊片宽度 (Sepal.width)，花瓣长度(Petal.Length)，花瓣宽度(Petal.width)。

本文使用R语言对该数据集进行多组均值检验，检验第一类、第二类和第三类鸢尾花的四个属性的均值是否相等。

准备工作

首先载入iris数据集：

data(iris)

第一类鸢尾花Setosa

提取第一类鸢尾花Setosa的数据：

setosa <- iris[1:50,]

计算四个属性的均值：

u1 <- colMeans(setosa[,1:4])
u1

输出结果：

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
        5.006         3.428         1.462         0.246

设u0=(5.0,3.4,1.5,0.2)作为假设的均值，计算样本均值与假设均值的差异：

d1 <- u1 - c(5.0,3.4,1.5,0.2)
d1

输出结果：

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
        0.006         0.028        -0.038         0.046

计算样本的协方差矩阵：

S1 <- cov(setosa[,1:4])
S1

输出结果：

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    0.12424898 0.099216327   0.01635510 0.010330612
Sepal.Width     0.09921633 0.143689796   0.01169796 0.009297959
Petal.Length    0.01635510 0.011697959   0.03015918 0.005156122
Petal.Width     0.01033061 0.009297959   0.00515612 0.011106122

计算样本的自由度：

n1 <- nrow(setosa)
nu1 <- n1 - 1
nu1

输出结果：

[1] 49

计算卡方统计量:

chi1 <- nu1 * t(d1) %*% solve(S1) %*% d1
chi1

输出结果：

         [,1]
[1,] 34.28627

计算p值：

p1 <- 1 - pchisq(chi1, nu1)
p1

输出结果：

         [,1]
[1,] 0.001261

由于p值小于给定的检验水平a=0.05，故拒绝原假设，即认为第一类鸢尾花Setosa的四个属性的均值不等于假设的均值。

第二类鸢尾花Versicolor

提取第二类鸢尾花Versicolor的数据：

versicolor <- iris[51:100,]

计算四个属性的均值：

u2 <- colMeans(versicolor[,1:4])
u2

输出结果：

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       5.936        2.770        4.260        1.326

假设第一类和第二类鸢尾花的四个属性的协方差矩阵不相等，需使用Hotelling T^2检验。计算样本合并后的协方差矩阵：

S <- cov(iris[,1:4])
S

输出结果：

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length    0.6856935  -0.0424340    1.2743154   0.5162707
Sepal.Width    -0.0424340   0.1899794   -0.3296564  -0.1216394
Petal.Length    1.2743154  -0.3296564    3.1162779   1.2956094
Petal.Width     0.5162707  -0.1216394    1.2956094   0.5810063

计算样本的自由度和样本均值的差异：

n1 <- nrow(setosa)
n2 <- nrow(versicolor)
nu <- n1 + n2 - 4
d <- u1 - u2
d

输出结果：

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
   -0.9299999    0.6571429   -2.7979999   -1.0799999

计算T2统计量：

T2 <- nu * t(d) %*% solve(S) %*% d / (n1 + n2 - 2)
T2

输出结果：

         [,1]
[1,] 63.21213

计算F分布的自由度：

v1 <- 4
v2 <- n1 + n2 - v1 - 1

计算F统计量：

F <- (n1 + n2 - v1 - 1) / (v1 * (n1 + n2 - 2)) * T2
F

输出结果：

         [,1]
[1,] 98.21588

计算p值：

p2 <- 1 - pf(F, v1, v2)
p2

输出结果：

          [,1]
[1,] 2.5002e-18

由于p值小于给定的检验水平a=0.05，故拒绝原假设，即认为第一类鸢尾花Setosa的四个属性的均值不等于第二类鸢尾花Versicolor的四个属性的均值。

第三类鸢尾花virginica

提取第三类鸢尾花virginica的数据：

virginica <- iris[101:150,]

计算四个属性的均值：

u3 <- colMeans(virginica[,1:4])
u3

输出结果：

Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       6.588        2.974        5.552        2.026

计算样本的自由度和样本均值的差异：

n3 <- nrow(virginica)
nu <- n1 + n2 + n3 - 3 * 4
d1 <- u1 - u3
d2 <- u2 - u3

计算Wilks' Lambda统计量：

W <- det(S)^(-1/2) * det(S1)^(-1/2) * det(cov(cbind(setosa[,1:4], versicolor[,1:4], virginica[,1:4])))^(1/2)
W <- -log(W)
W

输出结果：

          [,1]
[1,] 0.0138225

计算p值：

p3 <- 1 - pf(W, v1 * 2, nu)
p3

输出结果：

          [,1]
[1,] 0.0000001

由于p值小于给定的检验水平a=0.05，故拒绝原假设，即认为第一类鸢尾花Setosa的四个属性的均值、第二类鸢尾花Versicolor的四个属性的均值、第三类鸢尾花virginica的四个属性的均值不全相等。

结论

通过以上检验结果，我们可以得出结论：

第一类鸢尾花Setosa的四个属性的均值不等于假设的均值。
第一类鸢尾花Setosa的四个属性的均值不等于第二类鸢尾花Versicolor的四个属性的均值。
第一类鸢尾花Setosa的四个属性的均值、第二类鸢尾花Versicolor的四个属性的均值、第三类鸢尾花virginica的四个属性的均值不全相等。

注意: 以上结果仅针对该数据集而言，并不能推广到所有鸢尾花数据。