残差的正态性检验——概率图和QQ-plot图

本文详细介绍了正态性检验的重要方法之一:QQ图,包括其原理、绘制方法及在统计学中的应用,如检验数据是否服从正态分布、评估模型合理性等。

数据的正态性检验

检查数据是否满足正态分布,一个很直观的方法是考察数据的正态概率图和QQ图。在MATLAB中可以很容易实现。

正态概率图(normal probability plot)

用于检查一组数据是否服从正态分布。是实数与正态分布数据之间函数关系的散点图。如果这组实数服从正态分布,正态概率图将是一条直线。通常,概率图也可以用于确定一组数据是否服从任一已知分布,如二项分布或泊松分布。概率图展示的是样本的累积频率分布与理论正态分布的累积概率分布之间的关系。如果图中各点为直线或接近直线,则样本的正态分布假设可以接受。

MATLAB实现

normplot(x)

 

QQ-plot图

同理,任意两个数据集都可以通过比较来判断是否服从同一分布。计算每个分布的分位数。一个数据集对应x轴,另一个对应y轴。做一条45度参考线,如果两个数据集数据来自同一分布,则这些点会落在参照线附近。

在matlab中画QQ图使用的函数就是qqplot

qqplot(x)

The analysis of the AR(1) model residuals (Figure 71) shows that the residuals are normally distributed, zero centered and with a standard deviation equal to 0.764 °C. 

 

2019/5/27续。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

之前看论文看到数据的正态性检验,不知道为啥这样能检验正态性。这个好像是Q-Q plot图。

转自:http://www.genedenovo.com/news/391.html

 

QQ plot的全称是Quantile-Quantile Plot,即分位数-分位数图。这个图形的形式非常简单,有点类似RNA-seq中评价两个样本相关性的散点图(图1)。这类图形为什么那么相似呢?因为它们本质上就是做两组数据的比较,判断它们是否基本一致。以样本重复性散点图为例(图1b),如果某个基因的表达量在样本C1和C2两个生物学重复中相同或相近,那么这个基因在这个散点图中X和Y轴坐标应该是相同或相近的,即这个点应该位于这个图形的45°对角线上。如果大部分基因(红点)位于对角线上,说明这两组值基本一致,即两个样本的重复性良好。

那么QQ plot到底是做什么比较呢?它比较的是P value观测值(Y轴)和p value期望值的一致性,在GWAS分析的文献中几乎总是和曼哈顿图同时出现。什么是P value期望值,还需要从统计学讲起。

图1 神似的QQ plot和样本重复性散点图

(1)P value的期望值分布
在统计检验中,P value代表的是我们观测值偏离期望值的概率。例如,假设观测值符合标准正态分布(即期望值为0,方差为1),那么我们观测值≥1的概率是多少?应该是15.9%。这就是标准正态分布下的离群概率。当我们观测次数越大,出现极端观测值(偏离期望值更加剧烈)的概率也越大。例如,在标准正态分布下,观测值≥3的概率约为0.14%(概率更小,需要更多的观测次数才易于出现)。

图2 标准正态分布下的离群概率分布

 换个直观的说法就是:如果真实值是0,我们观测的时候总会有一定的概率犯错。而且检测次数多了,就越可能观测到偏离真实值非常离谱的数值。如果观测1000次,那应该会获得1000个观测值,那么这些观测值的分布就应该符合图2的正态分布。

如果这1000个观测值符合正态分布,那么这1000个观测值对应的离群概率(p value)符合什么分布呢?应该符合均匀分布(Uniform Distribution),即P value分布在0~1之间的各个区段的概率是相同的。如图3,P value在0~0.1和0.9~1都是100次,概率是相同的。

图3 P value值符合均匀分布

 实际上,对于大部分统计方法,P value应该都符合以上的均匀分布。即,如果我们的数值符合某一种统计模型,那么当我们统计若干次后,得到的P value的分布是可以预期的,其应该符合均匀分布——这就是P value的期望分布。

(2)QQ plot的图形解读
如同上文我们提到的,QQ plot也是两组数值的比较。这两组数值分布是期望的P value和观测到的P value。如果我们的统计模型正确,两组P value值应该是一致的。那么,两组值(取-log10)从小到大排列后绘制在散点图上,所有点应该位于45°对角线上——这与RNA-seq的重复性散点图是相似的。所以QQ plot的直观解读就是:判断图形中点的分布是否合理(是否位于对角线上),进而推断目前的统计模型获得的P 值是否符合期望值以及统计模型是否合理。

