R語(yǔ)言基礎(chǔ)培訓(xùn)第二講常用統(tǒng)計(jì)分析

1、基于基于R R的基本統(tǒng)計(jì)分析的基本統(tǒng)計(jì)分析內(nèi)容提要內(nèi)容提要 描述統(tǒng)計(jì)描述統(tǒng)計(jì) 頻數(shù)表分析頻數(shù)表分析 方差分析方差分析 t檢驗(yàn)檢驗(yàn) 卡方檢驗(yàn)卡方檢驗(yàn) 線性回歸線性回歸 相關(guān)分析相關(guān)分析描述分析描述分析(Descriptive statistics)描述統(tǒng)計(jì)就是把數(shù)據(jù)集所包含的信息加以簡(jiǎn)要地概況，描述統(tǒng)計(jì)就是把數(shù)據(jù)集所包含的信息加以簡(jiǎn)要地概況，如計(jì)算數(shù)據(jù)的數(shù)字特征、制作頻數(shù)表和頻數(shù)圖等等，如計(jì)算數(shù)據(jù)的數(shù)字特征、制作頻數(shù)表和頻數(shù)圖等等，用所獲得的統(tǒng)計(jì)量和圖表來(lái)描述數(shù)據(jù)集所反映的特用所獲得的統(tǒng)計(jì)量和圖表來(lái)描述數(shù)據(jù)集所反映的特征和規(guī)律，使得研究的問(wèn)題更加簡(jiǎn)單、直觀。征和規(guī)律，使得研究的問(wèn)題更加簡(jiǎn)單、直觀

2、。 描述性統(tǒng)計(jì)主要包括反映數(shù)據(jù)集中趨勢(shì)的特征值描述性統(tǒng)計(jì)主要包括反映數(shù)據(jù)集中趨勢(shì)的特征值( (比如比如平均數(shù)、中位數(shù)、眾數(shù)、分位數(shù)平均數(shù)、中位數(shù)、眾數(shù)、分位數(shù)) )、數(shù)據(jù)離散程度的、數(shù)據(jù)離散程度的特征值特征值( (比如方差、標(biāo)準(zhǔn)差、值域、變異系數(shù)比如方差、標(biāo)準(zhǔn)差、值域、變異系數(shù)) )和數(shù)和數(shù)據(jù)分布形態(tài)的特征值據(jù)分布形態(tài)的特征值( (比如偏度、峰度比如偏度、峰度) )。 標(biāo)準(zhǔn)差（std.dev）和標(biāo)準(zhǔn)誤（SE.mean）標(biāo)準(zhǔn)差（標(biāo)準(zhǔn)差（std.dev）真實(shí)均值真實(shí)均值SE樣本均值樣本均值當(dāng)樣本含量 n 足夠大時(shí)，標(biāo)準(zhǔn)差趨向穩(wěn)定；而標(biāo)準(zhǔn)誤隨n的增大而減小，甚至趨于0 。 標(biāo)準(zhǔn)差（std.dev）和

3、標(biāo)準(zhǔn)誤（SE.mean） 比如，某學(xué)校共有500名學(xué)生，現(xiàn)在要通過(guò)抽取樣本量為30的一個(gè)樣本，來(lái)推斷學(xué)生的身高。這時(shí)可以依據(jù)抽取的樣本信息，計(jì)算出樣本的均值與標(biāo)準(zhǔn)差。如果我們抽取的不是一個(gè)樣本，而是10個(gè)樣本，每個(gè)樣本30人，那么每個(gè)樣本都可以計(jì)算出均值，這樣就會(huì)有10個(gè)均值。也就是形成了一個(gè)10個(gè)數(shù)字的數(shù)列，然后計(jì)算這10個(gè)數(shù)字的標(biāo)準(zhǔn)差，此時(shí)的標(biāo)準(zhǔn)差就是標(biāo)準(zhǔn)誤。但是，在實(shí)際抽樣中我們不可能抽取10個(gè)樣本。所以，標(biāo)準(zhǔn)誤就由樣本標(biāo)準(zhǔn)差除以樣本量來(lái)表示。當(dāng)然，這樣的結(jié)論也不是隨心所欲，而是經(jīng)過(guò)了統(tǒng)計(jì)學(xué)家的嚴(yán)密證明的。 SE.mean= std.dev / n1/2在實(shí)際的應(yīng)用中，標(biāo)準(zhǔn)差主要有兩點(diǎn)作

