《統(tǒng)計分析軟件：使用R與Python》 課件 第6章-R 語言基本統(tǒng)計推斷

6.1R語言匯總統(tǒng)計量函數(shù)6.2R語言參數(shù)估計方法6.3參數(shù)假設(shè)檢驗6.4非參數(shù)假設(shè)檢驗第6章

R語言基本統(tǒng)計推斷R語言基本統(tǒng)計推斷R語言匯總統(tǒng)計量函數(shù)mean,median,sd,var,quantile,cor,cov非參數(shù)假設(shè)檢驗K-S檢驗Wilcoxon秩和檢驗

卡方檢驗Fisher檢驗R語言參數(shù)估計方法點估計（矩法、MLE）

區(qū)間估計參數(shù)假設(shè)檢驗假設(shè)檢驗的一般步驟、t-檢驗、二項分布、泊松分布、樣本比例檢驗、方差分析6.1R語言匯總統(tǒng)計量函數(shù)1.均值均值(mean)是數(shù)據(jù)的平均數(shù)，描述數(shù)據(jù)取值的平均位置。R語言中用mean()函數(shù)計算樣本的均值，基本語法是：mean(x,trim=0,na.rm=FALSE)其中x是對象(如向量、矩陣、數(shù)組或數(shù)據(jù)框)，trim參數(shù)是計算均值前去掉與均值差較大數(shù)據(jù)的比例，缺省值為0，即包括全部數(shù)據(jù)。當(dāng)na.rm=TRUE時，將剔除缺失數(shù)據(jù)后再求均值。選用參數(shù)trim減少輸入誤差對計算的影響如果第一個學(xué)生的體重少輸入一個點，變?yōu)?50千克

R語言中用median()函數(shù)計算樣本的均值median(x,na.rm=FALSE)median(w)##[1]63.5針對包含異常值的數(shù)據(jù)w1（即第一人體重錯誤輸入為750千克），計算中位數(shù)得：學(xué)生體重數(shù)據(jù)中的中位數(shù)median(w1)##[1]63.5學(xué)生體重中位數(shù)仍為63.5千克，基本不受異常值的影響,具有一定穩(wěn)健性。3.百分位數(shù)

在R軟件中,quantile()函數(shù)計算觀測量的百分位數(shù)quantile(x,probs=seq(0,1,0.25),na.rm=FALSE,names=TRUE,type=7,...)quantile(w)##0%25%50%75%100%##47.4057.8563.5066.7575.00quantile(w,probs=seq(0,1,0.2))##0%20%40%60%80%100%##47.4056.9862.2064.0067.3275.00改變probs參數(shù)的默認(rèn)值，可以給出不同位置的分位數(shù)：4.方差、標(biāo)準(zhǔn)差方差函數(shù)var()和標(biāo)準(zhǔn)差函數(shù)sd()的基本語法為var(x,y=NULL,na.rm=FALSE,use)sd(x,na.rm=FALSE)var(w)sd(w)##[1]56.47##[1]7.5155.極差

在R軟件中,計算極差可借助range()函數(shù)range(...,na.rm=FALSE)ran=range(w)#計算極差ran[2]-ran[1]##[1]27.6對于上述體重數(shù)據(jù),極差計算如下：6.相關(guān)系數(shù)與協(xié)方差對于來自多元總體的數(shù)據(jù)，除了分析各個分量的取值特點外，更重要的是分析各個分量之間的相關(guān)關(guān)系，如相關(guān)系數(shù)和協(xié)方差。R語言中，cor()和cov()函數(shù)可以計算相關(guān)系數(shù)和協(xié)方差，cor(x,y=NULL,use="everything",method=c("pearson","kendall","spearman"))cov(x,y=NULL,use="everything",method=c("pearson","kendall","spearman"))對于JS數(shù)據(jù)集，部分?jǐn)?shù)值變量的協(xié)方差矩陣如下：JS=read.csv("JSdata.csv")cor(JS[,6:8])##heightweightsalary##height1.00000.89480-0.13647##weight0.89481.00000-0.09135##salary-0.1365-0.091351.00000cor()函數(shù)在默認(rèn)情況下得到的結(jié)果是一個方陣（所有變量之間兩兩計算相關(guān)）。cor()函數(shù)也可以計算非方形的相關(guān)矩陣。#身高、體重兩列變量與工資變量之間的相關(guān)系數(shù)cor(JS[,6:7],JS[,8])##[,1]##height-0.13647##weight-0.09135注意：cor()只能計算出相關(guān)系數(shù)，無法給出顯著性水平p值。如果想計算p值，建議使用psych包中的corr.test函數(shù)或Hmisc包里的rcorr()函數(shù)。對于JS數(shù)據(jù)集，部分列的協(xié)方差矩陣如下：JS=read.csv("JSdata.csv")cov(JS[,6:8])##heightweightsalary##height78.8875.37-39.84##weight75.3789.94-28.48##salary-39.84-28.481080.536.2R語言參數(shù)估計方法

