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