数据集：两种形式的铁离子在不同剂量下在动物体内的存留量

本文主要讨论进行可以进行方差分析的先决条件以及等重复双因素方差分析。
链接：https://pan.baidu.com/s/1MUK4fSd_YGo1TD-20p5MXg
提取码：1234

1. 数据分析要求

为研究两种形式的铁离子（ Fe3+Fe^{3+}Fe3+和Fe2+Fe^{2+}Fe2+ ）在不同剂量下在动物体内的存留量是否有显著不同，进行了如下实验：将108只小白鼠随机地分为6组，每组均为18只，其中三组分别给以三种不同剂量（高剂量、中剂量、低剂量）的三价铁；另三组给以相应计量的二价铁。经过一段时间后，测量各小白鼠体内两种铁离子的残留量关于最初服用剂量的百分比，其数据如表2（略）所示：
（1）求出各组合水平上的观测值的样本均值和标准差。各水平组合上的标准差（从而样本方差）差异是否明显？你认为假定误差的等方差性是否合理？
（2）对观测数据作自然对数变换，再进行（1）中分析。此时，各组合水平上的标准差是否趋于一致？
（3）对变换后的数据进行方差分析，建立方差分析表。在显著水平α=0.05下，因素的交互效应是否显著？各因素的影响是否显著？
（4）根据（3）中分析，分别求各因素在其不同水平上的均值的置信度为95%置信区间以及两两均值之差的置信度不小于95%的Bonferroni同时置信区间，并解释其结果。

2. 数据预处理

# 导包
import pandas as pd
import numpy as np
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
from scipy import stats# 导入数据
df = pd.read_excel("3.2.xlsx")
print(df.head())

结果：

    A   B     C
0  a1  b1  0.71
1  a1  b1  1.66
2  a1  b1  2.01
3  a1  b1  2.16
4  a1  b1  2.42

3. 求均值与标准差并作分析

# 求均值与标准差
mean = []
std = []
j = 0
for i in range(len(df['C'])):mean.append(df['C'][j:j+18].mean())std.append(df['C'][j:j+18].mean())j = j + 18# 对数据进行对数变换
df["C"] = df["C"].apply(np.log)# 求解均值和方差
mean = []
std = []
j = 0
for i in range(len(df['C'])):mean.append(df['C'][j:j+18].mean())std.append(df['C'][j:j+18].mean())j = j + 18

原始数据：

由上表可知，(a1,b1)(a1,b1)(a1,b1)和(a1,b3)(a1,b3)(a1,b3)组合的标准差分别为2.030875、2.8067751与其他组合的标准差较大，所以我认为假定误差的等方差行不太合理。故不能直接进行方差分析。
作对数变换之后：

由以上结果可以看出，各组合水平上的标准差趋于一致，各组之间的标准差差异比较小。说明各组合的离散程度比较接近。故可以利用变换之后的数据在进行方差分析。

4. 等重复双因素方差分析

设A,BA,BA,B水平下各小白鼠体内两种铁离子的残留量关于最初服用剂量的百分比，yijky_{ijk}yijk，有：
yijk=μ+αi+βj+γij+εijk,i=1,2;j=1,2,3;k=1,2,...,18(1)y_{ijk}=\mu+\alpha_i+\beta_j+\gamma_{ij}+\varepsilon_{ijk},i=1,2;j=1,2,3;k=1,2,...,18 \tag 1 yijk=μ+αi+βj+γij+εijk,i=1,2;j=1,2,3;k=1,2,...,18(1)
检验假设为：
HA0:α1=α2=0↔HA1:αi≠0,∃iH_{A0}:\alpha_1=\alpha_2=0\leftrightarrow H_{A1}:\alpha_i\neq0,\exists i HA0:α1=α2=0↔HA1:αi=0,∃i
HB0:β1=β2=β3=0↔HB1:βj≠0,∃jH_{B0}:\beta_1=\beta_2=\beta_3=0\leftrightarrow H_{B1}:\beta_j \neq 0,\exists j HB0:β1=β2=β3=0↔HB1:βj=0,∃j
HAB0:γij=0,i=1,2;j=1,2,3↔HAB1:γij≠0,∃(i,j)H_{AB0}:\gamma_{ij}=0,i=1,2;j=1,2,3\leftrightarrow H_{AB1}:\gamma_{ij}\neq0,\exists\left(i,j\right) HAB0:γij=0,i=1,2;j=1,2,3↔HAB1:γij=0,∃(i,j)