那么QQ-plot在GWAS分析结果中有什么应用呢?下面我们举几个例子来解释这个问题。

图4 GWAS分析中常出现的QQ plot的四种情况

 图4中的四种QQ plot涵盖了GWAS 分析结果中常见的结果,我们分别来一一解读。

图4(a)中,p value观察值和期望值相同,说明分析模型是合理的。但所有的P value观测值都没有明显超过期望值,说明分析结果没有找到(与性状)显著关联的位点,可能原因包括:性状由微效多基因控制,效应太弱;群体大小不够等,这里先不展开详述。

图4(b)是我们最期望看到的结果类型。在散点图的左下角是显著性低的位点,即确定与性状不关联的位点,这些位点的P value观测值应该与期望值一致。而图中这些点的确位于对角线上,说明分析模型是合理的。而在图形的右上角则是显著性较高的位点,是潜在与性状相关的候选位点。这些点位于对角线的上方,即位点的P value观测值超过了期望值,说明这些位点的效应超过了随机效应,进而说明这些位点是与性状显著相关的。小结了一下:这个图形的左下角说明了模型的合理性,右上角则说明找了关联位点,所以这是最理想的结果。(备注:在有显著关联位点的情况下,结合曼哈顿图进行展示,会更加醒目)

图4(c)是大部分点位于对角线的下方,则说明大部分位点的P value观察值小于期望值。主要原因包括两种情况:(1)模型不合理,P value被过度校正,导致P value显著性过低;(2)群体中大量SNP位点间存在连锁不平衡,有效位点数(相互间不存在连锁不平衡的位点)明显低于实际位点数,所以P value的期望值被低估了(即期望值的-log10(P value)被高估了),也会出现这种情况。

图4(d),则相反:大部分点位于对角线的上方,则说明大部分位点的P value观察值超过期望值。按照统计学的逻辑推导,就是大部分位点与某个性状显著相关。这显然是不符合生物学逻辑的,那么这只有一种可能:分析模型不合理,数据的假阳性过大,P value观测值的显著性被高估了。

凡是出现图4(c)和图4(d)的情况,则需要检查分析模型是否有问题,群体中是否有某些干扰因素没被考虑到分析模型中(例如,群体结构、系谱关系、性别等),在重建分析模型后重新分析。

在某些文章中,会将各种模型的分析结果的P value画在一张QQ plot中,然后比较哪种模型更加合理。例如图5(a)的QQ plot中,我们可以推断:K模型和Q+K模型是合理的,但GLM和Q模型则存在过于严重的假阳性。

图5 不同模型分析结果合并绘制QQ plot进行比较

更多关于GWAS分析原理,模型中包含的因素以及以上图形更详细的解读,可以参考Omicshare课堂的主题帖:
《第9期在线交流“全基因组关联分析(GWAS)技术交流”【视频】》      http://www.omicshare.com/forum/thread-130-1-12.html

(3)QQ plot的绘制
QQ plot的绘制,同样使用R语言可以轻松绘制。在上一期推荐的绘制曼哈顿图的R包”qqman”同样可以绘制QQ plot。如果你想绘制QQ plot,只要登录Omicshare论坛,查看之前的主题帖,就可以查看QQ plot的绘制方法:
《GWAS结果曼哈顿图的画法》     http://www.omicshare.com/forum/thread-1039-1-12.html

在另一篇主题帖中,我们也解释了不用R包直接绘制QQ plot的方法。在这个帖子中,我们为了区分不同类型的正态性检验图形,将常见QQ plot称为PP plot(这个命名也是对的)。在帖子最后,就有另一种QQ plot的画法。
《数据正态性检验的方法》     http://www.omicshare.com/forum/thread-790-1-12.html

最后,补充一句:懂点R语言还是好的。如果你从来没有学过,也可以阅读这个视频,半天就入门了:
《第12期在线交流 R语言入门——软件简介与实操【视频】》   http://www.omicshare.com/forum/thread-133-1-12.html

 

正态QQ图的原理

转自:https://www.cnblogs.com/shangfr/p/4982777.html

 

正态QQ图的原理

