多元第四次作业上机

11.某项研究确定运动或膳食补充是否会减缓妇女的骨质流失,研究人员通过光子吸收法测量了实验前和实验1年后骨骼中的矿物质含量,表6.8是对参与该项目实验前25个个体和参与该项目实验1年后24个个体骨骼中的矿物质含量数据,记录了3个骨骼上主力侧和非主力侧上矿物质含量,其中$X_1$,表示主力侧的桡骨、$X_2$表示桡骨、$X_3$表示主力侧的肱骨、$X_4$,表示肱骨、$X_5$表示主力侧的尺骨、$X_6$表示尺骨中矿物质含量.假设$X=(X_1,\cdots ,X_6{)}^{{}^{\prime }}\sim N_6(\mu ,\Sigma )$ (1)分别绘制实验前数据和实验1年后数据的矩阵散点图; (2)给定显著性水平a=0.5,检验经过实验后骨骼中的矿物质含量是否有流失? (3)构造均值差95%的同时置信区间和 Bonferroni 同时置信区间; (4)给定显著性水平a=0.5,分别对实验前和实验后的数据进行独立性检验.首先对随机向量X 和协方差矩阵习进行如下剖分: $X=(\begin{array}{c}{\mathbf{X}}^{(1)}\\ {\mathbf{X}}^{(2)}\\ {\mathbf{X}}^{(3)}\end{array}),\mathbf{\mu }=\left(\begin{array}{c}{\mathbf{\mu }}^{(1)}\\ {\mathbf{\mu }}^{(2)}\\ {\mathbf{\mu }}^{(3)}\end{array}\right),\mathbf{\Sigma }=\left(\begin{array}{ccc}{\mathbf{\Sigma }}_{11}& {\mathbf{\Sigma }}_{12}& {\mathbf{\Sigma }}_{13}\\ {\mathbf{\Sigma }}_{21}& {\mathbf{\Sigma }}_{22}& {\mathbf{\Sigma }}_{23}\\ {\mathbf{\Sigma }}_{31}& {\mathbf{\Sigma }}_{32}& {\mathbf{\Sigma }}_{33}\end{array}\right)$ 其中${\mathbf{X}}^{(1)}=(X_1,X_2{)}^{\prime },{\mathbf{X}}^{(2)}=(X_3,X_4{)}^{\prime },{\mathbf{X}}^{(3)}=(X_5,X_6{)}^{\prime }.$考虑如下的假设检验问题： $H_0:{\mathbf{\Sigma }}_{12}=0,{\mathbf{\Sigma }}_{13}=0,{\mathbf{\Sigma }}_{23}=0,H_1:{\mathbf{\Sigma }}_{12},{\mathbf{\Sigma }}_{13},{\mathbf{\Sigma }}_{23}\text{不全为 0 矩阵}.$

首先导入数据：

# 实验前数据
data_0 <- data.frame(
  X1 = c(1.103, 0.842, 0.925, 0.857, 0.795, 0.787, 0.933, 0.945, 0.921, 0.792, 0.815, 0.755, 0.880, 0.900, 0.764, 0.733, 0.932, 0.856, 0.890, 0.688, 0.940, 0.493, 0.835, 0.915),
  X2 = c(1.052, 0.859, 0.873, 0.744, 0.809, 0.779, 0.880, 0.876, 0.906, 0.825, 0.751, 0.724, 0.866, 0.838, 0.757, 0.748, 0.898, 0.786, 0.950, 0.532, 0.850, 0.616, 0.752, 0.936),
  X3 = c(2.139, 1.873, 1.887, 1.739, 1.734, 1.509, 1.695, 1.811, 1.954, 1.624, 2.204, 1.508, 1.786, 1.902, 1.743, 1.863, 2.028, 1.390, 2.187, 1.650, 2.334, 1.037, 1.509, 1.971),
  X4 = c(2.238, 1.741, 1.809, 1.547, 1.715, 1.474, 1.656, 1.759, 2.009, 1.657, 1.846, 1.458, 1.811, 1.606, 1.794, 1.869, 2.032, 1.324, 2.087, 1.378, 2.225, 1.268, 1.422, 1.869),
  X5 = c(0.873, 0.590, 0.767, 0.706, 0.549, 0.782, 0.737, 0.853, 0.823, 0.686, 0.678, 0.662, 0.810, 0.723, 0.586, 0.672, 0.836, 0.578, 0.758, 0.533, 0.757, 0.546, 0.618, 0.869),
  X6 = c(0.872, 0.744, 0.713, 0.674, 0.654, 0.571, 0.803, 0.777, 0.765, 0.668, 0.546, 0.595, 0.819, 0.677, 0.541, 0.752, 0.805, 0.610, 0.718, 0.482, 0.731, 0.615, 0.664, 0.868)
)
 