# 等重复双因素方差分析
model = ols('C ~ A + B + A:B',df).fit()
anovat = anova_lm(model)
print(anovat)
# print(anovat['sum_sq'].sum())

结果：

在显著水平α=0.05\alpha=0.05α=0.05下，由于pA<αp_A<\alphapA<α，说明不同铁离子类型对百分比影响显著；
由于pB<αp_B<\alphapB<α，说明不同的剂量对百分比影响显著；
由于pA:B>αp_{A:B}>\alphapA:B>α，说明交互因素对百分比没有显著影响。

5. 各因素置信区间及差异置信区间

对A有：
各水平均值的置信区间：
(μAi−tα2(n−r)MSEnAi,μAi+tα2(n−r)MSEnAi),i=Fe3+,Fe2+(2)\left(\mu_{Ai}-t_\frac{\alpha}{2}\left(n-r\right)\sqrt{\frac{MS_E}{n_{Ai}}},\mu_{Ai}+t_\frac{\alpha}{2}\left(n-r\right)\sqrt{\frac{MS_E}{n_{Ai}}}\right),i=Fe^{3+},Fe^{2+} \tag 2 (μAi−t2α(n−r)nAiMSE,μAi+t2α(n−r)nAiMSE),i=Fe3+,Fe2+(2)
经查表之后有：
tscore=1.9825972t_{score}=1.9825972 tscore=1.9825972
各均值的置信区间为：
μFe3+:(1.62,1.94)μFe2+:(1.90,2.22)(3)\mu_{Fe^{3+}}:\left(1.62,1.94\right) \\ \mu_{Fe^{2+}}:\left(1.90,2.22\right) \tag 3 μFe3+:(1.62,1.94)μFe2+:(1.90,2.22)(3)
各对均值差异的置信区间：
(μAi−μAj−tα2(n−r)(1nAi+1nAj)MSE,μAi−μAj+tα2(n−r)(1nAi+1nAj)MSE)i,j=Fe3+,Fe2+,i≠j(4)\left(\mu_{Ai}-\mu_{Aj}-t_\frac{\alpha}{2}\left(n-r\right)\sqrt{\left(\frac{1}{n_{Ai}}+\frac{1}{n_{Aj}}\right)MS_E},\mu_{Ai}-\mu_{Aj}+t_\frac{\alpha}{2}\left(n-r\right)\sqrt{\left(\frac{1}{n_{Ai}}+\frac{1}{n_{Aj}}\right)MS_E}\right)i,j=Fe^{3+},Fe^{2+},i\neq j \tag 4 (μAi−μAj−t2α(n−r)(nAi1+nAj1)MSE,μAi−μAj+t2α(n−r)(nAi1+nAj1)MSE)i,j=Fe3+,Fe2+,i=j(4)
各对均值差异的置信区间为：
μFe3+−μFe2+:(−0.50,−0.05)(5)\mu_{Fe^{3+}}-\mu_{Fe^{2+}}:\left(-0.50,-0.05\right) \tag 5 μFe3+−μFe2+:(−0.50,−0.05)(5)
因为整个区间在0的左边，则我们可以95%的置信度断言\mu_i小于\mu_j，根据上述结果，可以得出以下结论：
μFe3+<μFe2+(6)\mu_{Fe^{3+}}<\mu_{Fe^{2+}} \tag 6 μFe3+<μFe2+(6)
因此可以认为，二价铁离子对百分比的影响更为显著。
对B有：
各水平均值的置信区间：
(μBi−tα2(n−s)MSEnBi,μBi+tα2(n−s)MSEnBi),i=H,M,L(7)\left(\mu_{Bi}-t_\frac{\alpha}{2}\left(n-s\right)\sqrt{\frac{MS_E}{n_{Bi}}},\mu_{Bi}+t_\frac{\alpha}{2}\left(n-s\right)\sqrt{\frac{MS_E}{n_{Bi}}}\right),i=H,M,L \tag 7 (μBi−t2α(n−s)nBiMSE,μBi+t2α(n−s)nBiMSE),i=H,M,L(7)
经查表之后有：
tscore=1.982815t_{score}=1.982815 tscore=1.982815
各均值的置信区间为：
μH:(1.23,1.61)μM:(1.80,2.19)μL:(2.15,2.54)(8)\mu_H:\left(1.23,1.61\right) \\ \mu_M:\left(1.80,2.19\right) \\ \mu_L:\left(2.15,2.54\right) \tag 8 μH:(1.23,1.61)μM:(1.80,2.19)μL:(2.15,2.54)(8)
各对均值差异的置信区间：
(μBi−μBj−tα2(n−s)(1nBi+1nBj)MSE,μBi−μBj+tα2(n−s)(1nBi+1nBj)MSE),i,j=H,M,L,i≠j(9)\left(\mu_{Bi}-\mu_{Bj}-t_\frac{\alpha}{2}\left(n-s\right)\sqrt{\left(\frac{1}{n_{Bi}}+\frac{1}{n_{Bj}}\right)MS_E},\mu_{Bi}-\mu_{Bj}+t_\frac{\alpha}{2}\left(n-s\right)\sqrt{\left(\frac{1}{n_{Bi}}+\frac{1}{n_{Bj}}\right)MS_E}\right),i,j=H,M,L,i\neq j \tag 9 (μBi−μBj−t2α(n−s)(nBi1+nBj1)MSE,μBi−μBj+t2α(n−s)(nBi1+nBj1)MSE),i,j=H,M,L,i=j(9)
各对均值差异的置信区间为：
μH−μM:(−0.85,−0.30)μM−μL:(−0.62,−0.07)μH−μL:(−1.20,−0.65)(10)\mu_H-\mu_M:\left(-0.85,-0.30\right) \\ \mu_M-\mu_L:\left(-0.62,-0.07\right) \\ {\mu_H-\mu}_L:\left(-1.20,-0.65\right) \tag {10} μH−μM:(−0.85,−0.30)μM−μL:(−0.62,−0.07)μH−μL:(−1.20,−0.65)(10)
因为整个区间在0的左边，则我们可以95%的置信度断言μi\mu_iμi小于μj\mu_jμj，根据上述结果，可以得出以下结论：
μH<μM<μL(11)\mu_H<\mu_M<\mu_L \tag {11} μH<μM<μL(11)
因此可以认为，剂量越小对百分比的影响越显著。