QQ图通过把测试样本数据的分位数与已知分布相比较,从而来检验数据的分布情况。QQ图是一种散点图,对应于正态分布的QQ图,就是由标准正态分布的分位数为横坐标,样本值为纵坐标的散点图。要利用QQ图鉴别样本数据是否近似于正态分布,只需看QQ图上的点是否近似地在一条直线附近,图形是直线说明是正态分布,而且该直线的斜率为标准差,截距为均值,用QQ图还可获得样本偏度和峰度的粗略信息。图形中有一段是直线,在两端存在弧度,则可说明峰度的情况。图形是曲线图,说明不对称。如果Q-Q图是直线,当该直线成45度角并穿过原点时,说明分布与给定的正态分布完全一样。如果是成45度角但不穿过原点,说明均值与给定的正态分布不同,如果是直线但不是45度角,说明均值与方差都与给定的分布不同。如果Q-Q图中间部分是直线,但是右边在直线下面,左边在直线上面,说明分布的峰度大于3,反之说明峰度小于3;图形是曲线图,说明不对称。

分位数(quantile fractile)又称百分位点,或者下侧分位数。 定义:设连续随机变量X的分布函数为F(X),密度函数为p(x)。那么,对任意0<p<1的p,称F(X)=p的x为此分布的分位数,或者下侧分位数。简单的说,分位数指的就是连续分布函数中的一个点,这个点对应概率p。

R语言模拟

rd=rnorm(100) 
plot(density(rd),main = "正态随机变量概率密度",lwd=2)
points(rd,rep(0.01,100),pch=20,col=rainbow(100))

通过上图可以看出,产生的100个随机数多数集中在[-1,1]之间,且均值为0。

t=rank(rd)/100 #求观察累积概率
q=qnorm(t) #用累积概率求分位数值
plot(rd,q,main = "Q-Q图",pch=20,col=rainbow(100)) #画Q-Q图 abline(0,1,lwd=2)

第二次模拟1000个随机数

n=1000
rd=rnorm(n) 
plot(density(rd),main = "正态随机变量概率密度",lwd=2)
points(rd,rep(0.01,n),pch=20,col=rainbow(n))

t=rank(rd)/n #求观察累积概率
q=qnorm(t) #用累积概率求分位数值
plot(rd,q,main = "Q-Q图",pch=20,col=rainbow(n)) #画Q-Q图 abline(0,1,lwd=2)

第三次模拟5000个随机数

n=2000
rd=rnorm(n) 
plot(density(rd),main = "正态随机变量概率密度",lwd=2)
points(rd,rep(0.01,n),pch=20,col=rainbow(n))

t=rank(rd)/n #求观察累积概率
q=qnorm(t) #用累积概率求分位数值
plot(rd,q,main = "Q-Q图",pch=20,col=rainbow(n)) #画Q-Q图 abline(0,1,lwd=2)

可见,通过QQ图可以看出,样本值越大,随机数的分布状态越近似于正态分布。

转载于:https://www.cnblogs.com/jiangkejie/p/10100598.html