# 实验1年后数据
data_1 <- data.frame(
  X1 = c(1.027, 0.857, 0.875, 0.873, 0.811, 0.640, 0.947, 0.991, 0.977, 0.825, 0.851, 0.770, 0.912, 0.905, 0.756, 0.765, 0.932, 0.843, 0.879, 0.673, 0.949, 0.463, 0.776),
  X2 = c(1.051, 0.817, 0.880, 0.698, 0.813, 0.734, 0.865, 0.923, 0.925, 0.826, 0.765, 0.730, 0.875, 0.826, 0.727, 0.764, 0.914, 0.782, 0.906, 0.537, 0.900, 0.637, 0.743),
  X3 = c(2.268, 1.718, 1.953, 1.668, 1.643, 1.396, 1.851, 1.931, 1.933, 1.609, 2.352, 1.470, 1.846, 1.842, 1.747, 1.923, 2.190, 1.242, 2.164, 1.573, 2.130, 1.041, 1.442),
  X4 = c(2.246, 1.710, 1.756, 1.443, 1.661, 1.378, 1.686, 1.776, 2.106, 1.651, 1.980, 1.420, 1.809, 1.579, 1.860, 1.941, 1.997, 1.228, 1.999, 1.330, 2.159, 1.265, 1.411),
  X5 = c(0.869, 0.602, 0.765, 0.761, 0.551, 0.753, 0.708, 0.844, 0.869, 0.654, 0.692, 0.670, 0.823, 0.746, 0.656, 0.693, 0.883, 0.577, 0.802, 0.540, 0.804, 0.570, 0.585),
  X6 = c(0.964, 0.689, 0.738, 0.698, 0.619, 0.515, 0.787, 0.656, 0.789, 0.726, 0.526, 0.580, 0.773, 0.729, 0.506, 0.740, 0.785, 0.627, 0.769, 0.498, 0.779, 0.634, 0.640)
)

(1)利用pairs函数画出变量两两之间的散点图

pairs(data_0, main = "实验前数据的散点图")
pairs(data_1, main = "实验1年后的数据散点图")

(2)为了探究矿物质是否流失，本质在检验$H_0:{\mu }_1={\mu }_2$其中${\mu }_1$为实验前的均值向量${\mu }_2$为实验后的均值向量，这里我们调用ICSNP的HotellingsT2函数进行假设检验

library(ICSNP)
delta0 = c(0,0,0,0,0,0);
HotellingsT2(data_0,data_1,mu=delta0)
 

结果输出： Hotelling’s two sample T2-test data: data_0 and data_1 T.2 = 0.071982, df1 = 6, df2 = 40, p-value = 0.9984 alternative hypothesis: true location difference is not equal to c(0,0,0,0,0,0) 由于p值为0.9984大于显著性水平0.5，则可接受原假设，认为矿物质未流失

(3)同时置信区间 首先计算一些需要用到的数据

#准备工作
p <- ncol(data_0)  
n <- nrow(data_0)  #24位患者
m <- nrow(data_1)  #23位患者
 
#计算系数
c_square <- (n + m - 2) * p / (n + m - p - 1) * qf(0.95, p, n + m - p - 1)
 
# 计算数据均值
X_b_bar <- colMeans(data_0)  # 实验前数据均值
X_a_bar <- colMeans(data_1)  # 实验后数据均值
 
# 计算协方差矩阵
S_0 <- cov(data_0)  # 实验前协方差
S_1 <- cov(data_1)  # 实验后协方差
 
# 计算联合协方差矩阵 S_pool
S_pool <- ((n - 1) * S_0 + (m - 1) * S_1) / (n + m - 2)
 
# 计算 c = sqrt(c_square)
c <- sqrt(c_square)

对于同时置信区间，选择$a$为沿轴方向的单位向量的特殊形式

# 计算各骨骼矿物质含量的置信区间
ci_results <- data.frame(Variable = character(), Lower = numeric(), Upper = numeric())
 
for (i in 1:6) {
  # 计算置信区间的左右两端
  mu_L <- X_b_bar[i] - X_a_bar[i] - c * sqrt((1/m + 1/n) * S_pool[i, i])
  mu_U <- X_b_bar[i] - X_a_bar[i] + c * sqrt((1/m + 1/n) * S_pool[i, i])
  
  # 将结果添加到数据框中
  ci_results <- rbind(ci_results, data.frame(Variable = paste("X", i, sep = ""), Lower = mu_L, Upper = mu_U))
}
 
print(ci_results)

结果如下