4、用，一是統(tǒng)計(jì)量樣本離散程度的表征；二是用來(lái)對(duì)樣本進(jìn)行標(biāo)準(zhǔn)化處理，即樣本觀察值減去樣本均值，然后除以標(biāo)準(zhǔn)差，這樣就變成了標(biāo)準(zhǔn)正態(tài)分布。標(biāo)準(zhǔn)誤的作用主要是用來(lái)做區(qū)間估計(jì)，常用的估計(jì)區(qū)間是均值加減n倍的標(biāo)準(zhǔn)誤（例如95%的置信區(qū)間是：均值+1.96*SE）標(biāo)準(zhǔn)差（std.dev）和標(biāo)準(zhǔn)誤（SE.mean）95% CI:假設(shè)上面這個(gè)隨機(jī)抽樣估計(jì)學(xué)生身高的例子，抽樣100次，每次抽10個(gè)學(xué)生測(cè)量身高，均值估計(jì)值及標(biāo)準(zhǔn)誤為 152cm12cm。但有時(shí)需要表示為估計(jì)量的95%的置信區(qū)間152cm-1.96*12cm，152cm+1.96*12cm。可以解釋為，如果從再?gòu)目傮w中抽樣100次（每次抽樣10個(gè)），

5、產(chǎn)生100個(gè)平均值，這100個(gè)平均值將有95次落在152cm-1.96*12cm，152cm+1.96*12cm這個(gè)范圍內(nèi)，5次落在這個(gè)范圍外，如果抽樣次數(shù)越多，這個(gè)推斷越準(zhǔn)確。這個(gè)來(lái)源于中心極限定理的應(yīng)用：任何分布(總體)抽樣n次，每次抽樣的和符合正態(tài)分布。通俗一點(diǎn)說(shuō), 不管是學(xué)校的學(xué)生身高是怎么分布, 每次隨機(jī)抽取10個(gè)求和, 抽取n次，這n個(gè)身高總和是符合正態(tài)分布的。平均身高為身高總和除于10，所以平均身高也是正態(tài)分布的。正態(tài)分布雙尾95%的分界點(diǎn)所對(duì)應(yīng)的值剛好是1.96。Histogram of observed datayieldDensity40608010012014016018

6、00.0000.0050.0100.015峰度（Kurtosis）峰度（Kurtosis）是描述某變量所有取值分布形態(tài)陡緩程度的統(tǒng)計(jì)量。 它是和正態(tài)分布相比較的。 Kurtosis=0 與正態(tài)分布的陡緩程度相同。 Kurtosis0 比正態(tài)分布的高峰更加陡峭尖頂 Kurtosis0 正偏差數(shù)值較大，為正偏或右偏。長(zhǎng)尾巴拖在右邊 Skewness0 負(fù)偏差數(shù)值較大，為負(fù)偏或左偏。長(zhǎng)尾巴拖在左邊計(jì)算公式：Skewness 越大，分布形態(tài)偏移程度越大頻數(shù)表頻數(shù)表(Frequency table)分析分析 頻數(shù)表分析是對(duì)數(shù)據(jù)集按數(shù)據(jù)范圍分成若干區(qū)間，頻數(shù)表分析是對(duì)數(shù)據(jù)集按數(shù)據(jù)范圍分成若干區(qū)間，即分成若

7、干組，求出每組組中值，各組數(shù)據(jù)用組中值代即分成若干組，求出每組組中值，各組數(shù)據(jù)用組中值代替，計(jì)算各組數(shù)據(jù)的頻數(shù)，并作出頻數(shù)表。替，計(jì)算各組數(shù)據(jù)的頻數(shù)，并作出頻數(shù)表。 頻數(shù)表分析例子頻數(shù)表分析例子summary(oats$yield)# 計(jì)算頻數(shù)計(jì)算頻數(shù)A - table(cut(oats$yield, breaks = 40 + 20 * (0:7) round(prop.table(A) * 100,2) # 計(jì)算頻數(shù)比例計(jì)算頻數(shù)比例# 畫(huà)頻數(shù)表畫(huà)頻數(shù)表hist(oats$yield, # breaks = 7, xlim = c(40,180), xlab = yield, main =