估計廢品率估計新生兒的體重估計湖中魚數(shù)估計降雨量

在參數(shù)估計問題中，假定總體分布形式已知，未知的僅僅是一個或幾個參數(shù).參數(shù)估計的兩種方法點估計區(qū)間估計

battery<-c(20.2,22.8,19.5,17.8,20.8)#續(xù)航時間樣本數(shù)據(jù)

6.2.1

點估計

尋求估計量的方法1.矩估計法2.極大似然法3.最小二乘法4.貝葉斯方法……這里我們只介紹前面兩種方法.其基本思想是用樣本矩估計總體矩.理論依據(jù):它是基于一種簡單的“替換”思想建立起來的一種估計方法.是英國統(tǒng)計學(xué)家K.皮爾遜最早提出的.大數(shù)定律矩估計法記總體k階矩為樣本k階矩為用相應(yīng)的樣本矩去估計總體矩的估計方法就稱為矩估計法.記總體k階中心矩為樣本k階中心矩為例：求伽馬分布的矩估計

set.seed(123)x.gam<-rgamma(200,rate=0.5,shape=3.5)##伽馬分布λ=0.5(scale參數(shù))α=3.5(shape參數(shù))#首先算出樣本矩med.gam<-mean(x.gam)##樣本均值var.gam<-var(x.gam)##樣本方差使用矩估計法，把總體矩用樣本矩表示出來lambda_hat<-med.gam/var.gam##lambda的矩估計量alpha_hat<-((med.gam)^2)/var.gam##alpha的矩估計量cat("lambda估計=",lambda_hat,"alpha估計=",alpha_hat,"\n")##lambda估計=0.544alpha估計=3.609優(yōu)點:是簡單易行，并不需要事先知道總體是什么分布缺點:當(dāng)總體類型已知時，沒有充分利用分布提供的信息，且矩估計量不具有唯一性。極大似然法

求極大似然估計量的步驟對數(shù)似然方程例如，對于正態(tài)分布總體，似然函數(shù)如下：

由于似然函數(shù)是正數(shù)的乘積（其中許多可能非常?。?，因此使用對數(shù)似然函數(shù)更方便、更穩(wěn)定

下面，我們生成一組服從正態(tài)分布的數(shù)據(jù)，假設(shè)我們并不知道參數(shù)的真實值，我們用MLE方法對其進(jìn)行估計##[1]4.967可以針對μ的許多候選值繪制似然函數(shù)，并查看曲線最高的位置。使得似然函數(shù)最大的mu值為4.97

引例1估計某明星年齡該明星的年齡35歲—點估計該明星的年齡，九成可能在30-50歲之間6.2.2

區(qū)間估計也就是說，我們希望確定一個區(qū)間，使我們能以比較高的可信度或可靠程度相信它包含真參數(shù)值.這里所說的“可靠程度”是用概率來度量的，稱為置信概率，置信度或置信水平.

求置信區(qū)間的一般步驟1.明確問題,是求什么參數(shù)的置信區(qū)間?置信水平

是多少?

P(a≤S(T,)≤b)=

則就是的100(

)％的置信區(qū)間.可見，確定區(qū)間估計很關(guān)鍵的是要尋找一個待估參數(shù)和估計量T的函數(shù)S(T,),且S(T,)的分布為已知,不依賴于任何未知參數(shù)(這樣我們才能確定一個大概率區(qū)間).而這與總體分布有關(guān)，所以，總體分布的形式是否已知，是怎樣的類型，至關(guān)重要.這里，我們主要討論總體分布為正態(tài)的情形.根據(jù)置信區(qū)間的求解過程自行編寫R函數(shù)利用R中相應(yīng)假設(shè)檢驗函數(shù)給出置信區(qū)間R

base包中沒有直接的函數(shù)做