$$\begin{array}{ccc}Variable& Lower& Upper\\ X1& -0.1341501& 0.1474835\\ X2& -0.1207036& 0.1339246\\ X3& -0.3446278& 0.3750735\\ X4& -0.3070870& 0.3479493\\ X5& -0.1319822& 0.1204170\\ X6& -0.1183431& 0.1359663\end{array}$$

对于Bonferroni置信区间有类似的代码

Bci_results <- data.frame(Variable = character(), Lower = numeric(), Upper = numeric())
t <- qt(0.05/(2*p),n+m-2,lower.tail =FALSE)
for (i in 1:6) {
  # 计算置信区间的左右两端
  mu_LB <- X_b_bar[i] - X_a_bar[i] - t * sqrt((1/m + 1/n) * S_pool[i, i])
  mu_UB <- X_b_bar[i] - X_a_bar[i] + t * sqrt((1/m + 1/n) * S_pool[i, i])
  
  # 将结果添加到数据框中
  Bci_results <- rbind(Bci_results, data.frame(Variable = paste("X", i, sep = ""), Lower = mu_LB, Upper = mu_UB))
}
 
# 输出结果
print(Bci_results)

$$\begin{array}{ccc}Variable& Lower& Upper\\ X1& -0.09120970& 0.10454303\\ X2& -0.08188064& 0.09510165\\ X3& -0.23489559& 0.26534124\\ X4& -0.20721419& 0.24807651\\ X5& -0.09349915& 0.08193393\\ X6& -0.07956880& 0.09719199\end{array}$$

(4)独立性检验 1.对实验前数据进行独立性检验

V <- (n - 1) * S_0
# 创建一个列表来存储 2x2 的子矩阵
block_list <- list()
 
# 分块
for (i in seq(1, 6, by = 2)) {
  for (j in seq(1, 6, by = 2)) {
    block_list[[length(block_list) + 1]] <- V[i:(i + 1), j:(j + 1)]
  }
}
 
# 计算行列式
det_block1 <- det(block_list[[1]])
det_block5 <- det(block_list[[5]])
det_block9 <- det(block_list[[9]])

A <- det(V) / (det_block1 * det_block5 * det_block9)
q <- 2  
rho <- 1 - (p^3 - 3*q^3) / (3 * (p^2 - 3 * q^2) + 3/2) / n  
f <- (p^2 - 3 * q^2) / 2  
T <- -2 * rho * (n - 1) / 2 * log(A)  
p.value <- 1 - pchisq(T, f)  

得到p值为6.61e-8远远小于0.05，故可以拒绝原假设，认为相互之间不独立 类似的有实验后的独立性检验

V <- (m - 1) * S_1
# 创建一个列表来存储 2x2 的子矩阵
block_list <- list()
 
# 分块
for (i in seq(1, 6, by = 2)) {
  for (j in seq(1, 6, by = 2)) {
    block_list[[length(block_list) + 1]] <- V[i:(i + 1), j:(j + 1)]
  }
}
 
# 计算行列式
det_block1 <- det(block_list[[1]])
det_block5 <- det(block_list[[5]])
det_block9 <- det(block_list[[9]])
 
A <- det(V) / (det_block1 * det_block5 * det_block9)
q <- 2  
rho <- 1 - (p^3 - 3*q^3) / (3 * (p^2 - 3 * q^2) + 3/2) / m  
f <- (p^2 - 3 * q^2) / 2  
T <- -2 * rho * (m - 1) / 2 * log(A)  
p.value <- 1 - pchisq(T, f)  

最后得到p值为2.61e-7，也可以拒绝原假设，认为不独立

14.自己设计一个模拟例子，编写程序，对6.5.2节介绍的球形检验问题进行模拟研究

# 加载必要的包
library(MASS)
 
# 设置参数
p <- 6
n <- 30
 
# 设置均值向量和协方差矩阵
mu <- rep(0, p) # 均值向量
sigma <- diag(p) # 协方差矩阵
 
# 生成多元正态分布数据
set.seed(123) # 设置随机种子以确保结果可重复
data <- mvrnorm(n, mu = mu, Sigma = sigma)
 
# 计算相关矩阵和统计量
V <- (n-1) * cov(data)
rho <- 1 - (2*p^2 + p + 2) / (6*p*(n-1))
T <- -2 * rho * ((n-1)/2 * log(det(V)) - (n-1)*p/2 * log(sum(diag(V))/p))
f <- p*(p+1)/2 - 1
p.value <- 1 - pchisq(T, f)
 
# 输出p值
print(p.value)

得到p值为0.89555>0.05,因此可以认为协方差矩阵与$I_p$成比例，与模拟数据情况保持一致。