对例题3.2.1人的出汗多少与人体内钠和钾的含量进行分析 X1表示成年女性的出汗量 X2为钠的含量 X3为钾的含量 接下来检验\(H0:\mu=\mu_0=(4,50,10)'\)，\(H1:\mu \neq \mu_0=(4,50,10)'\) 首先导入数据
data <- data.frame(
X1 = c(3.7, 4.7, 3.8, 3.2, 3.1, 4.6, 2.4, 7.2, 6.7, 5.4,
3.9, 4.5, 3.5, 4.5, 1.5, 8.5, 4.5, 6.5, 4.1, 5.5),
X2 = c(48.5, 65.1, 47.2, 53.2, 55.5, 36.1, 24.8, 33.1, 47.4, 54.1,
36.9, 58.8, 27.8, 40.2, 13.5, 56.4, 71.6, 52.8, 44.1, 40.9),
X3 = c(9.3, 8.0, 10.9, 12.0, 9.7, 7.9, 14.0, 7.6, 8.5, 11.3,
12.7, 12.3, 9.8, 8.4, 10.1, 7.1, 8.2, 10.9, 11.2, 9.4)
)
data#（和课本上有一点不一样，不改了）## X1 X2 X3
## 1 3.7 48.5 9.3
## 2 4.7 65.1 8.0
## 3 3.8 47.2 10.9
## 4 3.2 53.2 12.0
## 5 3.1 55.5 9.7
## 6 4.6 36.1 7.9
## 7 2.4 24.8 14.0
## 8 7.2 33.1 7.6
## 9 6.7 47.4 8.5
## 10 5.4 54.1 11.3
## 11 3.9 36.9 12.7
## 12 4.5 58.8 12.3
## 13 3.5 27.8 9.8
## 14 4.5 40.2 8.4
## 15 1.5 13.5 10.1
## 16 8.5 56.4 7.1
## 17 4.5 71.6 8.2
## 18 6.5 52.8 10.9
## 19 4.1 44.1 11.2
## 20 5.5 40.9 9.4记随机向量\(X=(X1,X2,X3)'\),假定\(X \sim N_3(\mu,\Sigma)\) 取检验统计量为\(F = \frac{n - p}{(n - 1)p} T^2 \quad (p = 3, n = 20)\)
Xmeans <-colMeans(data)
Xmeans## X1 X2 X3
## 4.590 45.400 9.965n <-nrow(data)
p <-ncol(data)
print(n)## [1] 20A<- (n-1)*cov(data) # 样本离差阵 = (n-1) * 协方差矩阵
A## X1 X2 X3
## X1 53.538 170.49 -32.4070
## X2 170.490 3795.98 -107.1600
## X3 -32.407 -107.16 68.9255#计算样本离差阵的逆
A_inv <-solve(A)
A_inv## X1 X2 X3
## X1 0.0291320829 -9.640614e-04 1.219831e-02
## X2 -0.0009640614 3.074329e-04 2.469584e-05
## X3 0.0121983096 2.469584e-05 2.028215e-02#计算hotellling统计量
u0 <-c(4,50,10)
D2 <- (n-1)*t(Xmeans-u0)%*%A_inv%*%(Xmeans-u0)
T2 <- n*D2
D2## [,1]
## [1,] 0.4067537T2## [,1]
## [1,] 8.135075得到上式D2和T2的值分别为0.4067537和8.135075，进一步地可以计算出F的值
F_stat <- (n-p)/((n-1)*p)*T2
F_stat## [,1]
## [1,] 2.42625得到F的值为2.42625之后，进一步对其进行假设检验，首先计算p值，在该情况下可知\(F \sim F(3,17)\)
#计算p值
p_value <-pf(F_stat,p,n-p,lower.tail = FALSE)
p_value## [,1]
## [1,] 0.1010612得到p值为\(0.1010612>0.05=\alpha\)，故接受H0
alpha <-0.05
F_cirt <-qf(1-alpha,p,n-p)
F_cirt## [1] 3.196777也可以用F分布的临界值表进行判断,得到\(F_{3,17}(0.05)=3.196777>F=2.42625\),所以接受H0.
最后对接受H0进行解释：认为意味着测量的这组数据在统计上没有显著偏离假设的均值，即这些女性的出汗量和体内钠、钾的含量与假设的值没有显著差异。可以认为人的出汗的多少与人体钾和钠的含量有一定的关系。
例题3.3.2 为了研究某种疾病，对一批人同时测量了 4 个指标：\(\beta\) 脂蛋白 (\(X_1\))，甘油三酯 (\(X_2\))，\(\alpha\) 脂蛋白 (\(X_3\))，前 \(\beta\) 脂蛋白 (\(X_4\))。按不同年龄、不同性别分为三组（20 至 35 岁的女性、20 至 25 岁的男性和 35 至 50 岁的男性），数据见表 3.3。试问这三个组的 4 项指标间有无显著性差异 (\(\alpha = 0.01\))？ 首先导入数据
# 创建数据框
body_data <- data.frame(
X1 = c(260, 200, 240, 170, 270, 205, 190, 200, 250, 200,
225, 210, 170, 270, 190, 280, 310, 270, 250, 260,
310, 310, 190, 225, 170, 210, 280, 210, 280, 200,
200, 280, 190, 295, 270, 280, 240, 280, 370, 280,
320, 260, 360, 295, 270, 380, 240, 260, 260, 295,
240, 310, 330, 345, 250, 260, 225, 345, 360, 250),
X2 = c(75, 72, 87, 65, 110, 130, 69, 46, 117, 107,
130, 125, 64, 76, 60, 81, 119, 57, 67, 135,
122, 60, 40, 65, 65, 82, 67, 38, 65, 76,
76, 94, 60, 55, 125, 120, 62, 69, 70, 40,
64, 59, 88, 100, 65, 114, 55, 55, 110, 73,
114, 103, 112, 127, 62, 59, 100, 120, 107, 117),
X3 = c(40, 34, 45, 39, 39, 34, 27, 45, 21, 28,
36, 26, 31, 33, 34, 20, 25, 31, 31, 39,
30, 35, 27, 34, 37, 31, 37, 36, 30, 40,
39, 26, 33, 30, 24, 32, 32, 29, 30, 37,
39, 37, 28, 36, 32, 36, 42, 34, 29, 33,
38, 32, 21, 24, 22, 21, 34, 36, 25, 36),
X4 = c(18, 17, 18, 17, 24, 23, 15, 15, 20, 20,
11, 17, 14, 13, 16, 18, 15, 8, 14, 29,
21, 18, 15, 16, 16, 17, 18, 17, 23, 17,
20, 11, 17, 16, 21, 18, 20, 20, 20, 17,
17, 11, 26, 12, 21, 21, 10, 20, 20, 21,
18, 18, 11, 20, 16, 19, 30, 18, 23, 16),
Group = c(rep(1, 20), rep(2, 20), rep(3, 20))
)
# 查看数据结构
str(body_data)## 'data.frame': 60 obs. of 5 variables:
## $ X1 : num 260 200 240 170 270 205 190 200 250 200 ...
## $ X2 : num 75 72 87 65 110 130 69 46 117 107 ...
## $ X3 : num 40 34 45 39 39 34 27 45 21 28 ...
## $ X4 : num 18 17 18 17 24 23 15 15 20 20 ...
## $ Group: num 1 1 1 1 1 1 1 1 1 1 ...# 查看前几行数据
head(body_data)## X1 X2 X3 X4 Group
## 1 260 75 40 18 1
## 2 200 72 34 17 1
## 3 240 87 45 18 1
## 4 170 65 39 17 1
## 5 270 110 39 24 1
## 6 205 130 34 23 1# 按组别查看数据
by(body_data, body_data$Group, head)## body_data$Group: 1
## X1 X2 X3 X4 Group
## 1 260 75 40 18 1
## 2 200 72 34 17 1
## 3 240 87 45 18 1
## 4 170 65 39 17 1
## 5 270 110 39 24 1
## 6 205 130 34 23 1
## ------------------------------------------------------------
## body_data$Group: 2
## X1 X2 X3 X4 Group
## 21 310 122 30 21 2
## 22 310 60 35 18 2
## 23 190 40 27 15 2
## 24 225 65 34 16 2
## 25 170 65 37 16 2
## 26 210 82 31 17 2
## ------------------------------------------------------------
## body_data$Group: 3
## X1 X2 X3 X4 Group
## 41 320 64 39 17 3
## 42 260 59 37 11 3
## 43 360 88 28 26 3
## 44 295 100 36 12 3
## 45 270 65 32 21 3
## 46 380 114 36 21 3为了比较三个组的四项指标间是否有差异问题，可以转化为对多总体均值向量是否相等的检验问题。设第 \(i\) 组为 4 元总体\((N_4(\mu^{(i)}, \Sigma) (i=1,2,3)\) 来自 3 个总体的样本容量 \(n_1=n_2=n_3=20\)。
检验:\(H_0: \mu^{(1)} = \mu^{(2)} = \mu^{(3)}\)
\(H_1: \mu^{(1)}, \mu^{(2)}, \mu^{(3)}\)至少有一对不等
means <-colMeans(body_data)
means## X1 X2 X3 X4 Group
## 259.08333 84.11667 32.36667 17.80000 2.00000接下来计算每一个组的样本均值
means1 <-colMeans(subset(body_data,Group==1))
means2 <-colMeans(subset(body_data,Group==2))
means3 <-colMeans(subset(body_data,Group==3))
means1## X1 X2 X3 X4 Group
## 231.0 89.6 32.9 17.1 1.0means2## X1 X2 X3 X4 Group
## 253.50 72.55 32.45 17.90 2.00means3## X1 X2 X3 X4 Group
## 292.75 90.20 31.75 18.40 3.00# 提取 Group 为 1 的子集
subset_group1 <- subset(body_data, Group == 1)
# 计算样本离差阵（协方差矩阵）
cov_matrix_group1 <- 19*cov(subset_group1[, 1:4]) # 假设 X1, X2, X3, X4 是数值列
# 打印样本离差阵
print(cov_matrix_group1)## X1 X2 X3 X4
## X1 30530 6298.0 -1078.0 198.0
## X2 6298 15736.8 -796.8 1387.8
## X3 -1078 -796.8 955.8 90.2
## X4 198 1387.8 90.2 413.8# 提取 Group 为 2的子集
subset_group2 <- subset(body_data, Group == 2)
# 计算样本离差阵（协方差矩阵）
cov_matrix_group2 <- 19*cov(subset_group2[, 1:4]) # 假设 X1, X2, X3, X4 是数值列
# 打印样本离差阵
print(cov_matrix_group2)## X1 X2 X3 X4
## X1 51705.0 7021.50 -1571.50 827.0
## X2 7021.5 12288.95 -807.95 321.1
## X3 -1571.5 -807.95 364.95 -5.1
## X4 827.0 321.10 -5.10 133.8# 提取 Group 为 3的子集
subset_group3 <- subset(body_data, Group == 3)
# 计算样本离差阵（协方差矩阵）
cov_matrix_group3 <- 19*cov(subset_group3[, 1:4]) # 假设 X1, X2, X3, X4 是数值列
# 打印样本离差阵
print(cov_matrix_group3)## X1 X2 X3 X4
## X1 43173.75 9959.0 -1301.25 723.0
## X2 9959.00 12441.2 -333.00 457.4
## X3 -1301.25 -333.0 761.75 -112.0
## X4 723.00 457.4 -112.00 476.8AA=cov_matrix_group1+cov_matrix_group2+cov_matrix_group3
AA## X1 X2 X3 X4
## X1 125408.75 23278.50 -3950.75 1748.0
## X2 23278.50 40466.95 -1937.75 2166.3
## X3 -3950.75 -1937.75 2082.50 -26.9
## X4 1748.00 2166.30 -26.90 1024.4继续计算总离差阵
# 将均值向量扩展为与数据相同维度的矩阵
means_matrix <- matrix(means[1:4],nrow=60,ncol=4,byrow=TRUE)
# 😦这里用as.matrix将数据框转化为矩阵
data_diff <- as.matrix(body_data[, 1:4] - means_matrix)
# 计算总离差阵
T <- t(data_diff) %*% data_diff
T## X1 X2 X3 X4
## X1 164474.583 25586.417 -4674.833 2534.0
## X2 25586.417 44484.183 -1973.567 2139.4
## X3 -4674.833 -1973.567 2095.933 -41.6
## X4 2534.000 2139.400 -41.600 1041.6因似然比统计量 \(\Lambda \sim \Lambda(p, n-k, k-1)\)，在比例中 \(k-1=2\)，可以利用 \(\Lambda\) 统计量与F 统计量的关系，取检验统计量为F统计量： $ F = \(\frac{(n-k)-p+1}{p} \frac{1-\sqrt{\Lambda}}{\sqrt{\Lambda}}\) (k=3, p=4, n=60), $
lambda<- det(AA)/det(T)
lambda## [1] 0.6621229f<- (54/4)*(1-sqrt(lambda))/sqrt(lambda)
f## [1] 3.090691对于给定的\(\alpha=0.01\),计算p值为
p_value <- 1 - pf(f,8,108)#自由度为8，108
p_value## [1] 0.003538437由于\(p=0.003538437<0.01=\alpha\),所以否定H0，表明三个组的指标之前有显著的差异。
进一步为了分析三个组指标间的差异是由哪几项引起的，也可以对四项指标逐项用一元方差分析方法进行检验： 首先对第一项指标X1进行检验，取\(f_1 =\frac{(T_{11} - AA_{11}) / (k - 1)}{AA_{11} / (n - k)}\)
f1=((T[1]-AA[1])/(3-1))/(AA[1]/(60-3))
f1## [1] 8.877979再由检验统计量\(F1 \sim F(2,57)\)得到p1值为
p1_value <- 1 - pf(f1,2,57)#自由度为2，57
p1_value## [1] 0.0004400761由于P1显著小于0.01，所以可以认为第一项指标X1在三个组之间存在显著差异。同理，也可以对X2，X3进行分析。