从原理和公式出发:python实现One_Way_ANOVA

2023-11-30 21:36

本文主要是介绍从原理和公式出发:python实现One_Way_ANOVA,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

文章目录

  • 目的:python实现one way ANOVA 单因素方差分析
    • 1. 代码流程
    • 2. python代码实现
      • 0 主要的函数
      • 1 加载数据
      • 2 查看数据统计结果
      • 3 数据处理及可视化
      • 4 方差分析
        • 4.1 模型拟合
        • 4.2 单因素方差分析
      • 5 Post Hoc t-test组间比较分析
      • 6 根据定义自行分解计算对比调用函数的结果
      • 7 获取F分布对应的P值
    • 3. 方差分析公式及原理参考

目的:python实现one way ANOVA 单因素方差分析

方差分析 (Analysis of Variance, ANOVA) 是一种用于比较两个或多个样本均值
是否有显著差异的统计方法。它通过比较组间变异组内变异的大小关系,来判
样本均值是否有显著差异
请添加图片描述

1. 代码流程

a. 通过调用python的包来进行方差分析
b. 根据公式进行方差分析
c. 对比两种方法(工具包与手算)的结果,结果发现:一致
d. 最后附上方差分析的原理和计算公式

2. python代码实现

0 主要的函数

在这里插入图片描述

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as snsfrom scipy import stats                       # 里面有方差齐性检验方差
from statsmodels.formula.api import ols       # 最小二乘法拟合
from statsmodels.stats.anova import anova_lm  # 方差分析import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning)

在这里插入图片描述

1 加载数据

dataD = {'treat1':[390,410,372,382,None],'treat2':[375,348,354,364,362],'treat3':[413,383,408,None,None]}
data = pd.DataFrame(dataD)
data
treat1treat2treat3
0390.0375413.0
1410.0348383.0
2372.0354408.0
3382.0364NaN
4NaN362NaN

2 查看数据统计结果

data.describe()
treat1treat2treat3
count4.0000005.0000003.000000
mean388.500000360.600000401.333333
std16.11417610.28591316.072751
min372.000000348.000000383.000000
25%379.500000354.000000395.500000
50%386.000000362.000000408.000000
75%395.000000364.000000410.500000
max410.000000375.000000413.000000

3 数据处理及可视化

我们为了方便计算,将所有的数据合成一列,并画一下箱线图看一下

data_melt = data.melt()
data_melt.columns = ['Treat', 'value']
print(data_melt)
     Treat  value