alpha = 0.05
n = len(df['C'])
MSE = anovat['mean_sq']['Residual']
# A的置信区间的求解
# print(df['A'].value_counts())
uA = [df['C'][0:54].mean(),df['C'][54:108].mean()]
r = len(df['A'].unique())
nA = [54,54]
t_score = stats.t.isf(alpha / 2, df = (n-r) )
for i in range(r):lower_limit = uA[i] - t_score*np.sqrt(MSE/nA[i])upper_limit = uA[i] + t_score*np.sqrt(MSE/nA[i])print(str((1-alpha)*100)+ '%% Confidence Interval: ( %.2f, %.2f)' % (lower_limit, upper_limit))# 各对均值差异的置信区间
for i in range(r):for j in range(i+1,r):lower_limit = uA[i] - uA[j] - t_score*np.sqrt(MSE*(1/nA[i]+1/nA[j]))upper_limit = uA[i] - uA[j] + t_score*np.sqrt(MSE*(1/nA[i]+1/nA[j]))print(str((1-alpha)*100)+ '%% Confidence Interval: ( %.2f, %.2f)' % (lower_limit, upper_limit))# B的置信区间的求解
# print(df['B'].value_counts())
uB = []
for i in list(df['B'].unique()):sum = 0for j in range(len(df['C'])):if df['B'][j] == i:sum += df['C'][j]uB.append(sum/36)
s = len(df['B'].unique())
nB = [36,36,36]
t_score = stats.t.isf(alpha / 2, df = (n-s) )
for i in range(s):lower_limit = uB[i] - t_score*np.sqrt(MSE/nB[i])upper_limit = uB[i] + t_score*np.sqrt(MSE/nB[i])print(str((1-alpha)*100)+ '%% Confidence Interval: ( %.2f, %.2f)' % (lower_limit, upper_limit))# 各对均值差异的置信区间
for i in range(s):for j in range(i+1,s):lower_limit = uB[i] - uB[j] - t_score*np.sqrt(MSE*(1/nB[i]+1/nB[j]))upper_limit = uB[i] - uB[j] + t_score*np.sqrt(MSE*(1/nB[i]+1/nB[j]))print(str((1-alpha)*100)+ '%% Confidence Interval: ( %.2f, %.2f)' % (lower_limit, upper_limit))