2已知時置信區(qū)間沒有只知道匯總數(shù)據(jù)的對應(yīng)函數(shù)單正態(tài)總體方差的區(qū)間估計也需自己編寫函數(shù)…R中求置信區(qū)間單正態(tài)總體均值的置信區(qū)間可用t.test()函數(shù)雙正態(tài)總體方差比的置信區(qū)間可用var.test()函數(shù)…假設(shè)方差未知，要估計JS數(shù)據(jù)中身高均值的置信區(qū)間，置信水平為95%，可以使用t.test()函數(shù):JS=read.csv("JSdata.csv")t.test(JS$height)####OneSamplet-test####data:JS$height##t=117,df=38,p-value<2e-16##alternativehypothesis:truemeanisnotequalto0##95percentconfidenceinterval:##163.6169.3##sampleestimates:##meanofx##166.4身高的95%置信區(qū)間為(163.5568,169.3150)當(dāng)然，也可以自己編寫R函數(shù)給出置信區(qū)間。下面我們以單正態(tài)總體均值的置信區(qū)間為例討論如何用R計算置信區(qū)間。

即找到一個隨機(jī)變量的臨界值，使得α/2的面積位于其右側(cè)。這個臨界值可用分位數(shù)函數(shù)qnorm()或qt()求出。

xbar=0.093;sigma=0.012;alpha=0.05;n=60z_half_alpha=qnorm(0.025,mean=0,sd=1,lower.tail=FALSE)cri_value=z_half_alpha*(sigma/sqrt(n))deg_upp=xbar+cri_valuedeg_low=xbar-cri_valuelibrary(scales)c(percent(deg_low,accuracy=0.001),percent(deg_upp,accuracy=0.001))置信區(qū)間為[8.996%,9.604%]例：R計算置信區(qū)間奶茶的容量按標(biāo)準(zhǔn)應(yīng)在標(biāo)稱容量480毫升和550毫升之間.

引例

500毫升奶茶實際370毫升？流水線上奶茶不斷地封裝，然后送到賣場

。

怎么知道這批奶茶的容量是否合格呢？把每一杯都打開倒入量杯,看看容量是否合于標(biāo)準(zhǔn).這樣做顯然不行！6.3假設(shè)檢驗通常的辦法是進(jìn)行抽樣檢查.如每隔1小時，抽查5杯，得5個容量的值X1，…，X5，根據(jù)這些值來判斷生產(chǎn)是否正常.如發(fā)現(xiàn)不正常，就應(yīng)停產(chǎn)，找出原因，排除故障，然后再生產(chǎn)；如沒有問題，就繼續(xù)按規(guī)定時間再抽樣，以此監(jiān)督生產(chǎn)，保證質(zhì)量.

很明顯，不能由5杯容量的數(shù)據(jù)，在把握不大的情況下就判斷生產(chǎn)不正常，因為停產(chǎn)的損失是很大的.

當(dāng)然也不能總認(rèn)為正常，有了問題不能及時發(fā)現(xiàn)，這也要造成損失.

如何處理這兩者的關(guān)系，假設(shè)檢驗面對的就是這種矛盾.

在正常生產(chǎn)條件下，由于種種隨機(jī)因素的影響，每罐奶茶的容量應(yīng)在500毫升上下波動.這些因素中沒有哪一個占有特殊重要的地位.因此，根據(jù)中心極限定理，假定每杯容量服從正態(tài)分布是合理的.它的對立假設(shè)是：稱H0為原假設(shè)（或零假設(shè)，解消假設(shè)）；稱H1為備選假設(shè)（或?qū)α⒓僭O(shè)）.在實際工作中，往往把不輕易否定的命題作為原假設(shè).H0：（=500）H1：

那么，如何判斷原假設(shè)H0

是否成立呢？問題歸結(jié)為根據(jù)小概率原理對差異作定量的分析，以確定其性質(zhì).

“系統(tǒng)誤差”或隨機(jī)誤差?假設(shè)檢驗的一般步驟

32

1檢驗方法Z檢驗t檢驗檢驗常見假設(shè)檢驗的拒絕域R語言中將假設(shè)檢驗存儲在稱為htest