8、Frequency chart of yield)方差分析方差分析ANOVA方差分析是一種在若干組能相互比較的試驗(yàn)數(shù)據(jù)中，把方差分析是一種在若干組能相互比較的試驗(yàn)數(shù)據(jù)中，把產(chǎn)生變異的原因加以區(qū)分的方法與技術(shù)，其主要用途是產(chǎn)生變異的原因加以區(qū)分的方法與技術(shù)，其主要用途是研究外界因素或試驗(yàn)條件的改變對(duì)試驗(yàn)結(jié)果影響是否顯研究外界因素或試驗(yàn)條件的改變對(duì)試驗(yàn)結(jié)果影響是否顯著。著。類型：?jiǎn)我蛩胤讲罘治鲱愋停簡(jiǎn)我蛩胤讲罘治?One-way ANOVA)、雙因素方差、雙因素方差分析分析(Two-way ANOVA) 。方差分析的基本模型是線性模型，并假設(shè)隨機(jī)變量是獨(dú)方差分析的基本模型是線性模型，并假設(shè)隨機(jī)變量

9、是獨(dú)立、正態(tài)和等方差的。立、正態(tài)和等方差的。方差分析是根據(jù)平方和的加和原理，利用方差分析是根據(jù)平方和的加和原理，利用 F 檢驗(yàn)，進(jìn)而檢驗(yàn)，進(jìn)而判斷試驗(yàn)因素對(duì)試驗(yàn)結(jié)果的影響是否顯著。判斷試驗(yàn)因素對(duì)試驗(yàn)結(jié)果的影響是否顯著。單因素方差分析單因素方差分析# Tukey HSD 方法方法 #install.packages(multcomp)library(multcomp)tuk - glht(fit, linfct = mcp(Treat = Tukey)summary(tuk) # standard displaytuk.cld - cld(tuk) # letter-based displayo

10、par chisq.test(freq,p=probs) Chi-squared test for given probabilitiesdata: freqX-squared = 6.7, df = 5, p-value = 0.2423卡方檢驗(yàn)28x = c(100,110,80,55,14)probs = c(29, 21, 17, 17, 16)/100 chisq.test(x,p=probs) chisq.test(x,p=probs) Chi-squared test for given probabilitiesdata: xX-squared = 55, df = 4, p-

11、value = 2.685e-11卡方檢驗(yàn)（列聯(lián)表）29yesbelt = c(12813,647,359,42)nobelt = c(65963,4000,2642,303)chisq.test(data.frame(yesbelt,nobelt) chisq.test(data.frame(yesbelt,nobelt) Pearsons Chi-squared testdata: data.frame(yesbelt, nobelt)X-squared = 59, df = 3, p-value = 8.61e-13練習(xí)四練習(xí)四以數(shù)據(jù)以數(shù)據(jù)stu.data.csv為例，為例，試對(duì)體重做頻

12、數(shù)分析。試對(duì)體重做頻數(shù)分析。請(qǐng)分析身高是否符合正態(tài)分布？請(qǐng)分析身高是否符合正態(tài)分布？試分析性別對(duì)體重有無(wú)影響。試分析性別對(duì)體重有無(wú)影響。問(wèn)題問(wèn)題4：請(qǐng)檢驗(yàn)總體平均體重與：請(qǐng)檢驗(yàn)總體平均體重與60kg有無(wú)顯著差有無(wú)顯著差異？男生和女生的平均體重有無(wú)顯著差異？異？男生和女生的平均體重有無(wú)顯著差異？問(wèn)題問(wèn)題5：男女生比例是否符合：男女生比例是否符合 1.2 : 1.0？練習(xí)四練習(xí)四 答案答案df-read.csv(file=stu.data.csv,header=T)#問(wèn)題問(wèn)題1A - table(cut(df$weight, breaks = 40 + 15 * (0:7) round(prop.

13、table(A) * 100,2) # 計(jì)算頻數(shù)比例計(jì)算頻數(shù)比例hist(df$weight, breaks = 7, xlim = c(40,140), xlab = weight, main = Frequency chart of weight)#問(wèn)題問(wèn)題2shapiro.test(df$height)#問(wèn)題問(wèn)題3fit-aov(weight Sex,data=df)summary(fit)library(agricolae)duncan.test(fit, Sex, alpha=0.05)$groups#問(wèn)題問(wèn)題4t.test(df$weight, mu = 60, alternati

14、ve = two.sided)wt.m-subset(df$weight,df$Sex=男男)wt.f-subset(df$weight,df$Sex=女女)var.test(wt.m,wt.f) #等方差檢驗(yàn)等方差檢驗(yàn)t.test(wt.m,wt.f, paired=F) #問(wèn)題問(wèn)題5summary(df$Sex)ct-c(87, 33)pt 0）,公式是有效的。公式是有效的。多項(xiàng)式回歸多項(xiàng)式回歸fit - lm( weight N + I(N2), data = df ) # 二次項(xiàng)回歸模型二次項(xiàng)回歸模型 summary( fit ) # 回歸分析結(jié)果回歸分析結(jié)果 plot( df$N,

