统计分布的检验有很多种，例如KS检验、卡方检验等，从图形的角度来说，我们也可以用QQ图（Quantile-Q...
统计分布的检验有很多种，例如KS检验、卡方检验等，从图形的角度来说，我们也可以用QQ图（Quantile-Quantile Plots）来检查数据是否服从某种分布。QQ图的原理并不复杂：如果一批数据x1, x2, . . . , xn 服从某种理论分布，那么将排序后的数据x(1), x(2), . . . , x(n)和理论分布的分位数q1/n, q2/n, . . . , qn/n去画散点图，得到的n个点应该大致排列在对角线上，因为这两批数字应该大致相等。从另一个角度来看，检验一批数据是否服从某种理论分布，也就是看其经验分布和理论分布是否一致，而排序后的数据x(1), x(2), . . . , x(n)可以看作是经验分布的1/n, 2/n, · · · , n/n分位数，若这些分位数和理论分位数一致，也就说明了经验分布和理论分布相似。
quantile(rnorm(1000),probs = seq(0.1,1,0.2))
#真是分位数
qnorm(seq(0.1,1,0.2))
以上数据的5个分位数和理论分位数都比较接近，读者可以模拟其它分布，例如从卡方分布中生成随机数，看其分位数是多少，与正态分布分位数差异如何。
qqnorm(y, ...)
## Default S3 method:
qqnorm(y, ylim, main = &Normal Q-Q Plot&,
xlab = &Theoretical Quantiles&, ylab = &Sample Quantiles&,
plot.it = TRUE, datax = FALSE, ...)
qqline(y, datax = FALSE, distribution = qnorm,
probs = c(0.25, 0.75), qtype = 7, ...)
qqplot(x, y, plot.it = TRUE, xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)), ...) qqplot()检验的是两批数据的分布是否相同，所以它需要两个数据参数x和y，qqnorm()只需要一个数据参数x，其它设置标签和标题等元素的图形参数此处不再赘述。[b]library(MASS)
par(mfrow=c(1,2))
手机浏览器扫一扫下载神器
安全&& 正版
微信扫一扫关注公众号
要回复文章请先或
微信扫一扫 进一步关注数据正态性检验画图的4种方法
我的图书馆
数据正态性检验画图的4种方法
由于有人问如何使用R进行数据正态性检验，所以周老师干脆写个主题帖解释一下。如果恰好解决了你的问题，请读完后给个好评哟~正态性检验，是很多数据分析前要做的准备性工作。例如，你有组数量性状的表型值，你想先判断其是否符合正态分布，再开展后续的数据分析。最简单的检验方法正态性检验，最简单的方法是使用R语言的shapiro.test命令。如果P value > 5%，则说明数据分布近似正态分布。图形化的比较当然，你还期望有图形化的比较，以便在文章中展示。那么有4种画法。1QQ-plot分位数图功能和原理：检验样本的概率分布是否服从某种理论分布。PP概率图的原理是检验实际累积概率分布与理论累积概率分布是否吻合，若吻合，则散点应围绕在一条直线周围，或者实际概率与理论概率之差分布在对称于以0为水平轴的带内。QQ概率图的原理是检验实际分位数与理论分位数之差分布是否吻合，若吻合，则散点应围绕在一条直线周围，或者实际分位数与理论分位数之差分布在对称于以0为水平轴的带内。QQ概率图以样本的分位数为横轴，以指定理论分布的分位数为纵轴绘制散点图。#install.packages('DAAG')library(DAAG)data(possum)attach(possum) &# 数据准备fpossum # 只分析这些样本中的雌性个体x#将totlngth这个表型均一化，即 标准正态化n plot(qnorm((1:n-0.5)/n),sort(x),col=2,type = 'p',main = 'QQ plot',xlab='Theoretical Quantiles',ylab='Studentized Quantiles' )abline(a=0,b=1,lty=3)图形表示，数据与正态性略有差异，特别是中部区域。2与正态密度函数直接比较library(DAAG)data(possum)attach(possum)fpossum dens xlim ylim mean = mean(totlngth)sd = sd(totlngth)par(mfrow=c(1,2))hist(totlngth,& & breaks=72.5+(0:5)*5,& & xlim = xlim ,& & ylim = ylim ,& & probability = T ,& & xlab = 'total length',& & main = 'A:Breaks at 72.5...')lines(dens,& & & &col = par('fg'),& & & &lty = 2)curve( dnorm(x, mean, sd),& & & & &col = 'red',& & & & &add = T)hist(totlngth,& & & &breaks = 75 + (0:5) * 5 ,& & & &xlim = xlim,& & & &ylim = ylim,& & & &probability = T,& & & &xlab='total length',& & & &main = 'B:Breaks at 75')lines(dens,& & & &col = par('fg'),& & & &lty = 2)curve(dnorm(x,mean,sd),& & & &col = 'red',& & & &add = T)看图直接看和正态密度函数的差异度。这张图在数量性状的文章里最常出现。3经验分布与正态分布函数对比library(DAAG)data(possum)attach(possum)fpossum mean = mean(totlngth)sd = sd(totlngth)x n y plot(x,y,& & & &type = 's',& & & &main = 'Empirical CDF of ')curve(pnorm(x, mean, sd),& & & &col = 'red',& & & &lwd = 2,& & & &add = T)4P-P plot图使用p value画图，常用于比较GWAS分析结果中，观测的P value和理论p value间的差异，代码大概如下：r=read.table('temp_name.txt')& && &# 含位点p值的文件。一行1个位点，p值在第六列；o=-log10(sort(r$V6,decreasing=F))&&# 观测到的p值，对p value列排序，假设p value在第六列e=-log10(ppoints(length(r$V6)))& & # 生成对应的平均分布的p值，为期望值；plot(e,o,pch=20,xlab='Expected~-log10(p)',ylab='Observed~-log10(p)',main='QQ plot',col= 'blue',xlim=c(0,max(e)+0.1),ylim=c(0,max(o)+0.1),bty= 'l ',yaxs= 'i ',xaxs= 'i ',cex=2)&abline(0,1,col= 'red ')#结果如下4种图形化展示方式不知道你心水哪种呢？可以都拿走噢~哈哈~明天见~
TA的最新馆藏[转]&苹果/安卓/wp
积分 9, 距离下一级还需 1 积分
道具: 彩虹炫, 涂鸦板, 雷达卡, 热点灯下一级可获得
道具: 金钱卡
购买后可立即获得
权限: 隐身
道具: 金钱卡, 彩虹炫, 雷达卡, 热点灯, 涂鸦板
请大家都来帮忙，我在学习这个。。。老师要求我会画二维正态分布图。。。。。。
载入中......
满仓中石油
x&-y&-seq(-4,4,length=20)f&-function(x,y){(exp(-0.5*x^2-0.5*y^2))/(2*pi)}z&-outer(x,y,f)persp(x,y,z,theta=45,phi=25,col='lightblue')怎么老师都会留同样的作业。这是我以前用的方法，是个最简单的二元正态分布期望为0，方差为1，相关系数为0如果你需要，可以自己更改函数表达式
十分感谢^_^
heavenicefox
那如果有两个二维正态分布呢？要想把它们的分布图在一张图里画出来，可以么？谢谢~~~
画二维正态分布还可以这样：
x=seq(-5,5,by=0.1) #步长很小时画的图就是黑色的了，因为都是画格子的黑线的颜色
x1=dnorm(x,0,1)&&#dnorm()为正态分布密度函数
z=outer(x1,x1)
persp(x,y,z,theta =30,phi = 25,expand = 0.5,col = &Blue2&)
根据楼上两位的代码做出的图形，画在一张图上了。
(10.91 KB)
21:33:37 上传
本帖最后由 hanson.don 于
21:03 编辑
你完全可以这样做
x&-seq(-10,10,0.1)
y&-dnorm(x,0,1)
y1&-dnorm(x,0,2)
plot(x,y,type=&l&);lines(x,y1)
21:03:59 上传
欲说还休！！！
谢谢楼上。正需要呢
这是用R语言作的，用sas不能吗？我也正在学习中，谢谢楼主
谢谢！！！！
无限扩大经管职场人脉圈！每天抽选10位免费名额，现在就扫& 论坛VIP& 贵宾会员& 可免费加入
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
如有投资本站或合作意向，请联系（010-）；
邮箱：service@pinggu.org
投诉或不良信息处理：（010-）
京ICP证090565号
论坛法律顾问：王进律师R语言（12）
前言：本次建模过程是基于RedHat6.8或者CentOS6.8，R3.1.2,Rstudio-server
关于R3.1.2，Rstudio-server的整个配置，原始数据（已经脱敏处理，不涉及泄密，如有侵权，请随时联系）以及本分析的源码均放置在GitHub上,通过访问
数据导入：
install.packages("rJava",dependencies=TRUE)
install.packages("xlsx",dependencies=TRUE)
install.packages("corrplot",dependencies=TRUE)
install.packages("leaps",dependencies=TRUE)
install.packages("lmtest",dependencies=TRUE)
library(rJava)
setwd("/home/steven/workstation")
library(xlsx)
src &- read.xlsx("/usr/local/workstation/test.xls",1,encoding="UTF-8")
数据处理：
src &- src[,-c(1)]
src &- src[,c(1:6,8:10,7)]
attach(src)
names(src)
mean=sapply(src,mean)
max=sapply(src,max)
min=sapply(src,min)
median=sapply(src,median)
sd=sapply(src,sd)
cbind(mean,max,min,median,sd)
建立回归模型步骤：
#1、参数全部默认情况下的相关系数图
#混合方法之上三角为圆形，下三角为黑色数字
library(corrplot)
corr &- cor(src[,1:10])
corrplot(corr = corr,order="AOE",type="upper",tl.pos="tp")
corrplot(corr = corr,add=TRUE, type="lower",
method="number",order="AOE", col="black",
diag=FALSE,tl.pos="n", cl.pos="n")
输出结果如下：
980.00000 1.
1.00000 1.
7.00000 3.
2.00000 4.
3.00000 9.
1.00000 4.
95.00923 3.
#1、参数全部默认情况下的相关系数图
#混合方法之上三角为圆形，下三角为黑色数字
library(corrplot)
corr &- cor(src[,1:10])
corrplot(corr = corr,order="AOE",type="upper",tl.pos="tp")
corrplot(corr = corr,add=TRUE, type="lower",
method="number",order="AOE", col="black",
diag=FALSE,tl.pos="n", cl.pos="n")
结果如下：
#2.画相关图选择回归方程的形式
plot(target~cg);abline(lm(target~cg))
plot(target~h);abline(lm(target~h))
plot(target~type);abline(lm(target~type))
plot(target~tjg);abline(lm(target~tjg))
plot(target~ggjg);abline(lm(target~ggjg))
plot(target~area);abline(lm(target~area))
plot(target~cs);abline(lm(target~cs))
plot(target~yyjg);abline(lm(target~yyjg))
plot(target~rs);abline(lm(target~rs))
#3.do regression and check results
dim(src)[1]
lm.test&-lm(target~rs+cs+area+type+cg+h+tjg+ggjg+yyjg,data=src)
summary(lm.test)
结果如下：
Residuals:
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
6.507 1.22e-10 ***
-5.617e+02
-3.414e+04
-7.923 6.20e-15 ***
-1.846e+03
& 2e-16 ***
-1.198e+01
-4.868 1.31e-06 ***
6.625 5.66e-11 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12980 on 993 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 138.4 on 9 and 993 DF,
p-value: & 2.2e-16
#4.delete variable which is not significant(rs,area)
lm.test&-lm(target~cs+type+cg+h+tjg+ggjg+yyjg,data=src)
summary(lm.test)
Residuals:
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
6.991 4.99e-12 ***
-7.947 5.14e-15 ***
& 2e-16 ***
-5.279 1.59e-07 ***
7.068 2.96e-12 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12980 on 995 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 177.5 on 7 and 995 DF,
p-value: & 2.2e-16
#4.1.use residual analysis delete outlier points
plot(lm.test,which=1:4)
src = src[-c(12,765,788,790),]
dim(src)[1]
结果如下：
得到的四个图依次为：
4.1普通残差与拟合值的残差图
4.2正态QQ的残差图（若残差是来自正态总体分布的样本，则QQ图中的点应该在一条直线上）
4.3标准化残差开方与拟合值的残差图（对于近似服从正态分布的标准化残差，应该有95%的样本点落在[-2,2]的区间内。这也是判断异常点的直观方法）
4.4cook统计量的残差图（cook统计量值越大的点越可能是异常值，但具体阀值是多少较难判别）
从图中可见，12,765,788,790三个样本存在异常，需要剔除。
Residuals:
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
6.991 4.99e-12 ***
-7.947 5.14e-15 ***
& 2e-16 ***
-5.279 1.59e-07 ***
7.068 2.96e-12 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12980 on 995 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 177.5 on 7 and 995 DF,
p-value: & 2.2e-16
#5.异方差检验
#5.1GQtest，H0（误差平方与自变量，自变量的平方和其交叉相都不相关），
#p值很小时拒绝H0，认为上诉公式有相关性，存在异方差
src.test&-residuals(lm.test)
library(lmtest)
gqtest(lm.test)
#5.2BPtest,H0(同方差),p值很小时认为存在异方差
bptest(lm.test)
结果如下：
Goldfeld-Quandt test
studentized Breusch-Pagan test
BP = 93.9696, df = 7, p-value & 2.2e-16
两个检验的p值都很小时认为存在异方差，需要修正异方差
#6.修正异方差
#修正的方法选择FGLS即可行广义最小二乘
#6.1修正步骤
#6.1.1将y对xi求回归，算出res--u
#6.1.2计算log(u^2)
#6.1.3做log(u^2)对xi的辅助回归 log(u^2),得到拟合函数g=b0+b1x1+..+b2x2
#6.1.4计算拟合权数1/h=1/exp(g)，并以此做wls估计
lm.test2&-lm(log(resid(lm.test)^2)~cs+type+cg+h+tjg+ggjg+yyjg,data=src)
lm.test3&-lm(target~cs+type+cg+h+tjg+ggjg+yyjg,weights=1/exp(fitted(lm.test2)),data=src)
summary(lm.test3)
结果如下：
Weighted Residuals:
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
4.607 4.62e-06 ***
-5.571 3.26e-08 ***
& 2e-16 ***
-6.638 5.20e-11 ***
7.548 9.96e-14 ***
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.868 on 995 degrees of freedom
Multiple R-squared:
Adjusted R-squared:
F-statistic: 361.6 on 7 and 995 DF,
p-value: & 2.2e-16
#7.1计算解释变量相关稀疏矩阵的条件数，&100多重共线性程度很小，100&&1000较强，&1000严重
结果如下：
[1] 160.9175
K&100 and K&1000，说明共线性较强，接下来找出共线性强的变量
[1] 3.4192494 1.3030312 1.1563047 0.9199288 0.9041041 0.7562702 0.4205589 0.1011783
[9] 0.0193744
0.2026552 .
0.4116265 . . .
0. .1293286 .
0. . .1278988
0. . .2352263
0. . .2471633 . .
0.4469450 .
#8.修正多重共线性---逐步回归法
#ps2：step中可进行参数设置：
#direction=c("both","forward","backward")来选择逐步回归
#的方向，默认both，forward时逐渐增加解释变两个数，backward则相反。
step(lm.test)
结果如下：
AIC=19007.25
target ~ cs + type + cg + h + tjg + ggjg + yyjg
1 1.2269e+08 1.6778e+11 19006
1 6.6200e+08 1.6832e+11 19009
1 1.5794e+09 1.6924e+11 19015
1 4.6960e+09 1.7235e+11 19033
1 8.4171e+09 1.7608e+11 19054
1 1.0643e+10 1.7830e+11 19067
1 2.2201e+10 1.8986e+11 19130
AIC=19005.98
target ~ cs + type + h + tjg + ggjg + yyjg
1 5.9084e+08 1.6837e+11 19008
1 1.5277e+09 1.6931e+11 19013
1 5.3779e+09 1.7316e+11 19036
1 1.2972e+10 1.8075e+11 19079
1 2.2083e+10 1.8986e+11 19128
1 1.5438e+11 3.2216e+11 19658
lm(formula = target ~ cs + type + h + tjg + ggjg + yyjg, data = src)
Coefficients:
(Intercept)
所以最终入选的变量是：formula = target ~ cs + type + h + tjg + ggjg + yyjg
&&相关文章推荐
参考知识库
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：21423次
排名：千里之外
原创：16篇
转载：141篇
(8)(4)(16)(7)(8)(22)(6)(27)(25)(6)(9)(14)(8)(1)(2)