的特殊對象類中。htest對象包含來自假設(shè)檢驗、檢驗統(tǒng)計量（例如，t檢驗的t統(tǒng)計量或相關(guān)性檢驗的相關(guān)系數(shù)）、p值、置信區(qū)間的所有主要結(jié)果?？梢允褂?符號訪問htest對象的特定元素##[1]"statistic""parameter""p.value"##[4]"""estimate""null.value"##[7]"stderr""alternative""method"##[10]""6.3.1t-檢驗t.test(x,y=NULL,alternative=c("two.sided","less","greater"),mu=0,paired=FALSE,var.equal=FALSE,conf.level=0.95,...)若僅出現(xiàn)數(shù)據(jù)x,則進(jìn)行單樣本t檢驗;若出現(xiàn)數(shù)據(jù)x和y,則進(jìn)行雙樣本的t檢驗alternative=c("two.sided","less","greater")用于指定所求置信區(qū)間的類型;alternative="two.sided"是缺省值,表示求置信區(qū)間,alternative="less"表示求置信上限;alternative="greater"表示求置信下限.mu表示均值,它僅在假設(shè)檢驗中起作用,默認(rèn)值為零.單正態(tài)總體參數(shù)t檢驗例：有原始數(shù)據(jù)的t檢驗汽車廠商聲稱其發(fā)動機(jī)排放標(biāo)準(zhǔn)的一個指標(biāo)平均低于20個單位。在抽查了10臺發(fā)動機(jī)之后，得到下面的排放數(shù)據(jù)：17.0、21.7、17.9、22.9、20.7、22.4、17.3、21.8、24.2、25.4。該樣本均值為21.13。究竟能否由此認(rèn)為該指標(biāo)均值超過20？這次我們的假設(shè)檢驗問題就是檢驗統(tǒng)計量為>x=c(17.0,21.7,17.9,22.9,20.7,22.4,17.3,21.8,24.2,25.4)>mu=20>t.test(x,mu=20,alternative="greater")OneSamplet-testdata:xt=1.2336,df=9,p-value=0.1243alternativehypothesis:truemeanisgreaterthan2095percentconfidenceinterval:19.4508Infsampleestimates:meanofx21.13因此，我們沒有證據(jù)否定零假設(shè)(如果顯著性水平小于它)。單正態(tài)總體參數(shù)t檢驗也可以自己寫代碼實現(xiàn)上述功能：>x=c(17.0,21.7,17.9,22.9,20.7,22.4,17.3,21.8,24.2,25.4)>xbar=mean(x);sigma=sd(x)

>mu=20

>n=length(x)

>t<-(xbar-mu)/(sigma/sqrt(n))

>t

[1]1.233576

>p1<-pt(t,df=n-1,lower.tail=FALSE)

>p1

[1]0.1243025健康成年男子脈搏均數(shù)為72次/分。某醫(yī)生在某山區(qū)隨機(jī)抽查健康成年男子25人，其脈搏均數(shù)為74.2次/分，標(biāo)準(zhǔn)差為6.5次/分。根據(jù)這個資料能否認(rèn)為某山區(qū)健康成年男子脈搏數(shù)與一般健康成年男子的不同？檢驗統(tǒng)計量為例：無原始數(shù)據(jù)的t檢驗#根據(jù)公式算出t值x<-74.2;mu<-72;sigma<-6.5;n<-25t<-(x-mu)/(sigma/sqrt(n))#或者用n-1代替n#用pt()函數(shù)，輸入t值和自由度df（n-1），得到p值#雙側(cè)p=2*pt(-abs(t),df=24)檢驗結(jié)果為t=1.692308，顯著性P值=0.1035319>0.05，接受原假設(shè)，說明該法測定的成年男子脈搏數(shù)與總體無顯著差異，認(rèn)為某山區(qū)健康成年男子脈搏數(shù)與一般健康成年男子的相同。思考題

t.test(X,alternative=

“greater”,mu=

225)兩獨立樣本t檢驗兩獨立樣本t檢驗用于檢驗兩組相互獨立數(shù)據(jù)均值間有無顯著差別：先假定兩組相應(yīng)的總體均值相等，看兩組均值實際相差與此假設(shè)是否靠近。

