从原理和公式出发: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 单因素方差分析

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_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个水平
print("Residual sum_sq:", 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个水平
print("C(Treat):", SS_A)
Residual sum_sq: 1718.8700000000003
C(Treat): 3536.007500999999
## 计算F
MS_A = SS_A / df_A
MS_e = SS_e / df_eF = 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. 方差分析公式及原理参考

SS_e: 随机误差的影响,又被称为组内偏差;
SS_A: 另一部分表示因素A的各水平之间的差异带来的影响,又被称为组间偏差。
原理参考:https://zhuanlan.zhihu.com/p/33357167 (方差分析的好文章)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/234477.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

从零开始:打造自己的抖音核销工具小程序

对于商家而言&#xff0c;如何高效核销活动中的抖音优惠券成为一项挑战。在这篇文章中&#xff0c;我们将探讨如何从零开始&#xff0c;打造一个个性化、高效的抖音核销工具小程序。 第一步&#xff1a;明确需求和目标 在动手之前&#xff0c;我们需要明确自己的需求和目标。…

CAD精品Eyeshot Fem 2023.3.630 -2023-11-05 Crack

2023.3.630 更新25天前 分享 跟随还没有人关注 改进的 Brep.TransformBy() 方法修复了工具栏内存泄漏修复了 glTF 材质导出期间的异常改进了 glTF 材质金属粗糙度设置修复了渐进式绘图和剪辑平面的错误在 Workspace.UseShaders 属性设置器中添加了缺少的 RenderContext.MakeCur…

JavaScript 的发展史你真的了解吗?

在不断发展的软件开发领域中&#xff0c;很少有编程语言像 JavaScript 一样产生深远的影响。它起初只是一种简单的脚本语言&#xff0c;但如今已成为现代 Web 的驱动力量&#xff0c;改变了应用构建和体验的方式。本文将带你沿着时间线&#xff0c;穿越 JavaScript 的演进历程&…

ffmpeg开发 环境配置

ffmpeg开发简图 1 下载ffmpeg开发包 https://ffmpeg.org/download.html 包含三个版本&#xff1a;Static、Shared以及Dev Static --- 包含3个应用程序&#xff1a;ffmpeg.exe , ffplay.exe , ffprobe.exe&#xff0c;体积都很大&#xff0c;相关的DLL已经被编译到exe里面去…

ClassCMS2.4漏洞复现

ClassCMS2.4漏洞复现 环境搭建 任意文件下载漏洞复现 漏洞成因 ClassCMS2.4漏洞复现 CMS源码在附件中 环境搭建 使用phpstudy2016搭建web环境&#xff0c;php版本为5.5 安装CMS 这里选择Mysql数据库进行安装 用户名和密码都写默认的admin方便记忆 输入完成后点击安装 点…

Redis学习文档

目录 一、概念1、特征2、关系型数据库和非关系型数据库的区别3、键的结构4、Redis的Java客户端5、缓存更新策略5.1、概念5.2、代码 6、缓存穿透6.1、含义6.2、解决办法6.3、缓存空值代码举例6.4、布隆过滤器代码举例 7、缓存击穿7.1、概念7.2、解决办法7.3、互斥锁代码举例7.4、…

P28 C++ 对象的生存周期(栈的作用域生存周期)

前言 本期的主题是栈作用域中对象的生存期&#xff0c;通俗来讲&#xff0c;就是讨论对象是如何在栈上生存的。 这章内容整体分为两部分。 第一部分是&#xff0c;你必须理解栈上的东西是如何存在的&#xff0c;这样你才能真正写出能正常工作的代码。第二部分是&#xff0c;一…

防火墙补充NAT

目录 1.iptables保存规则 2.自定义链 3.NAT NAT的实现分为下面类型&#xff1a; SNAT实验操作 DNAT实验操作 1.iptables保存规则 永久保存方法一&#xff1a; iptables -save > /data/iptables_rule //输出重定向备份 iptables -restore < /data/iptables_r…

修复尺寸显示操作说明

有时我们在绘制草图/工程图中标注尺寸时会显示<MOD-DIAM>符号&#xff0c;如下图所示&#xff1a; 此符号通常在SolidWorks软件更新后会显示出来&#xff0c;Windows更新后可能也会出现此符号。 解决方法&#xff1a; 1. 打开软件-系统选项-文件位置-符号图库文件 2. 选…

MybtisPlus快速开发(从controller到mapper)

创建新项目 写好配置文件 server:port: 8905#配置MP控制台打印日志 mybatis-plus:configuration:log-impl: org.apache.ibatis.logging.stdout.StdOutImplspring:datasource:type: com.zaxxer.hikari.HikariDataSourcedriver-class-name: com.mysql.cj.jdbc.Driverurl: jdbc:m…

行情分析——加密货币市场大盘走势(11.30)

大饼已经形成了底背离&#xff0c;即MACD往下走&#xff0c;而价格还在往上走&#xff0c;这种往往后续会大跌。我们继续把空单拿好&#xff0c;已经持仓的无需加仓。多次上涨都一直不能突破&#xff0c;也说明多空和空军力量都很强&#xff0c;等待后续出方向&#xff0c;而笔…

四大排序方法(java版)

四大排序方法 &#x1f4d1;前言 本文主要是【算法】——常用的排序方法的文章&#xff0c;如果有什么需要改进的地方还请大佬指出⛺️ &#x1f3ac;作者简介&#xff1a;大家好&#xff0c;我是听风与他&#x1f947; ☁️博客首页&#xff1a;CSDN主页听风与他 &#x1f3…