15、df$weight ) # 繪制散點(diǎn)圖繪制散點(diǎn)圖 lines( df$N, fitted(fit) ) # 添加回歸線添加回歸線 添加光滑曲線添加光滑曲線require(graphics)plot(cars, main = lowess(cars)lines(cars)lines(lowess(cars), col = 2)多元線性逐步回歸多元線性逐步回歸最優(yōu)線性回歸方程為：最優(yōu)線性回歸方程為： y= -649.779 + 14.592y= -649.779 + 14.592x1 + 6.841x1 + 6.841x2 + 9.329x2 + 9.329x3 x3 回歸方程表明：對(duì)于川農(nóng)回歸方

16、程表明：對(duì)于川農(nóng) 16 16 號(hào)小麥而言，當(dāng)號(hào)小麥而言，當(dāng) x2 x2 和和 x3 x3 固定時(shí)，穗數(shù)固定時(shí)，穗數(shù) x1 x1 每增加每增加 1 1萬(wàn)萬(wàn)/ /畝，產(chǎn)量畝，產(chǎn)量 y y 將平均增加將平均增加 14.592Kg/14.592Kg/畝；當(dāng)畝；當(dāng) x1 x1 和和 x3 x3 固定時(shí)，每穗粒數(shù)固定時(shí)，每穗粒數(shù) x2 x2 每增加每增加 1 1粒，粒，產(chǎn)量產(chǎn)量 y y 將平均增加將平均增加 6.841Kg/6.841Kg/畝；當(dāng)畝；當(dāng) x1 x1 和和 x2 x2 固定時(shí)，千粒固定時(shí)，千粒重重 x3 x3 每增加每增加 1g1g，產(chǎn)量，產(chǎn)量y y 將平均增加將平均增加 9.329 Kg/

17、9.329 Kg/畝。畝。自變量的作用主次自變量的作用主次第一種方法：計(jì)算通徑系數(shù)第一種方法：計(jì)算通徑系數(shù)(path coefficient, p)(path coefficient, p)。df.2 - as.data.frame( scale(df) ) #df.2 - as.data.frame( scale(df) ) #對(duì)數(shù)據(jù)集做標(biāo)準(zhǔn)化處理對(duì)數(shù)據(jù)集做標(biāo)準(zhǔn)化處理 lmfit - lm(y x1 + x2 + x3, data = df.2) lmfit round(coef(lmfit),3) round(coef(lmfit),3) (Intercept) x1 x2 x3 (Int

18、ercept) x1 x2 x3 0.000 0.777 0.410 0.609 0.000 0.777 0.410 0.609 第二種方法：根據(jù)第二種方法：根據(jù) x1x1、x2x2、x3 x3 回歸系數(shù)的回歸系數(shù)的 t t 值大小，也可判斷自變值大小，也可判斷自變量量 x x 的主次順序，凡是的主次順序，凡是 t t 值較大者就是較重要的因子。值較大者就是較重要的因子。結(jié)論：結(jié)論：自變量自變量 x x 的主次順序?yàn)椋核霐?shù)的主次順序?yàn)椋核霐?shù) x1 x1 千粒重千粒重 x3 x3 每穗粒數(shù)每穗粒數(shù) x2x2。交互作用多元線性回歸交互作用多元線性回歸練習(xí)四練習(xí)四以數(shù)據(jù)以數(shù)據(jù)stu.data.csv

19、為例，為例，試對(duì)體重做頻數(shù)分析。試對(duì)體重做頻數(shù)分析。請(qǐng)分析身高是否符合正態(tài)分布？請(qǐng)分析身高是否符合正態(tài)分布？試分析性別對(duì)體重有無(wú)影響。試分析性別對(duì)體重有無(wú)影響。問(wèn)題問(wèn)題4：請(qǐng)檢驗(yàn)總體平均體重與：請(qǐng)檢驗(yàn)總體平均體重與60kg有無(wú)顯著差有無(wú)顯著差異？男生和女生的平均體重有無(wú)顯著差異？異？男生和女生的平均體重有無(wú)顯著差異？問(wèn)題問(wèn)題5：男女生比例是否符合：男女生比例是否符合 1.2 : 1.0？練習(xí)四練習(xí)四 答案答案df-read.csv(file=stu.data.csv,header=T)#問(wèn)題問(wèn)題1A - table(cut(df$weight, breaks = 40 + 15 * (0:7)