利用t檢驗法檢驗要求:具有相同方差的兩正態(tài)總體均值差的假設(shè).對應(yīng)正態(tài)性檢驗，可以使用shapiro.test()函數(shù)；方差齊性檢驗，則可以使用var.test()函數(shù)或bartlett.test()函數(shù)。JSdata數(shù)據(jù)中性別不同的兩組教師體重均值是否相同？例：有原始數(shù)據(jù)的兩樣本t檢驗第一步：正態(tài)性檢驗兩組數(shù)據(jù)Shapiro-Wilk正態(tài)性檢驗的p值均遠(yuǎn)大于0.05，說明可以認(rèn)為兩組數(shù)據(jù)服從正態(tài)分布。第二步：方差齊性檢驗var.test方差齊性檢驗中的結(jié)果p值大于0.05，說明兩組數(shù)據(jù)基本滿足方差齊性要求。第三步：t檢驗t檢驗中p值大于0.05，不能拒絕原假設(shè)，即不能認(rèn)為兩組教師中體重存在顯著差異。如果數(shù)據(jù)滿足正態(tài)性假設(shè)，但不滿足方差齊性要求，則進(jìn)行修正的WelchT檢驗。var.equal=FALSE例如，測量某兩個地區(qū)水中碳酸鈣的含量，分別從兩個地區(qū)隨機(jī)抽取20份樣品進(jìn)行碳酸鈣檢測，分別得到兩個地區(qū)碳酸鈣含量的均數(shù)（20.95，21.79）和標(biāo)準(zhǔn)差（5.89，3.43）。試判斷兩個地區(qū)水中碳酸鈣的含量是否有差異？例：無原始數(shù)據(jù)的兩樣本t檢驗#輸入對照組實驗組均值x1，x2；組數(shù)n1，n2；方差s1，s2x1<-20.95;x2<-21.79;n1<-20;n2<-20;s1<-5.89;s2<-3.43#計算兩獨立樣本共同的標(biāo)準(zhǔn)差sc<-sqrt((1/n1+1/n2)*((n1-1)*s1**2+(n2-1)*s2**2)/(n1+n2-2))#t值，自由度df，p值t<-(x2-x1)/scdf<-n1+n2-2p<-2*pt(-abs(t),df=df)t;pt=0.5511486，p值0.5847581>0.05，不拒絕原假設(shè)，不能認(rèn)為兩個地區(qū)水中碳酸鈣的含量有差異思考題

t.test(X,Y,var.equal=TRUE,alternative=

“l(fā)ess“)配對t檢驗在醫(yī)學(xué)研究中，常用配對設(shè)計。配對設(shè)計主要有4種情況：同一受試對象處理前后的數(shù)據(jù)，同一受試對象兩個部位的數(shù)據(jù)，同一樣品用兩種方法檢驗的結(jié)果，配對的兩個受試對象分別接受兩種處理后的數(shù)據(jù)。例1.原始數(shù)據(jù)的配對t檢驗判斷簡便法和常規(guī)法測定尿鉛含量的差別有無統(tǒng)計意義，對12份人尿同時用兩種方法進(jìn)行測定，所得結(jié)果如下表所示，請分析兩種測定方法的測量結(jié)果是否不同？配對t檢驗的結(jié)果為：t=0.16232，顯著性p值=0.874>0.05，不能拒絕原假設(shè)H0，說明不能認(rèn)為兩種方法測定尿鉛含量的結(jié)果不相同。#輸入兩組值x<-c(2.41,2.90,2.75,2.23,3.67,4.49,5.16,5.45,2.06,1.64,1.06,0.77)y<-c(2.80,3.04,1.88,3.43,3.81,4.00,4.44,5.41,1.24,1.83,1.45,0.92)#配對樣本t檢驗t.test(x,y,paired=T)慢性支氣管炎病人血中膽堿酯酶活性常常偏高。某校藥理教研室將同性別同年齡的病人與健康人配成8對，測量該值加以比較，配對兩組人差值的均值為0.625，標(biāo)準(zhǔn)差為0.78.問可否通過這一資料得出較明確的結(jié)論？例:無原始數(shù)據(jù)的配對t檢驗配對t檢驗的結(jié)果為：t=2.266368，顯著性p值=0.05778616>0.05，不能拒絕原假設(shè)H0，說明不能認(rèn)為慢性支氣管炎病人與健康人血中膽堿酯酶活性不同。#依次輸入配對樣本的差值d、標(biāo)準(zhǔn)差s、配對數(shù)nd<-0.625;s<-0.78;n<-8#算t值t<-d/(s/sqrt(n))#輸入自由度n-1，pt()函數(shù)得到p值df<-n-1p=2*pt(-abs(t),df=df)t;p

6.3.2二項分布的總體檢驗使用binom.test(x,n,p=0.5,alternative=c("two.sided","less","greater"),conf.level=0.95)，x為指定成功的次數(shù)，n為實驗的總次數(shù)，p為每次實驗成功的概率。