0   treat1  390.0
1   treat1  410.0
2   treat1  372.0
3   treat1  382.0
4   treat1    NaN
5   treat2  375.0
6   treat2  348.0
7   treat2  354.0
8   treat2  364.0
9   treat2  362.0
10  treat3  413.0
11  treat3  383.0
12  treat3  408.0
13  treat3    NaN
14  treat3    NaN
import seaborn as sns
# plt.figure(dpi=600)
sns.boxplot(data = data_melt, x = 'Treat', y = 'value', palette = 'pastel',  # 控制箱子颜色)

请添加图片描述

4 方差分析

4.1 模型拟合
from statsmodels.formula.api import ols                    # 最小二乘法拟合
from statsmodels.stats.anova import anova_lm               # 方差分析
from statsmodels.stats.multicomp import pairwise_tukeyhsd  # post Hoc t_testmodel = ols('value ~C(Treat)', data = data_melt).fit()  # 最小二乘法拟合
# ols模型拟合的参数
print(model.params)# 模型拟合的均值:
# mean_treat1 = 388.5
# mean_treat2 = 388.5 - 27.9 = 360.6
# mean_treat3 = 388.5 + 12.833 = 401.33## 实际的均值:
# mean_treat1 = 388.500000	
# mean_treat2 = 360.600000	
# mean_treat3 = 401.333333
Intercept             388.500000
C(Treat)[T.treat2]    -27.900000
C(Treat)[T.treat3]     12.833333
dtype: float64

无论是普通线性模型还是广义线性模型,预测的都是自变量x取特定值时因变量y的平均值。

print(model.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  value   R-squared:                       0.673
Model:                            OLS   Adj. R-squared:                  0.600
Method:                 Least Squares   F-statistic:                     9.257
Date:                Thu, 30 Nov 2023   Prob (F-statistic):            0.00655
Time:                        16:43:18   Log-Likelihood:                -46.814
No. Observations:                  12   AIC:                             99.63
Df Residuals:                       9   BIC:                             101.1
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
======================================================================================coef    std err          t      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------------
Intercept            388.5000      6.910     56.224      0.000     372.869     404.131
C(Treat)[T.treat2]   -27.9000      9.271     -3.010      0.015     -48.871      -6.929
C(Treat)[T.treat3]    12.8333     10.555      1.216      0.255     -11.044      36.710
==============================================================================
Omnibus:                        0.469   Durbin-Watson:                   2.839
Prob(Omnibus):                  0.791   Jarque-Bera (JB):                0.509
Skew:                           0.080   Prob(JB):                        0.775
Kurtosis:                       2.004   Cond. No.                         3.80
==============================================================================Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
4.2 单因素方差分析
anova_table = anova_lm(model, type = 2)                 # 方差分析
pd.DataFrame(anova_table)                               # 查看方差分析结果
dfsum_sqmean_sqFPR(>F)
C(Treat)2.03536.0500001768.0250009.2573930.006547
Residual9.01718.866667190.985185NaNNaN

参数解释:

  1. df: degree of freedom 自由度:自由度是你现有数据中包含的可能性。

    • 例1:为什么当你求全班50个同学身高总和的时候自由度是49?
      因为如果把学号1到学号49的同学身高全都固定下来, 比方说总和为83米,而全班50个人的身高总和,为84.7米,
      那么第50个同学的身高必须为1.7米,没有任何自由变化的余地。
    • 例2:假设有三个数字,A+B+C=10,假设A为任意数取8,B为任意数取2,那么c就是0 这个等式才能成立,这个等式当中有三个未知量,
      但是可以有自由变换的数值只有两个,所以这个式子自由度就是3-1=2,所以样本当中能自由变化的数据个数叫做自由度。

    对于本案例:

    • K = 3: 共有3组,treat1,treat2,treat3。
    • n = 12: 3组共有12个测量值
    • df_C(Treat): K-1 = 3-1=2
    • df_Residual: n-k = 12-3 = 9
  2. sum_sq: error sum of square 误差平方和

    • 表示实验误差大小的偏差平方和,在相同的条件下各次测定值xi对测定平均值x的偏差平方后再加和∑(xi-x)2
  3. mean_sq: Mean Squared Error 均方误差

    • 均方误差是指参数估计值与参数真值之差平方的期望值
  • 可以发现,treat组间存在显著性差异,p=0.006547 < 0.5
    因为对比的组别超过三个,并且呈现出显著性差异,所以考虑使用事后检验(post hoc)进一步对比具体两两组别间的差异情况。

5 Post Hoc t-test组间比较分析

print(pairwise_tukeyhsd(data_melt['value'], data_melt['Treat']))
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===============================================
group1 group2 meandiff p-adj lower upper reject
-----------------------------------------------
treat1 treat2      nan   nan   nan   nan  False
treat1 treat3      nan   nan   nan   nan  False
treat2 treat3      nan   nan   nan   nan  False
-----------------------------------------------

这里我们发现:

  • Q:单因素方差分析结果显著,但事后t检验两两比较均不显著,这样的结果合理吗?
  • A:合理,方差分析结果显著只说明组间可能存在显著差异,到底有无显著差异还要看事后比较

6 根据定义自行分解计算对比调用函数的结果

# 全部的算数平均数为:380.083
mean_all = ((390+410+372+382+375+348+354+364+362+413+383+408)/12)
# 3个品种的算数平均数分别为:388.5,360.6,401.33
print('总平均:\n', mean_all)
mean_k = data.mean(axis = 0)
print('组平均:\n', mean_k)
总平均:380.0833333333333
组平均:treat1    388.500000
treat2    360.600000
treat3    401.333333
dtype: float64
SS_A: 表示因素A的各水平之间的差异带来的影响,又被称为组间偏差。
SS_e: 表示随机误差的影响,又被称为组内偏差。
''' 组间方差 ''' 
SS_A = 4*(388.5-380.083)**2 + 5*(360.6-380.083)**2 + 3*(401.333-380.083)**2
df_A = 3-1  # 共3个水平
MS_A = SS_A/df_A
print("C(Treat) sum_sq:", SS_A)
print("C(Treat) mean_sq:", MS_A)''' 组内方差 ''' 
SS_e = (390-388.5)**2 + (410-388.5)**2 + (372-388.5)**2+(382-388.5)**2+\(375-360.6)**2+(348-360.6)**2+(354-360.6)**2+(364-360.6)**2+(362-360.6)**2+\(413-401.3)**2+(383-401.3)**2+(408-401.3)**2
df_e = 12-3  # 自由度:共12个样本,3个水平
MS_e = SS_e/df_e
print("Residual sum_sq:", SS_e)
print("Residual mean_sq:", MS_e)
C(Treat) sum_sq: 3536.007500999999
C(Treat) mean_sq: 1768.0037504999996Residual sum_sq: 1718.8700000000003
Residual mean_sq: 190.9855555555556
## 计算F
F = MS_A / MS_e
print("MS_A:", MS_A)
print("MS_e:", MS_e)
print("F:", F)
MS_A: 1768.0037504999996
MS_e: 190.9855555555556
F: 9.257264222716081

7 获取F分布对应的P值

from scipy.stats import f                           #导入f
PR = f.sf(F, df_A, df_e)
print(PR)m = df_A      #设置自由度m
n = df_e      #设置自由度n
alpha=0.05    #设置alphaa=f.ppf(q=alpha, dfn=m, dfd=n)                      #单侧左分位点
b=f.isf(q=alpha, dfn=m, dfd=n)                      #单侧右分位点
print('单侧左、右分位点:a=%.4f, b=%.4f'%(a, b))a1, b1=f.interval(1-alpha, dfn=m, dfd=n)              #双侧分位点
print('双侧左、右分位点:a=%.4f, b=%.4f'%(a1, b1))## 可以发现 MS_A/MS_e = 9.25 > Fm,n(0.05) = 5.7147  拒绝原假设H0,即存在组间差异
0.0065472957497462216
单侧左、右分位点:a=0.0516, b=4.2565
双侧左、右分位点:a=0.0254, b=5.7147

3. 方差分析公式及原理参考

请添加图片描述
在这里插入图片描述

原理参考:https://zhuanlan.zhihu.com/p/33357167 (方差分析的好文章)

这篇关于从原理和公式出发:python实现One_Way_ANOVA的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/438649

相关文章

pandas中位数填充空值的实现示例

《pandas中位数填充空值的实现示例》中位数填充是一种简单而有效的方法,用于填充数据集中缺失的值,本文就来介绍一下pandas中位数填充空值的实现,具有一定的参考价值,感兴趣的可以了解一下... 目录什么是中位数填充?为什么选择中位数填充?示例数据结果分析完整代码总结在数据分析和机器学习过程中,处理缺失数

Golang HashMap实现原理解析

《GolangHashMap实现原理解析》HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持高效的插入、查找和删除操作,:本文主要介绍GolangH... 目录HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持

Pandas使用AdaBoost进行分类的实现

《Pandas使用AdaBoost进行分类的实现》Pandas和AdaBoost分类算法,可以高效地进行数据预处理和分类任务,本文主要介绍了Pandas使用AdaBoost进行分类的实现,具有一定的参... 目录什么是 AdaBoost?使用 AdaBoost 的步骤安装必要的库步骤一:数据准备步骤二:模型

使用Pandas进行均值填充的实现

《使用Pandas进行均值填充的实现》缺失数据(NaN值)是一个常见的问题,我们可以通过多种方法来处理缺失数据,其中一种常用的方法是均值填充,本文主要介绍了使用Pandas进行均值填充的实现,感兴趣的... 目录什么是均值填充?为什么选择均值填充?均值填充的步骤实际代码示例总结在数据分析和处理过程中,缺失数

Java对象转换的实现方式汇总

《Java对象转换的实现方式汇总》:本文主要介绍Java对象转换的多种实现方式,本文通过实例代码给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要的朋友参考下吧... 目录Java对象转换的多种实现方式1. 手动映射(Manual Mapping)2. Builder模式3. 工具类辅助映

如何使用 Python 读取 Excel 数据

《如何使用Python读取Excel数据》:本文主要介绍使用Python读取Excel数据的详细教程,通过pandas和openpyxl,你可以轻松读取Excel文件,并进行各种数据处理操... 目录使用 python 读取 Excel 数据的详细教程1. 安装必要的依赖2. 读取 Excel 文件3. 读

Go语言开发实现查询IP信息的MCP服务器

《Go语言开发实现查询IP信息的MCP服务器》随着MCP的快速普及和广泛应用,MCP服务器也层出不穷,本文将详细介绍如何在Go语言中使用go-mcp库来开发一个查询IP信息的MCP... 目录前言mcp-ip-geo 服务器目录结构说明查询 IP 信息功能实现工具实现工具管理查询单个 IP 信息工具的实现服

Python的time模块一些常用功能(各种与时间相关的函数)

《Python的time模块一些常用功能(各种与时间相关的函数)》Python的time模块提供了各种与时间相关的函数,包括获取当前时间、处理时间间隔、执行时间测量等,:本文主要介绍Python的... 目录1. 获取当前时间2. 时间格式化3. 延时执行4. 时间戳运算5. 计算代码执行时间6. 转换为指

利用Python调试串口的示例代码

《利用Python调试串口的示例代码》在嵌入式开发、物联网设备调试过程中,串口通信是最基础的调试手段本文将带你用Python+ttkbootstrap打造一款高颜值、多功能的串口调试助手,需要的可以了... 目录概述:为什么需要专业的串口调试工具项目架构设计1.1 技术栈选型1.2 关键类说明1.3 线程模

SpringBoot基于配置实现短信服务策略的动态切换

《SpringBoot基于配置实现短信服务策略的动态切换》这篇文章主要为大家详细介绍了SpringBoot在接入多个短信服务商(如阿里云、腾讯云、华为云)后,如何根据配置或环境切换使用不同的服务商,需... 目录目标功能示例配置(application.yml)配置类绑定短信发送策略接口示例:阿里云 & 腾