|
由于有人問如何使用R進行數(shù)據(jù)正態(tài)性檢驗，所以周老師干脆寫個主題帖解釋一下。如果恰好解決了你的問題，請讀完后給個好評喲~
正態(tài)性檢驗，是很多數(shù)據(jù)分析前要做的準備性工作。例如，你有組數(shù)量性狀的表型值，你想先判斷其是否符合正態(tài)分布，再開展后續(xù)的數(shù)據(jù)分析。
正態(tài)性檢驗，最簡單的方法是使用R語言的shapiro.test命令。如果P value > 5%，則說明數(shù)據(jù)分布近似正態(tài)分布。
當然，你還期望有圖形化的比較，以便在文章中展示。那么有4種畫法。
功能和原理：檢驗樣本的概率分布是否服從某種理論分布。PP概率圖的原理是檢驗實際累積概率分布與理論累積概率分布是否吻合，若吻合，則散點應(yīng)圍繞在一條直線周圍，或者實際概率與理論概率之差分布在對稱于以0為水平軸的帶內(nèi)。QQ概率圖的原理是檢驗實際分位數(shù)與理論分位數(shù)之差分布是否吻合，若吻合，則散點應(yīng)圍繞在一條直線周圍，或者實際分位數(shù)與理論分位數(shù)之差分布在對稱于以0為水平軸的帶內(nèi)。QQ概率圖以樣本的分位數(shù)為橫軸，以指定理論分布的分位數(shù)為縱軸繪制散點圖。
#install.packages('DAAG') library(DAAG) data(possum) attach(possum) # 數(shù)據(jù)準備 fpossum <- possum[possum$sex="='f',]" =""># 只分析這些樣本中的雌性個體 x<- scale(fpossum$totlngth)="" ="">#將totlngth這個表型均一化，即 標準正態(tài)化 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)
圖形表示，數(shù)據(jù)與正態(tài)性略有差異，特別是中部區(qū)域。
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)
看圖直接看和正態(tài)密度函數(shù)的差異度。這張圖在數(shù)量性狀的文章里最常出現(xiàn)。
3經(jīng)驗分布與正態(tài)分布函數(shù)對比
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分析結(jié)果中，觀測的P value和理論p value間的差異，代碼大概如下：
r=read.table('temp_name.txt') # 含位點p值的文件。一行1個位點，p值在第六列；
o=-log10(sort(r$V6,decreasing=F)) # 觀測到的p值，對p value列排序，假設(shè)p value在第六列
e=-log10(ppoints(length(r$V6))) # 生成對應(yīng)的平均分布的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 ')
#結(jié)果如下
4種圖形化展示方式不知道你心水哪種呢？可以都拿走噢~哈哈~明天見~
培訓班詳情>>
|