例：二項分布的總體檢驗>binom.test(445,500,p=0.85)Exactbinomialtestdata:445and500numberofsuccesses=445,numberoftrials=500,p-value=0.01207alternativehypothesis:trueprobabilityofsuccessisnotequalto0.8595percentconfidenceinterval:0.85923420.9160509sampleestimates:probabilityofsuccess0.89p-value=0.01207<0.05，故拒絕原假設(shè)，即認(rèn)為種衣劑對種子發(fā)芽率有顯著效果p-value=0.09048>0.05，故不能接受原假設(shè)，即不能認(rèn)為該地區(qū)新生兒染色體低于一般水平例：按照以往經(jīng)驗，新生兒染色體異常率一般為1%，某醫(yī)院觀察了當(dāng)?shù)?00名新生兒，有一例染色體異常，問該地區(qū)新生兒染色體是否低于一般水平？>binom.test(1,400,p=0.01,alternative="less")Exactbinomialtestdata:1and400numberofsuccesses=1,numberoftrials=400,p-value=0.09048alternativehypothesis:trueprobabilityofsuccessislessthan0.0195percentconfidenceinterval:0.00000000.0118043sampleestimates:probabilityofsuccess0.00256.3.3泊松分布的總體檢驗poisson.test(x,T=1,r=1,alternative=c("two.sided","less","greater"),conf.level=0.95)#x-事件數(shù)量,T-事件計數(shù)的時間基準(zhǔn),r-假設(shè)的比率對Poisson總體參數(shù)的檢驗，完全類似于二項分布總體參數(shù)的檢驗。對泊松分布中的速率參數(shù)或兩個速率參數(shù)之間的比率的精確檢驗。

例：泊松分布的總體檢驗

>poisson.test(x=60,T=50,r=1,alternative="greater",conf.level=0.95)ExactPoissontestdata:60timebase:50numberofevents=60,timebase=50,p-value=0.09227alternativehypothesis:trueeventrateisgreaterthan195percentconfidenceinterval:0.9570464Infsampleestimates:eventrate1.2p-value=0.09227>0.05，在0.95的置信區(qū)間下接受原假設(shè)H0，后一小時注冊人數(shù)與前一小時無差異。6.3.4樣本比例的檢驗

prop.test(x,n,p=NULL,alternative="two.sided",correct=TRUE)#x為樣本成功次數(shù)，n是樣本量，p是成功的比例

檢驗統(tǒng)計量

例：樣本比率的檢驗

p-value=0.5477>0.05，無法拒絕原假設(shè)，故認(rèn)為該大學(xué)的男、女學(xué)生家中擁有計算機(jī)的比率一致。>success<-c(23,25)>total<-c(102,135)>prop.test(success,total)2-sampletestforequalityofproportionswithcontinuitycorrectiondata:successoutoftotalX-squared=0.36148,df=1,p-value=0.5477alternativehypothesis:two.sided95percentconfidenceinterval:-0.072564760.15317478sampleestimates:prop1

prop20.2254902

0.18518526.3.5方差分析方差分析（AnalysisofVariance,簡寫為ANOVA）又稱變異數(shù)分析或F檢驗,用于兩個及兩個以上樣本均值差別的顯著性檢驗,目的是推斷兩組或多組數(shù)據(jù)的總體均值是否相同，檢驗兩個或多個樣本均值的差異是否有統(tǒng)計學(xué)意義。方差分析假設(shè)：各樣本須是相互獨立的隨機(jī)樣本各樣本來自正態(tài)分布總體各總體方差相等，即方差齊