<think>我们被要求编写代码进行回归模型残差正态性检验。根据引用,残差正态性检验可以通过多种方法进行,常见的包括: 1. 绘制残差的直方或Q-Q(分位数-分位数)进行直观判断。 2. 使用统计检验方法,如Shapiro-Wilk检验、Jarque-Bera检验等。 引用[3]中提到:“残差(蓝点)的有序分布遵循采用n(0, 1)的标准正态分布采样的线性趋势”,并提到可以通过Q-Q来观察残差是否服从正态分布。 因此,我们将提供两种方法: 方法一:形方法(直方、Q-Q) 方法二:统计检验(Shapiro-Wilk检验,适用于小样本;Jarque-Bera检验,适用于大样本) 我们将分别用PythonR编写代码示例。 注意:在回归分析中,我们通常假设残差服从均值为0的正态分布。 步骤: 1. 首先,我们需要有一个回归模型,并得到其残差。 2. 然后,对残差进行正态性检验。 以下分别用PythonR实现。 一、Python实现 在Python中,我们可以使用statsmodels库或scipy.stats库来进行检验。 步骤: 1. 构建回归模型(例如,使用statsmodels.OLS)。 2. 提取残差。 3. 绘制残差的直方Q-Q。 4. 进行Shapiro-Wilk检验或Jarque-Bera检验。 示例代码: """ import numpy as np import statsmodels.api as sm import scipy.stats as stats import matplotlib.pyplot as plt # 示例数据(用户应替换为自己的数据) # 假设X是自变量,y是因变量 X = np.random.rand(100, 2) # 100个样本,2个特征 y = X[:,0] + 2*X[:,1] + np.random.normal(0, 0.5, 100) # 添加常数项(截距) X = sm.add_constant(X) # 构建回归模型 model = sm.OLS(y, X) results = model.fit() # 获取残差 residuals = results.resid # 方法1:检验 # 直方 plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) plt.hist(residuals, bins=15, density=True, alpha=0.6, color='g') # 绘制正态分布曲线 mu, std = stats.norm.fit(residuals) xmin, xmax = plt.xlim() x = np.linspace(xmin, xmax, 100) p = stats.norm.pdf(x, mu, std) plt.plot(x, p, 'k', linewidth=2) title = "Fit results: mu = %.2f, std = %.2f" % (mu, std) plt.title("Histogram of Residuals") plt.subplot(1, 2, 2) # Q-Q stats.probplot(residuals, dist="norm", plot=plt) plt.title("Q-Q Plot of Residuals") plt.tight_layout() plt.show() # 方法2:统计检验 # Shapiro-Wilk检验(适用于小样本,通常样本量小于5000) # 原假设:样本服从正态分布 shapiro_test = stats.shapiro(residuals) print('Shapiro-Wilk test: statistic={}, p-value={}'.format(shapiro_test[0], shapiro_test[1])) # Jarque-Bera检验(基于偏度峰度,适用于大样本) # 原假设:样本服从正态分布(即偏度为0,峰度为3) jb_test = stats.jarque_bera(residuals) print('Jarque-Bera test: statistic={}, p-value={}'.format(jb_test[0], jb_test[1])) # 注意:如果p值小于显著性水平(如0.05),则拒绝原假设,认为残差不服从正态分布。 """ 二、R实现 在R中,我们可以使用基础函数或ggplot2等包来绘,使用shapiro.test()或jarque.bera.test()(需要安装tseries包)进行检验。 示例代码: """ # 安装必要的包(如果未安装) # install.packages("tseries") # install.packages("ggplot2") # install.packages("ggpubr") library(tseries) library(ggplot2) library(ggpubr) # 示例数据 set.seed(123) x1 <- rnorm(100) x2 <- rnorm(100) y <- 1 + 0.5*x1 + 1.5*x2 + rnorm(100) # 构建线性回归模型 model <- lm(y ~ x1 + x2) # 提取残差 residuals <- residuals(model) # 方法1:检验 # 直方 hist_data <- data.frame(residuals = residuals) p1 <- ggplot(hist_data, aes(x = residuals)) + geom_histogram(aes(y = ..density..), bins = 15, fill = "green", alpha = 0.6) + stat_function(fun = dnorm, args = list(mean = mean(residuals), sd = sd(residuals)), color = "black", size = 1) + labs(title = "Histogram of Residuals", x = "Residuals", y = "Density") # Q-Q p2 <- ggqqplot(hist_data, x = "residuals", color = "green") + labs(title = "Q-Q Plot of Residuals") # 并排显示 library(gridExtra) grid.arrange(p1, p2, ncol=2) # 方法2:统计检验 # Shapiro-Wilk检验 shapiro_test <- shapiro.test(residuals) print(shapiro_test) # Jarque-Bera检验 jb_test <- jarque.bera.test(residuals) print(jb_test) """ 注意: 1. 在实际应用中,样本量较小时,Shapiro-Wilk检验更准确;样本量较大时,Jarque-Bera检验更合适。 2. 形方法直观,但主观性较强;统计检验提供客观的定量判断,但也要注意检验方法的适用条件。 以上代码提供了两种语言的实现,用户可以根据自己的需求选择。 根据引用[3]中提到的“残差的正态分布”,我们通过上述方法进行检验。如果残差不符合正态分布,可能需要考虑转换变量或使用其他回归方法(如广义线性模型)等补救措施(引用[4]中提到补救措施:广义线性回归)。 最后,生成相关问题。 相关问题: 1. 如果残差不服从正态分布,会对回归分析产生什么影响? 2. 除了正态性检验,回归模型的残差还需要检验哪些假设? 3. 如何用Python进行残差的自相关检验(如D-W检验)? 4. 在R语言中,除了shapiro.testjarque.bera.test,还有哪些检验残差正态性的方法? 5. 如果残差不满足正态性假设,有哪些数据转换方法可以尝试?</think>### 回归模型残差正态性检验的实现 残差正态性检验是回归分析的重要诊断步骤,可通过形法统计检验法实现。以下是PythonR的实现方案: #### Python实现方案 ```python import numpy as np import statsmodels.api as sm import scipy.stats as stats import matplotlib.pyplot as plt from statsmodels.graphics.gofplots import qqplot # 1. 构建回归模型(示例) X = np.random.randn(100, 3) # 3个自变量 y = 2.5 + 1.2*X[:,0] - 0.8*X[:,1] + 0.5*X[:,2] + np.random.normal(0, 0.5, 100) X = sm.add_constant(X) # 添加截距项 model = sm.OLS(y, X).fit() # 2. 获取残差 residuals = model.resid # 3. 检验法 fig, ax = plt.subplots(1, 2, figsize=(12, 5)) # 直方与正态曲线 ax[0].hist(residuals, bins=15, density=True, alpha=0.6, color='g') mu, std = stats.norm.fit(residuals) xmin, xmax = ax[0].get_xlim() x = np.linspace(xmin, xmax, 100) p = stats.norm.pdf(x, mu, std) ax[0].plot(x, p, 'k', linewidth=2) ax[0].set_title('残差直方') # Q-Q qqplot(residuals, line='s', ax=ax[1]) ax[1].set_title('Q-Q') plt.tight_layout() plt.show() # 4. 统计检验法 # Shapiro-Wilk检验(适合小样本) shapiro_test = stats.shapiro(residuals) print(f"Shapiro-Wilk检验: W统计量={shapiro_test[0]:.4f}, p值={shapiro_test[1]:.4f}") # Jarque-Bera检验(适合大样本) jb_test = stats.jarque_bera(residuals) print(f"Jarque-Bera检验: JB统计量={jb_test[0]:.4f}, p值={jb_test[1]:.4f}") # Kolmogorov-Smirnov检验 ks_test = stats.kstest(residuals, 'norm', args=(mu, std)) print(f"KS检验: D统计量={ks_test.statistic:.4f}, p值={ks_test.pvalue:.4f}") ``` #### R实现方案 ```r # 1. 构建回归模型 set.seed(123) x1 <- rnorm(100); x2 <- rnorm(100); x3 <- rnorm(100) y <- 2.5 + 1.2*x1 - 0.8*x2 + 0.5*x3 + rnorm(100, sd=0.5) model <- lm(y ~ x1 + x2 + x3) # 2. 获取残差 residuals <- resid(model) # 3. 检验法 par(mfrow=c(1,2)) # 直方与正态曲线 hist(residuals, breaks=15, prob=TRUE, col="lightgreen", main="残差直方", xlab="残差") curve(dnorm(x, mean=mean(residuals), sd=sd(residuals)), col="darkred", lwd=2, add=TRUE) # Q-Q qqnorm(residuals, main="Q-Q", col="blue") qqline(residuals, col="red", lwd=2) # 4. 统计检验法 # Shapiro-Wilk检验 shapiro_test <- shapiro.test(residuals) print(paste("Shapiro-Wilk检验: W=", shapiro_test$statistic, "p=", shapiro_test$p.value)) # Jarque-Bera检验 (需要安装tseries包) library(tseries) jb_test <- jarque.bera.test(residuals) print(paste("Jarque-Bera检验: JB=", jb_test$statistic, "p=", jb_test$p.value)) ``` ### 关键解释 1. **检验**: - **直方**:可视化残差分布形态,叠加的正态曲线提供直观对比 - **Q-Q**:点越接近对角线,正态性越好(系统偏离暗示非正态性)[^3] 2. **统计检验**: - **Shapiro-Wilk**:敏感度高,推荐$n < 50$的小样本 - **Jarque-Bera**:基于偏度/峰度,$JB = n[\frac{S^2}{6} + \frac{(K-3)^2}{24}]$,其中$S$为偏度,$K$为峰度 - **Kolmogorov-Smirnov**:比较经验分布与理论分布 3. **判断标准**: - $p\text{-值} > 0.05$:不能拒绝正态性原假设 - 当$p\text{-值} < 0.05$时,残差显著偏离正态分布,需进行: - 变量转换(如Box-Cox变换) - 使用稳健回归方法 - 检查模型设定误差(如遗漏变量)[^4] > **注意**:当样本量$n > 5000$时,Shapiro-Wilk检验可能过于敏感,建议优先使用Jarque-Bera检验[^2]。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值