python实现two way ANOVA

＾哪来的＆永远～

已于 2023-12-05 17:53:01 修改

阅读量1.7k

点赞数 17

文章标签： python 概率论

于 2023-11-30 18:56:33 首次发布

本文链接：https://blog.csdn.net/qq_45538220/article/details/134719298

版权

文章目录

目的：用python实现two way ANOVA 双因素方差分析

目的：用python实现two way ANOVA 双因素方差分析

1. python代码实现

在这里插入图片描述

0 双因素方差分析计算公式与结果分析表格

在这里插入图片描述

1 加载python库

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats   

from statsmodels.formula.api import ols                     # 最小二乘法拟合
from statsmodels.stats.anova import anova_lm                # 方差分析
from statsmodels.stats.multicomp import pairwise_tukeyhsd   # post Hoc t_test

2 加载数据

value（month）：治疗2周的各类被试的血糖值
group（PhysicalTherapy）：表示有两种治疗方案，sham组（1）和rTMS组（2）
people（PsychiatricTreatment）：表示有3种被试，年轻健康组（3）、老年健康组（1）、老年患病组（2）

df = pd.read_excel('data//TMS_demoData1.xlsx')
data = pd.DataFrame(df)
data.head(24)

	value	group	people
0	11.0	1	1
1	9.4	1	2
2	12.5	1	3
3	9.6	1	1
4	9.6	1	2
5	11.5	1	3
6	10.8	1	1
7	9.6	1	2
8	10.5	1	3
9	10.5	1	1
10	10.8	1	2
11	12.5	1	3
12	10.5	2	1
13	10.8	2	2
14	10.5	2	3
15	11.5	2	1
16	10.5	2	2
17	11.8	2	3
18	12.0	2	1
19	10.5	2	2
20	11.5	2	3
21	11.8	2	1
22	10.2	2	2
23	11.5	2	3

3 统计样本重复次数，均值和方差，绘制箱线图

data.describe()

	value	group	people
count	24.000000	24.000000	24.000000
mean	10.891667	1.500000	2.000000
std	0.889512	0.510754	0.834058
min	9.400000	1.000000	1.000000
25%	10.500000	1.000000	1.000000
50%	10.800000	1.500000	2.000000
75%	11.500000	2.000000	3.000000
max	12.500000	2.000000	3.000000

fig, ax = plt.subplots(1,2,figsize=(12,6),dpi=600)  # 1行2列的子图
sns.boxplot(x = 'group', y = 'value', data = data, ax = ax[0])
sns.boxplot(x = 'people', y = 'value', data = data, ax = ax[1])

## 可以看出rTMS组的血糖水平sham组的高，因此我们得看这是由于治疗方案引起的还是由于随机误差引起的
## 即从试验结果推断，因素 group 对试验结果有无显著影响，即当 group 取不同水平时试验结果有无显著差别
## 第5步方差分析的结果显示，因素 group 对试验结果无显著影响（p = 0.149），即当 group 取不同水平时试验结果无显著差别

请添加图片描述

4 查看people和group是否存在交互效应

主效应：一个自变量变化时，因变量所出现的变化。
交互效应：反应的是两个或多个自变量对因变量的联合影响，这种影响不能简单的通过自变量的主效应相加获得。

fig, ax = plt.subplots(1,2,figsize=(12,6),dpi=600)  # 1行2列的子图
sns.lineplot(y='value', x = 'people', hue = 'group', palette="tab10", data=data, ax = ax[0])
sns.lineplot(y='value', x = 'group', hue = 'people', palette="tab10", data=data, ax = ax[1])

请添加图片描述

从上图可以看出，people 2和3之间是存在交互效应的，下面可以通过方差分析来检验

5 模型拟合与Two Way ANOVA：双因素方差分析

model = ols('value ~C(group) + C(people) + C(group):C(people)', data = data).fit()
anova_table = anova_lm(model, type = 2)
pd.DataFrame(anova_table)

	df	sum_sq	mean_sq	F	PR(>F)
C(group)	1.0	0.960000	0.960000	2.272189	0.149064
C(people)	2.0	7.485833	3.742917	8.858974	0.002096
C(group):C(people)	2.0	2.147500	1.073750	2.541420	0.106623
Residual	18.0	7.605000	0.422500	NaN	NaN

根据上述结果可以发现：

people组是小于0.05的，存在显著性差异。即people因素对指标value有显著性影响。
group和两者的交互效应是大于0.05的，接受假设，不存在显著性差异，不存在交互效应。
因为people对value存在显著性差异，我们得进行进一步的T检验，查看是那两组之间存在显著性差异。