aov()函數(shù)的語法為aov(formula,data=dataframe)formula可使用的特殊符號如下，其中y為因變量，A、B、C為自變量。符號用法：～分隔符，左邊為因變量，右邊為自變量。例y～A+B+C+分隔自變量：表示交互項，如y～A+B+A:B*表示所有可能的交互項，如y～A*B*C等價于y～A+B+C+A:B+A:C+B:C+A:B:C^表示交互項達(dá)到的某個次數(shù)，如y～(A+B+C)^2等價于y～A+B+C+A:B+A:C+B:C.表示包含除因變量以外的所有變量。如y～.例：R語言方差分析某醫(yī)院欲研究A、B、C三種降血脂藥物對家兔血清腎素血管緊張素轉(zhuǎn)化酶（ACE）的影響，將家兔隨機(jī)分為三組，均喂以高脂飲食，分別給予不同的降血脂藥物。一定時間后測定家兔血清ACE濃度（u/ml），A組（4544434748444644404542404346474546454344），B組（45484743464748464349464347464746454644454644434245），C組（4748454646444548495049484744454645434445464342），問三組家兔血清ACE濃度是否相同？>a<-c(45,44,43,47,48,44,46,44,40,45,42,40,43,46,47,45,+46,45,43,44)>b<-c(45,48,47,43,46,47,48,46,43,49,46,43,47,46,47,46,+45,46,44,45,46,44,43,42,45)>c<-c(47,48,45,46,46,44,45,48,49,50,49,48,47,44,45,46,+45,43,44,45,46,43,42)>dfCRp<-data.frame(value=c(a,b,c),group=factor(c(rep(“1”,20),rep("2",25),rep("3",23))))方差分析需要一定的假設(shè)，即數(shù)據(jù)集應(yīng)該符合正態(tài)和各組的方差相等，可以分別用shapiro.test和bartlett.test檢驗從P值觀察到這兩個假設(shè)是符合的。對于不符合假設(shè)的情況，我們就要用到非參數(shù)方法，例如Kruskal-Wallis秩和檢驗>shapiro.test(dfCRp$value)#檢查是否正態(tài)分布Shapiro-Wilknormalitytestdata:dfCRp$valueW=0.97397,p-value=0.1654>bartlett.test(value~group,data=dfCRp)Bartletttestofhomogeneityofvariancesdata:valuebygroupBartlett'sK-squared=0.76406,df=2,p-value=0.6825正態(tài)性檢驗和方差齊性檢驗P值均大于0.05，可以認(rèn)為數(shù)據(jù)滿足正態(tài)性和方差齊性的要求。>aovCRp=aov(value~group,data=dfCRp)>summary(aovCRp)

Df

SumSqMeanSq

Fvalue

Pr(>F)

group

2

26.29

13.146

3.2440.0454*Residuals65263.40

4.052---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘

’1p-value=0.0454<0.05，拒絕原假設(shè)，故認(rèn)為三組家兔血清ACE濃度不同。>oneway.test(value~group,data=dfCRp,var.equal=TRUE)#與aov()結(jié)果基本相同。One-wayanalysisofmeansdata:valueandgroupF=3.2442,numdf=2,denomdf=65,p-value=0.0454方差分析只告訴我們這三組之間是不同的，但沒有告訴哪兩組之間有明顯差別，此時需要使用TukeyHSD等函數(shù)進(jìn)行均值的多重比較分析做出推斷。>TukeyHSD(aovCRp)Tukeymultiplecomparisonsofmeans95%family-wiseconfidencelevelFit:aov(formula=value~group,data=dfCRp)$group

diff

lwr

upr

p-adj2-11.1300000

-0.31850529

2.578505

0.15526733-11.5195652

0.04333482

2.995796

0.04224953-20.3895652

-1.00547115

1.784602

0.7817904a組和c組的差異有顯著意義。boxplot(value~group,data=dfCRp,xlab="group",ylab="value")6.4非參數(shù)假設(shè)檢驗參數(shù)檢驗方法推定或假定了總體分布的具體形式（如正態(tài)分布），只是其中含有若干未知參數(shù)。但在許多實際問題中，人們往往對總體的分布知之甚少，很難對總體的分布形式做出正確的假定，最多只能對總體的分布做出諸如連續(xù)型分布、關(guān)于某點對稱分布等一般性的假定。這種不假定總體分布的具體形式，盡量從數(shù)據(jù)本身來獲得所需要的信息的檢驗方法稱為非參數(shù)檢驗方法。非參數(shù)假設(shè)檢驗仍然要求各樣本須是相互獨立的隨機(jī)樣本。非參數(shù)檢驗中常用到的一個概念是秩(rank)，秩指的是樣本數(shù)據(jù)按照升序排列之后，每個觀測值的位置。比如：159183178513719觀測值357891315171819秩12345678910如果有重復(fù)，就取平均，比如1，2，2，3。那么秩就為1，2.5，2.5，4。R語言非參數(shù)檢驗函數(shù)K-S檢驗Wilcoxon秩和檢驗卡方檢驗Fisher檢驗6.4.1K-S檢驗Kolmogorov-Smirnov分布一致性檢驗：檢驗的目的：對于單樣本，檢驗其是否符合某種分布。對于雙樣本，檢驗其是否屬于同一分布。ks檢驗，理論上可以檢驗任何分布.原假設(shè)：X具有分布F(單樣本)或X與Y同分布(雙樣本)ks.text(x,y,…,alternative=c(“two.sided”,“l(fā)ess”,“greater”),exact=NULL,tol=1e-8,

simulate.p.value=FALSE,B=2000)x：數(shù)據(jù)向量y：數(shù)據(jù)向量或用于命名分布函數(shù)的字符串?！河蓎值定義的參數(shù)例：K-S檢驗(單樣本和雙樣本)>x<-c(420,500,920,1380,1510,1650,1760,2100,2300,2350)>ks.test(x,"pexp",1/1500)One-sampleKolmogorov-Smirnovtestdata:xD=0.30148,p-value=0.2654alternativehypothesis:two-sided>set.seed(123)>x=runif(100)>y=runif(100)>ks.test(x,y)Two-sampleKolmogorov-Smirnovtestdata:xandyD=0.09,p-value=0.8127alternativehypothesis:two-sidedP值均大于0.05，可以認(rèn)為數(shù)據(jù)和給定分布一致或兩組數(shù)據(jù)分布無明顯差異。6.4.2Wilcoxon符號秩檢驗、秩和檢驗

