|
由于有人问如何使用R进行数据正态性检验，所以周老师干脆写个主题帖解释一下。如果恰好解决了你的问题，请读完后给个好评哟~
正态性检验，是很多数据分析前要做的准备性工作。例如，你有组数量性状的表型值，你想先判断其是否符合正态分布，再开展后续的数据分析。
正态性检验，最简单的方法是使用R语言的shapiro.test命令。如果P value > 5%，则说明数据分布近似正态分布。
当然，你还期望有图形化的比较，以便在文章中展示。那么有4种画法。
功能和原理：检验样本的概率分布是否服从某种理论分布。PP概率图的原理是检验实际累积概率分布与理论累积概率分布是否吻合，若吻合，则散点应围绕在一条直线周围，或者实际概率与理论概率之差分布在对称于以0为水平轴的带内。QQ概率图的原理是检验实际分位数与理论分位数之差分布是否吻合，若吻合，则散点应围绕在一条直线周围，或者实际分位数与理论分位数之差分布在对称于以0为水平轴的带内。QQ概率图以样本的分位数为横轴，以指定理论分布的分位数为纵轴绘制散点图。
#install.packages('DAAG') library(DAAG) data(possum) attach(possum) # 数据准备 fpossum <- possum[possum$sex="='f',]" =""># 只分析这些样本中的雌性个体 x<- scale(fpossum$totlngth)="" ="">#将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)
图形表示，数据与正态性略有差异，特别是中部区域。
library(DAAG) data(possum) attach(possum) fpossum <- possum[possum$sex=""> 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)
看图直接看和正态密度函数的差异度。这张图在数量性状的文章里最常出现。
library(DAAG) data(possum) attach(possum) fpossum <- possum[possum$sex=""> 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)
使用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种图形化展示方式不知道你心水哪种呢？可以都拿走噢~哈哈~明天见~
培训班详情>>
|