20、 round(prop.table(A) * 100,2) # 計(jì)算頻數(shù)比例計(jì)算頻數(shù)比例hist(df$weight, breaks = 7, xlim = c(40,140), xlab = weight, main = Frequency chart of weight)#問(wèn)題問(wèn)題2shapiro.test(df$height)#問(wèn)題問(wèn)題3fit-aov(weight Sex,data=df)summary(fit)library(agricolae)duncan.test(fit, Sex, alpha=0.05)$groups#問(wèn)題問(wèn)題4t.test(df$weight, mu = 6

21、0, alternative = two.sided)wt.m-subset(df$weight,df$Sex2=1)wt.f-subset(df$weight,df$Sex2=2)var.test(wt.m,wt.f) #等方差檢驗(yàn)等方差檢驗(yàn)t.test(wt.m,wt.f, paired=F) #問(wèn)題問(wèn)題5summary(df$Sex)ct-c(87, 33)pt corr.test(df, use = complete)Call:corr.test(x = df, use = complete)Correlation matrix h dbh v cpro wd wpro tl tw l

22、rt h dbh v cpro wd wpro tl tw lrth 1.00 0.86 0.90 -0.43 -0.29 0.32 -0.20 -0.04 -0.17h 1.00 0.86 0.90 -0.43 -0.29 0.32 -0.20 -0.04 -0.17dbh 0.86 1.00 0.98 -0.41 -0.35 0.39 -0.18 -0.15 -0.07dbh 0.86 1.00 0.98 -0.41 -0.35 0.39 -0.18 -0.15 -0.07v 0.90 0.98 1.00 -0.43 -0.38 0.43 -0.20 -0.11 -0.12v 0.90 0

23、.98 1.00 -0.43 -0.38 0.43 -0.20 -0.11 -0.12cpro -0.43 -0.41 -0.43 1.00 0.22 -0.25 0.18 -0.08 0.20cpro -0.43 -0.41 -0.43 1.00 0.22 -0.25 0.18 -0.08 0.20wd -0.29 -0.35 -0.38 0.22 1.00 -0.98 -0.09 -0.08 0.00wd -0.29 -0.35 -0.38 0.22 1.00 -0.98 -0.09 -0.08 0.00wpro 0.32 0.39 0.43 -0.25 -0.98 1.00 0.02 -

24、0.02 0.01wpro 0.32 0.39 0.43 -0.25 -0.98 1.00 0.02 -0.02 0.01tl -0.20 -0.18 -0.20 0.18 -0.09 0.02 1.00 0.16 0.79tl -0.20 -0.18 -0.20 0.18 -0.09 0.02 1.00 0.16 0.79tw -0.04 -0.15 -0.11 -0.08 -0.08 -0.02 0.16 1.00 -0.47tw -0.04 -0.15 -0.11 -0.08 -0.08 -0.02 0.16 1.00 -0.47lrt -0.17 -0.07 -0.12 0.20 0.

25、00 0.01 0.79 -0.47 1.00lrt -0.17 -0.07 -0.12 0.20 0.00 0.01 0.79 -0.47 1.00Sample Size Sample Size 1 301 30Probability values (Entries above the diagonal are adjusted for multiple tests.) h dbh v cpro wd wpro tl tw lrt h dbh v cpro wd wpro tl tw lrth 0.00 0.00 0.00 0.53 1.00 1.00 1.00 1.00 1.00h 0.0

26、0 0.00 0.00 0.53 1.00 1.00 1.00 1.00 1.00dbh 0.00 0.00 0.00 0.63 1.00 0.85 1.00 1.00 1.00dbh 0.00 0.00 0.00 0.63 1.00 0.85 1.00 1.00 1.00v 0.00 0.00 0.00 0.53 0.92 0.53 1.00 1.00 1.00v 0.00 0.00 0.00 0.53 0.92 0.53 1.00 1.00 1.00cpro 0.02 0.02 0.02 0.00 1.00 1.00 1.00 1.00 1.00cpro 0.02 0.02 0.02 0.

27、00 1.00 1.00 1.00 1.00 1.00wd 0.12 0.06 0.04 0.24 0.00 0.00 1.00 1.00 1.00wd 0.12 0.06 0.04 0.24 0.00 0.00 1.00 1.00 1.00wpro 0.08 0.03 0.02 0.18 0.00 0.00 1.00 1.00 1.00wpro 0.08 0.03 0.02 0.18 0.00 0.00 1.00 1.00 1.00tl 0.28 0.35 0.29 0.35 0.65 0.92 0.00 1.00 0.00tl 0.28 0.35 0.29 0.35 0.65 0.92 0.00 1.00 0.00tw 0