例：R語言Wilcoxon秩檢驗

對來自一個總體樣本的檢驗假定某電池廠宣稱該廠生產(chǎn)的某種型號電池壽命的中位數(shù)為140安培小時，為了檢驗該廠生產(chǎn)的電池是否符合其規(guī)定的標(biāo)準(zhǔn)，現(xiàn)從新近生產(chǎn)的一批電池中抽取20個隨機(jī)樣本，并對這個20個電池的壽命進(jìn)行了測試，試用Wilcoxon符號秩檢驗分析該廠生產(chǎn)的電池是否符合其標(biāo)準(zhǔn)。假設(shè)如下：H0:電池中位數(shù)M>=140安培小時H1:電池中位數(shù)M<140安培小時X<-

c(137,140,138.3,139,144.3,139.1,141.7,137.3,133.5,138.2,141.1,139.2,136.5,136.5,135.6,138,140.9,140.6,136.3,134.1)

wilcox.test(X,mu=140,alternative="less",exact=FALSE,correct=FALSE,+=TRUE)Wilcoxonsignedranktestdata:XV=34,p-value=0.007034alternativehypothesis:truelocationislessthan14095percentconfidenceinterval:-Inf139.2sampleestimates:(pseudo)

median138.2p-value=0.007034<0.05，故拒絕原假設(shè)，認(rèn)為中位數(shù)達(dá)不到140今測得10名非鉛作業(yè)工人和7名鉛作業(yè)工人的血鉛值，如下表所示。試用Wilcoxon秩和檢驗分析兩組工人血鉛值有無差異。>nonlead<-c(24,26,29,34,43,58,63,72,87,101)>lead<-c(82,87,97,121,164,208,213)>wilcox.test(nonlead,lead,alternative="less",exact=FALSE,correct=FALSE)Wilcoxonranksumtestdata:nonleadandleadW=4.5,p-value=0.001449alternativehypothesis:truelocationshiftislessthan0P值小于0.05，拒絕原假設(shè)，即鉛作業(yè)工人血鉛值高于非作業(yè)工人。

非成對樣本的秩次和檢驗6.4.4Pearson卡方檢驗“Pearson卡方檢驗”可用于兩種情境的變項比較：擬合優(yōu)度(GoodnessofFit

)檢驗，和獨立性檢驗?！皵M合優(yōu)度”驗證一組觀察值的次數(shù)分配是否異于理論上的分配?！蔼毩⑿詸z驗”驗證從兩個變量抽出的配對觀察值組是否互相獨立（例如：每次都從A國和B國各抽一個人，看他們的反應(yīng)是否與國籍無關(guān)）。

擬合優(yōu)度(總體分布)的卡方檢驗chisq.test(x,y=NULL,correct=TRUE,p=rep(1/length(x),length(x)),rescale.p=FALSE,simulate.p.value=FALSE,B=2000)#p默認(rèn)取均勻分布的概率

因為p=0.0123小于0.05，因此我們不能認(rèn)為這批數(shù)據(jù)服從均勻分布。>y=c(12,5,19,7,7)>chisq.test(y)Chi-squaredtestforgivenprobabilitiesdata:yX-squared=12.8,df=4,p-value=0.0123例：總體分布的卡方檢驗(均勻分布)用Pearson擬合優(yōu)度卡方檢驗檢驗學(xué)生成績是否服從正態(tài)分布:X<-c(25,45,50,54,55,61,64,68,72,75,75,78,79,81,83,84,84,84,85,86,86,86,87,89,89,89,90,91,91,92,100)A<-

table(cut(X,br=c(0,69,79,89,100)))#分組計數(shù)例：總體分布的卡方檢驗(正態(tài)分布)>A