In [2]:
Copied!
import os # 修改工作目录
import numpy as np
import pandas as pd
import scipy.stats as stats # 统计函数
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from scipy.stats import t
import math
#import matplotlib.pyplot as plt
#from plotnine import * # ggplot 绘图
#from plotnine.data import mpg
from jupyterquiz import display_quiz # Quiz
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import os # 修改工作目录
import numpy as np
import pandas as pd
import scipy.stats as stats # 统计函数
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from scipy.stats import t
import math
#import matplotlib.pyplot as plt
#from plotnine import * # ggplot 绘图
#from plotnine.data import mpg
from jupyterquiz import display_quiz # Quiz
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
设置数据目录，如下：
In [3]:
Copied!
os.chdir("/Users/lyuni/ECNU_DaSE/Courses/Stat_ML/Experiment/Data")
os.chdir("/Users/lyuni/ECNU_DaSE/Courses/Stat_ML/Experiment/Data")
背景¶
回顾 Tutorial One，我们利用了 One-way ANOVA 模型分析三种饮食方式对减重的影响是否存在差异。而后我们所得到的结论如下：
- 拒绝原假设：我们认为三种饮食方式对减重的影响是存在差异的。检验统计量为$4.8139$，其对应的$p$值为0.011。
- One-way ANOVA模型是符合数据假定的，具体来说，数据满足（1）独立性，（2）方差齐性，以及（3）正态性。
- 三种饮食方式下，体重变化量的估计分别为$3.3, 3.113, 5.037$，而方差的估计为$5.6168$。
我们需要进一步分析哪些饮食方式对减重的影响是有差异的。
数据¶
我们仍使用数据集Data_1，具体形式如下所示。
In [4]:
Copied!
print('Data 1 is shown as follows: \n', pd.read_csv("Data_1.csv"))
print('Data 1 is shown as follows: \n', pd.read_csv("Data_1.csv"))
Data 1 is shown as follows:
gender age height diet.type initial.weight final.weight
0 Female 22 159 A 58 54.2
1 Female 46 192 A 60 54.0
2 Female 55 170 A 64 63.3
3 Female 33 171 A 64 61.1
4 Female 50 170 A 65 62.2
.. ... ... ... ... ... ...
67 Male 35 183 C 83 80.2
68 Male 49 177 C 84 79.9
69 Male 28 164 C 85 79.7
70 Male 40 167 C 87 77.8
71 Male 51 175 C 88 81.9
[72 rows x 6 columns]
任务¶
试用Tukey比较 3 个饮食方式水平下的减重程度两两是否存在差异.
In [5]:
Copied!
## Settings
alpha =0.05
a = 3
m =24
## Load Data
Data = pd.read_csv("Data_1.csv")
print(Data.columns)# Print the column names of Data_1
## Construct a New Dataset
Data = Data[['diet.type','initial.weight','final.weight']] # select some columns from a dataset
Data['weight.loss'] = Data['initial.weight'] - Data['final.weight']
Data = Data.drop(labels=['initial.weight', 'final.weight'], axis = 1) # delete some columns from a dataset
Data.columns = ['Diet_type','Weight_loss']
print(Data)
## Settings
alpha =0.05
a = 3
m =24
## Load Data
Data = pd.read_csv("Data_1.csv")
print(Data.columns)# Print the column names of Data_1
## Construct a New Dataset
Data = Data[['diet.type','initial.weight','final.weight']] # select some columns from a dataset
Data['weight.loss'] = Data['initial.weight'] - Data['final.weight']
Data = Data.drop(labels=['initial.weight', 'final.weight'], axis = 1) # delete some columns from a dataset
Data.columns = ['Diet_type','Weight_loss']
print(Data)
Index(['gender', 'age', 'height', 'diet.type', 'initial.weight',
'final.weight'],
dtype='object')
Diet_type Weight_loss
0 A 3.8
1 A 6.0
2 A 0.7
3 A 2.9
4 A 2.8
.. ... ...
67 C 2.8
68 C 4.1
69 C 5.3
70 C 9.2
71 C 6.1
[72 rows x 2 columns]
具体分析¶
在One-way ANOVA模型中，我们假定了体重的变化量$y_{ij} \sim N(\mu+\alpha_i,\sigma^2)$。欲检验 $$ H_0: \alpha_1 = \alpha_2 = \alpha_3 =0 \quad \text{vs} \quad H_1: \alpha_1,\alpha_2,\alpha_3\text{不全相等}. $$
除了采用第一周代码可以实现One-way ANOVA模型之外，这里我们介绍另一种实现方法。
In [6]:
Copied!
data1 = Data.values
list_type = ["A","B","C"]
groups = [data1[data1[:,0] == x, 1] for x in list_type]
F_stat, F_pVal = stats.f_oneway(groups[0],groups[1],groups[2])
print("The test statistic in One-way ANOVA model is", round(F_stat,4))
print("The p value in One-way ANOVA model is", round(F_pVal,4))
data1 = Data.values
list_type = ["A","B","C"]
groups = [data1[data1[:,0] == x, 1] for x in list_type]
F_stat, F_pVal = stats.f_oneway(groups[0],groups[1],groups[2])
print("The test statistic in One-way ANOVA model is", round(F_stat,4))
print("The p value in One-way ANOVA model is", round(F_pVal,4))
The test statistic in One-way ANOVA model is 4.8139 The p value in One-way ANOVA model is 0.011
无论采用“临界值法”还是根据$p$值来判断，我们均拒绝原假设，即认为这三种饮食方式对体重的变化量有不同的影响。
接下来，我们需要判断哪几对是存在差异的，并且需要量化这些差异是否是显著的？这里我们采用的是理论课程里介绍过的方法——Tukey。在Tukey方法中，需要计算每组样本均值的差，并将其与临界值 $$ c = q_{1-\alpha}(a,df) \cdot \hat{\sigma}/\sqrt{m}. $$ 进行比较。这里所计算的分位数，实际上是$t$化极差统计量的分位数，其中$t$化极差统计量定义为 $$ q(a,df) = \max_i \frac{\bar{y}_{i\cdot} - \mu}{\hat{\sigma}/\sqrt{m}} - \min_i \frac{\bar{y}_{i\cdot} - \mu}{\hat{\sigma}/\sqrt{m}}. $$ 这里$a$表示因子水平数，$m$表示一组内的数据个数，$\hat{\sigma}$是标准差的估计，即$\hat{\sigma}=\sqrt{\hat{\sigma}^2}$。
值得注意的是，这里$q_{1-\alpha}(a,df)$表示$t$化极差统计量的分位数，通常需要查表而得。
备注：上表中的$r$对应于所定义的$a$，而上表中的$f$对应于所定义的df。
我们判断方式是
- 如果两组样本均值的绝对差大于临界值，那么认为这两组样本的均值是不同的；
- 如果两组样本均值的绝对差小于等于临界值，那么认为这两组样本的均值是相同的。
In [7]:
Copied!
Tukey= pairwise_tukeyhsd(endog = Data["Weight_loss"], groups=Data["Diet_type"],alpha=alpha)
print(Tukey)
Tukey= pairwise_tukeyhsd(endog = Data["Weight_loss"], groups=Data["Diet_type"],alpha=alpha)
print(Tukey)
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj lower upper reject
---------------------------------------------------
A B -0.1875 0.9 -1.8263 1.4513 False
A C 1.7375 0.0352 0.0987 3.3763 True
B C 1.925 0.0173 0.2862 3.5638 True
---------------------------------------------------
从上表中，我们可以发现：饮食方式A和C，B和C对体重的改变量是有显著的区别的，而饮食方式A和B对体重的改变量没有显著差异。
可以注意到，除了可以判断每两组均值差的同时置信区间。根据所构建的同时置信区间，我们也可以判断两组样本的均值是否存在显著差异。请回答以下的问题。
In [8]:
Copied!
display_quiz("/Users/lyuni/ECNU_DaSE/Courses/Stat_ML/Experiment/Question/T2_Q1.json")
display_quiz("/Users/lyuni/ECNU_DaSE/Courses/Stat_ML/Experiment/Question/T2_Q1.json")
Q1【点击提示】
请回顾：假设检验与区间估计的关系。Q2¶
In [9]:
Copied!
display_quiz("/Users/lyuni/ECNU_DaSE/Courses/Stat_ML/Experiment/Question/T2_Q2.json")
display_quiz("/Users/lyuni/ECNU_DaSE/Courses/Stat_ML/Experiment/Question/T2_Q2.json")
Q2【点击提示】
请计算出Bonferroni方法所构造的同时置信区间的上下限。