从原理和公式出发: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

相关文章

Python将博客内容html导出为Markdown格式

《Python将博客内容html导出为Markdown格式》Python将博客内容html导出为Markdown格式,通过博客url地址抓取文章,分析并提取出文章标题和内容,将内容构建成html,再转... 目录一、为什么要搞?二、准备如何搞?三、说搞咱就搞!抓取文章提取内容构建html转存markdown

Python获取中国节假日数据记录入JSON文件

《Python获取中国节假日数据记录入JSON文件》项目系统内置的日历应用为了提升用户体验,特别设置了在调休日期显示“休”的UI图标功能,那么问题是这些调休数据从哪里来呢?我尝试一种更为智能的方法:P... 目录节假日数据获取存入jsON文件节假日数据读取封装完整代码项目系统内置的日历应用为了提升用户体验,

SpringBoot3实现Gzip压缩优化的技术指南

《SpringBoot3实现Gzip压缩优化的技术指南》随着Web应用的用户量和数据量增加,网络带宽和页面加载速度逐渐成为瓶颈,为了减少数据传输量,提高用户体验,我们可以使用Gzip压缩HTTP响应,... 目录1、简述2、配置2.1 添加依赖2.2 配置 Gzip 压缩3、服务端应用4、前端应用4.1 N

Java编译生成多个.class文件的原理和作用

《Java编译生成多个.class文件的原理和作用》作为一名经验丰富的开发者,在Java项目中执行编译后,可能会发现一个.java源文件有时会产生多个.class文件,从技术实现层面详细剖析这一现象... 目录一、内部类机制与.class文件生成成员内部类(常规内部类)局部内部类(方法内部类)匿名内部类二、

SpringBoot实现数据库读写分离的3种方法小结

《SpringBoot实现数据库读写分离的3种方法小结》为了提高系统的读写性能和可用性,读写分离是一种经典的数据库架构模式,在SpringBoot应用中,有多种方式可以实现数据库读写分离,本文将介绍三... 目录一、数据库读写分离概述二、方案一:基于AbstractRoutingDataSource实现动态

Python FastAPI+Celery+RabbitMQ实现分布式图片水印处理系统

《PythonFastAPI+Celery+RabbitMQ实现分布式图片水印处理系统》这篇文章主要为大家详细介绍了PythonFastAPI如何结合Celery以及RabbitMQ实现简单的分布式... 实现思路FastAPI 服务器Celery 任务队列RabbitMQ 作为消息代理定时任务处理完整

Python Websockets库的使用指南

《PythonWebsockets库的使用指南》pythonwebsockets库是一个用于创建WebSocket服务器和客户端的Python库,它提供了一种简单的方式来实现实时通信,支持异步和同步... 目录一、WebSocket 简介二、python 的 websockets 库安装三、完整代码示例1.

揭秘Python Socket网络编程的7种硬核用法

《揭秘PythonSocket网络编程的7种硬核用法》Socket不仅能做聊天室,还能干一大堆硬核操作,这篇文章就带大家看看Python网络编程的7种超实用玩法,感兴趣的小伙伴可以跟随小编一起... 目录1.端口扫描器:探测开放端口2.简易 HTTP 服务器:10 秒搭个网页3.局域网游戏:多人联机对战4.

Java枚举类实现Key-Value映射的多种实现方式

《Java枚举类实现Key-Value映射的多种实现方式》在Java开发中,枚举(Enum)是一种特殊的类,本文将详细介绍Java枚举类实现key-value映射的多种方式,有需要的小伙伴可以根据需要... 目录前言一、基础实现方式1.1 为枚举添加属性和构造方法二、http://www.cppcns.co

使用Python实现快速搭建本地HTTP服务器

《使用Python实现快速搭建本地HTTP服务器》:本文主要介绍如何使用Python快速搭建本地HTTP服务器,轻松实现一键HTTP文件共享,同时结合二维码技术,让访问更简单,感兴趣的小伙伴可以了... 目录1. 概述2. 快速搭建 HTTP 文件共享服务2.1 核心思路2.2 代码实现2.3 代码解读3.