print(model.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  value   R-squared:                       0.582
Model:                            OLS   Adj. R-squared:                  0.466
Method:                 Least Squares   F-statistic:                     5.015
Date:                Thu, 30 Nov 2023   Prob (F-statistic):            0.00473
Time:                        17:55:19   Log-Likelihood:                -20.264
No. Observations:                  24   AIC:                             52.53
Df Residuals:                      18   BIC:                             59.60
Df Model:                           5                                         
Covariance Type:            nonrobust                                         
================================================================================================
                                   coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------------------
Intercept                       10.4750      0.325     32.231      0.000       9.792      11.158
C(group)[T.2]                    0.9750      0.460      2.121      0.048       0.009       1.941
C(people)[T.2]                  -0.6250      0.460     -1.360      0.191      -1.591       0.341
C(people)[T.3]                   1.2750      0.460      2.774      0.013       0.309       2.241
C(group)[T.2]:C(people)[T.2]    -0.3250      0.650     -0.500      0.623      -1.691       1.041
C(group)[T.2]:C(people)[T.3]    -1.4000      0.650     -2.154      0.045      -2.766      -0.034
==============================================================================
Omnibus:                        1.196   Durbin-Watson:                   2.279
Prob(Omnibus):                  0.550   Jarque-Bera (JB):                1.081
Skew:                          -0.461   Prob(JB):                        0.582
Kurtosis:                       2.518   Cond. No.                         9.77
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

print(model.params)
# 拟合公式
# Yij = 10.475*G(1)*P(1) + 0.9150*G(2) - 0.625*P(2) + 1.275*P(3) - 0.325*G(2)*P(2) -1.4*G(2)*P(3)

Intercept                       10.475
C(group)[T.2]                    0.975
C(people)[T.2]                  -0.625
C(people)[T.3]                   1.275
C(group)[T.2]:C(people)[T.2]    -0.325
C(group)[T.2]:C(people)[T.3]    -1.400
dtype: float64

曲线拟合出来的实为每一种组合的均值：拟合参数验算
*** 无论是普通线性模型还是广义线性模型，预测的都是自变量x取特定值时因变量y的平均值。
因变量y的实际取值与其平均值之差被称为误差项，而误差的分布很大程度上决定了使用什么模型。

# Yij = 10.475*G(1)*P(1) + 0.9150*G(2) - 0.625*P(2) + 1.275*P(3) - 0.325*G(2)*P(2) -1.4*G(2)*P(3)
# GPxx：实际值
# Y_GPxx：预测值
''' Group = 1，People = 1 这个作为截距，后面的每一种组合要加上Intercept'''
Intercept = (11+9.6+10.8+10.5)/4 

''' Group = 1, People = 2 '''
GP12 = (9.4+9.6+9.6+10.8)/4
Y_GP12 = 10.475 - 0.625

''' Group = 1, People = 3 '''
GP13 = (12.5+11.5+10.5+12.5)/4
Y_GP13 = 10.475 + 1.275

''' Group = 2, People = 1 '''
GP21 = (10.5+11.5+12+11.8)/4
Y_GP21 = 10.475 + 0.9105

''' Group = 2, People = 2 '''
GP22 = (10.8+10.5+10.5+10.2)/4
Y_GP22 = 10.475 + 0.9105 - 0.625 - 0.325

''' Group = 2, People = 3 '''
GP23 = (10.5+11.8+11.5+11.5)/4
Y_GP23 = 10.475 + 0.9105 + 1.275 - 1.4

print('Intercept:', Intercept)
print('GP12:',GP12, 'Y_GP12:',Y_GP12)
print('GP13:',GP13, 'Y_GP13:',Y_GP13)
print('GP21:',GP21, 'Y_GP21:',Y_GP21)
print('GP22:',GP22, 'Y_GP22:',Y_GP22)
print('GP23:',GP23, 'Y_GP23:',Y_GP23)

Intercept: 10.475000000000001
GP12: 9.850000000000001 Y_GP12: 9.85
GP13: 11.75 Y_GP13: 11.75
GP21: 11.45 Y_GP21: 11.3855
GP22: 10.5 Y_GP22: 10.435500000000001
GP23: 11.325 Y_GP23: 11.2605

6 单因素多水平 --> 多重比较 post hoc t-tests

pairwise_tukeyhsd is just a wrapper around tukeyhsd method of MultiComparison
print(pairwise_tukeyhsd(data['value'], data['people']))
等价于
from statsmodels.stats.multicomp import MultiComparison
print(MultiComparison(data['value'], data['people']).tukeyhsd())
不选Tudey法，自定义多重比较方法 testfunc
from scipy.stats import mannwhitneyu # Wilcoxon秩和检验，两样本t检验
Tres = MultiComparison(data['value'], data['people']).allpairtest(testfunc = mannwhitneyu)
Tres[0]

print("people因子不同水平的比较结果：", pairwise_tukeyhsd(data['value'], data['people']))
print("###########################\n")
print("group 因子不同水平的比较结果：", pairwise_tukeyhsd(data['value'], data['group']))

print("结果说明: reject=True，说明两组之间有显著性差异。")

people因子不同水平的比较结果： Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
     1      2  -0.7875 0.0935 -1.6876 0.1126  False
     1      3    0.575 0.2635 -0.3251 1.4751  False
     2      3   1.3625 0.0028  0.4624 2.2626   True
---------------------------------------------------
###########################

group 因子不同水平的比较结果： Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
     1      2      0.4 0.2803 -0.3495 1.1495  False
---------------------------------------------------
结果说明: reject=True，说明两组之间有显著性差异。

7 多因素多水平 --> 多重比较 post hoc t-tests

# 这块是因为我开始命名为数字，int相加的结果是int，我要转成str，str相加变成字符串连接
data['group'] = data['group'].replace({1: 'g1'})
data['group'] = data['group'].replace({2: 'g2'})

data['people'] = data['people'].replace({1: 'p1'})
data['people'] = data['people'].replace({2: 'p2'})
data['people'] = data['people'].replace({3: 'p3'})

print(pairwise_tukeyhsd(data['value'], data['people']+data['group']))

Multiple Comparison of Means - Tukey HSD, FWER=0.05
====================================================
group1 group2 meandiff p-adj   lower   upper  reject
----------------------------------------------------
p1g1   p1g2    0.975 0.3202 -0.4857  2.4357  False
p1g1   p2g1   -0.625 0.7491 -2.0857  0.8357  False
p1g1   p2g2    0.025    1.0 -1.4357  1.4857  False
p1g1   p3g1    1.275 0.1082 -0.1857  2.7357  False
p1g1   p3g2     0.85 0.4615 -0.6107  2.3107  False
p1g2   p2g1     -1.6 0.0272 -3.0607 -0.1393   True
p1g2   p2g2    -0.95 0.3462 -2.4107  0.5107  False
p1g2   p3g1      0.3  0.985 -1.1607  1.7607  False
p1g2   p3g2   -0.125 0.9998 -1.5857  1.3357  False
p2g1   p2g2     0.65 0.7185 -0.8107  2.1107  False
p2g1   p3g1      1.9  0.007  0.4393  3.3607   True
p2g1   p3g2    1.475  0.047  0.0143  2.9357   True
p2g2   p3g1     1.25 0.1195 -0.2107  2.7107  False
p2g2   p3g2    0.825 0.4928 -0.6357  2.2857  False
p3g1   p3g2   -0.425 0.9349 -1.8857  1.0357  False

p1g1 p1g2 ：在主效应people的水平为p1的情况下，g1与g2的配对检验结果

8 计算效应量Correlation family： η^2、ω2 （适用于 Correlational data）

神奇的发现：计算方式不同，但是η^2 = ω^2

# η^2 = (F_A X df_A)/(F_A X df_A +df_e)
yitaG = (2.272 * 1)/(2.272 * 1 + 18)
yitaP = (8.86 * 2)/(8.86 * 2 + 18)
yitaGP = (2.5414 * 2)/(2.5414 * 2 + 18)
print('Group的效应量', yitaG)
print('People的效应量',yitaP)
print('GroupxPeople的效应量',yitaGP)

Group的效应量 0.11207576953433307
People的效应量 0.49608062709966405
GroupxPeople的效应量 0.22019858942589288

# ω^2 = sq_A /(sq_A + sq_e)
oumigaG = 0.96/(0.96 + 7.605)
oumigaP = 7.486/(7.486 + 7.605)
oumigaGP = 2.147/(2.147 + 7.605)
print('Group的效应量', oumigaG)
print('People的效应量',oumigaP)
print('GroupxPeople的效应量',oumigaGP)

Group的效应量 0.11208406304728544
People的效应量 0.49605725266715256
GroupxPeople的效应量 0.22015996718621816

2. 双因素方差分析理论和公式

参考：https://zhuanlan.zhihu.com/p/33357167

3. 效应量分析

参考：https://zhuanlan.zhihu.com/p/137779235