r語言卡方檢驗和似然比檢驗_r語言似然比檢驗代碼 - CSDN

本文對應《R語言實戰》第9章：方差分析；第10章：功效分析

 ====================================================================

方差分析：

回歸分析是通過量化的預測變量來預測量化的響應變量，而解釋變量裡含有名義型或有序型因子變量時，我們關注的重點通常會從預測轉向組別差異的分析，這種分析方法就是方差分析（ANOVA）。因變量不只一個時，稱為多元方差分析（MANOVA）。有協變量時，稱為協方差分析（ANCOVA）或多元協方差分析（MANCOVA）。

#基本格式aov(formula, data = dataframe)

　　

基本表達式符號參考回歸中的表格

研究設計的表達式

下表中，小寫字母表示定量變量，大寫字母表示組別因子，Subject是對被試者獨有的標識變量

設計	表達式
單因素ANOVA	y ~ A
含單個協變量的單因素ANCOVA	y ~ x + A
雙因素ANOVA	y ~ A * B
含兩個協變量的ANCOVA	y ~ x1 + x2 + A*B
隨機化區組	y ~ B + A （B是區組因子）
單因素組內ANOVA	y ~ A + Error(Subject/A)
含單個組內因子（W）和單個組間因子（B）的重複測量ANOVA	y ~ B * W + Error(Subject/W)

 

表達式中的各項順序：

有兩種情況會造成影響：（1）因子不止一個，並且是非平衡設計；（2）存在協變量。出現任意一種情況時，等式右邊的變量都與其他每個變量相關，此時我們無法清晰地劃分它們對因變量的影響。

例如，對於雙因素方差分析，若不同處理方式中的觀測數不同，那麼模型y ~ A * B與模型y ~ B * A的結果不同

R默認類型1（序貫型）方法計算ANOVA效應。第一個模型可以這樣寫：y ~ A + B + A : B

R中的ANOVA表的結果將評價：

1. A對y的影響
2. 控制A時，B對y的影響
3. 控制A和B的主效應時，A與B的交互效應

順序很重要

當自變量與其他自變量或者協變量相關時，沒有明確的方法可以評價自變量對因變量的貢獻。例如，含因子A、B和因變量y的雙因素不平衡因子設計，有三種效應：A和B的主效應，A和B的交互效應。假設你正使用如下表達式對數據進行建模：

Y ~ A + B + A : B

有三種類型的方法可以分解等式右邊各效應對y所解釋的方差。

類型1（序貫型）

效應根據表達式中先出現的效應進行調整。A不做調整，B根據A調整，A:B交互項根據A和B調整。

類型2（分層型）

效應根據同水平或低水平的效應做調整。A根據B調整，B依據A調整，A:B交互項同時根據A和B調整

類型3（邊界型）

每個效應根據模型其他各效應做相應調整。A根據B和A:B做調整，A:B交互項根據A和B調整。

 

R默認調用類型1方法，其他軟體（比如SAS或SPSS）默認調用類型3方法。

 

樣本大小越不平衡，效應項的順序對結果的影響越大。一般來說，越基礎性的效應越需要放在表達式前面。

具體來講，首先是協變量，然後是主效應，接著是雙因素的交互項，再接著是三因素的交互項，以此類推。

對於主效應，越基礎性的變量越應放在表達式前面，比如性別要放在處理方式之前。有一個基本準則：若研究設計不是正交的（即因子與協變量相關），一定要謹慎設置效應的順序。

需要注意的是，car包中的Anova()函數提供了類型2和類型3方法的選項，而aov()函數使用的是類型1方法。若想使結果與其他軟體提供的結果保持一致，可以使用Anova()函數。

 

 

單因素方差分析：

多重比較：

#各組均值差異的成對檢驗TukeyHSD(fit)#glht()函數提供了更為全面的方法，適用於線性模型和廣義線性模型library(multcomp)tuk <- glht(fit, linfct = mcp(trt = 「Tukey」))plot(cld(tuk, level = 0.05), col = 「lightgrey」)

　　

評估檢驗的假設條件： 

單因素方差分析的假設條件：因變量服從正態分布；各組方差相等。

#Q-Q圖檢驗正態性假設library(car)qqPlot(lm(response ~ trt, data = cholesterol), simulate = TRUE, labels = FALSE)#數據落在95%置信區間範圍內，說明滿足正態性假設#方差齊性檢驗（Bartlett檢驗）bartlett.test(response ~ trt, data = cholesterol)#p大於0.05說明沒有顯著不同 #方差齊性分析對離群點非常敏感，需要補做一次離群點的檢測library(car)outlierTest(fit)#若沒有離群點，說明上述結果較為可信

　　

單因素協方差分析： 

#單因素ANCOVAdata(litter, package = 「multcomp」)attach(litter)table(dose)aggregate(weight, by = list(dose), FUN = mean)fit <- aov(weight ~ gesttime + dose)summary(fit) #獲取調整的組均值（即去除協變量效應後的組均值）library(effects)effect(「dose」, fit)

　　

評估檢驗的假設條件： 

單因素協方差分析的假設條件：正態性、同方差性。

與單因素方差分析的假設條件判斷相同。

結果可視化：

library(HH)ancova(weight ~ gesttime + dose, data = litter)

　　

雙因素方差分析：

#table()函數觀察是否是均衡設計attach(ToothGrowth)table(supp, dose)aggregate(len, by = list(supp, dose), FUN = mean)aggregate(len, by = list(supp, dose), FUN = sd)fit <- aov(len ~ supp*dose)summary(fit) #可視化處理interaction.plot(dose, supp, len, type = 「b」,col = c(「red」, 「blue」), pch = c(16, 18),main = 「Interaction between Dose and Supplement Type」)#或者是library(gplots)plotmeans(len ~ interaction(supp, dose, sep = 「 」),connect = list(c(1, 3, 5), c(2, 4, 6)),col = c(「red」, 「darkgreen」),main = 「Interaction Plot with 95% CIs」,xlab = 「Treatment and Dose Combination」)#也可以（推薦）library(HH)interaction2wt(len ~ supp * dose)

　　

重複測量方差分析：受試者被測量不止一次，重點關注含一個組內和一個組間因子的重複測量方差分析。 

混合模型設計概述：

由於傳統的重複測量方差分析假設任意組內因子的協方差矩陣為球形，並且任意組內因子兩水平間的方差之差都相等。但在現實中這種假設不可能滿足，於是衍生了一系列備選方法：

使用lme4包中的lmer()函數擬合線性混合模型；

使用car包中的Anova()函數調整傳統檢驗統計量以彌補球形假設的不滿足（例如Geisser-Greenhouse校正）；

使用nlme包中的gls()函數擬合給定方差-協方差結構的廣義最小二乘模型；

用多元方差分析對重複測量數據進行建模。

 

 

多元方差分析：

當因變量（結果變量）不止一個時，可用多元方差分析（MANOVA）對它們同時進行分析。

評估假設檢驗：單因素多元方差分析的假設前提，一個是多元正態性，一個是方差-協方差矩陣同質性。

多元正態性假設即指因變量組合成的向量服從一個多元正態分布，可以用Q-Q圖來檢驗該假設條件。方差-協方差矩陣同質性即指各組的協方差矩陣相同，通常可用Box’s M檢驗來評估該假設。

最後可使用mvoutlier包中的aq.plot()函數來檢驗多元離群點。

library(mvoutlier)outliers <- aq.plot(y)ouliers

　　

穩健多元方差分析： 

如果多元正態性或者方差-協方差均值假設都不滿足，又或者擔心多元離群點，那麼可以考慮用穩健或非參數版本的MANOVA檢驗。

穩健單因素MANOVA可通過rrcov包中的Wilks.test()函數實現。vegan包中的adonis()函數則提供了非參數MANOVA的等同形式。

#Wilks.test()函數應用示例library(rrcov)Wilks.test(y, shelf, method = 「mcd」)

　　

用回歸做ANOVA： 

事實上，ANOVA和回歸都是廣義線性模型的特例。因此可以用lm()函數來分析。

這部分看不大懂，以後再回頭看吧。

 

=========================================================================

功效分析：

功效分析可以幫助在給定置信度的情況下，判斷檢測到給定效應值時所需的樣本量。反過來，它也可以幫助你在給定置信度水平情況下，計算在某樣本量內能檢測到給定效應值的概率。如果概率低得難以接受，修改或者放棄這個實驗將是一個明智的選擇。

本章中將學習如何對多種統計檢驗進行功效分析，包括比例檢驗、t檢驗、卡方檢驗、平衡的單因素ANOVA、相關性分析，以及線性模型分析。由於功效分析針對的是假設檢驗，我們將首先簡單回顧零假設顯著性檢驗（NHST）過程，然後學習如何用R進行功效分析，主要關注於pwr包。

 

假設檢驗回顧：

首先對總體分布參數作零假設H0, 從總體分布中抽樣，通過樣本計算統計量對總體參數進行推斷。假定H0為真，如果計算獲得的觀測樣本的統計量的概率非常小，便可以拒絕原假設，接受備擇假設H1。

Ⅰ型錯誤：H0為真卻拒絕H0

Ⅱ型錯誤：H0為假卻不拒絕H0

 

研究過程中，研究者通常關注四個量：

樣本大小：實驗設計中每種條件/組中觀測的數目

顯著性水平（即alpha，也就是概率p的閾值）：Ⅰ型錯誤的概率

功效：（1-P(Ⅱ型錯誤)）可以看作是真是效應發生的概率

效應值：指的是在備擇或研究假設下效應的量。效應值的表達式依賴於假設檢驗中使用的統計方法

 

這四個量緊密相關，給定其中任意三個量，便可推算第四個量。本章就是利用這一點進行各種各樣的功效分析。

 

pwr包中的函數

函數	功效計算的對象
pwr.2p.test()	兩比例（n相等）
pwr.2p2n.test()	兩比例（n不相等）
pwr.anova.test()	平衡的單因素ANOVA
pwr.chisq.test()	卡方檢驗
pwr.f2.test()	廣義線性模型
pwr.p.test()	比例（單樣本）
pwr.r.test()	相關係數
pwr.t.test()	t檢驗（單樣本、兩樣本、配對）
pwr.t2n.test()	t檢驗（n不相等的兩樣本）

 

t檢驗

pwr.t.test(n = , d = , sig.level = , power = , alternative = )

　　

n為樣本大小 

 

sig.level表示顯著性水平，默認0.05

power為功效水平

type指檢驗類型：雙樣本t檢驗（two.sample）、單樣本檢驗（one.sample）或相依樣本t檢驗（paired）默認為雙樣本t檢驗。

alternative指統計檢驗是雙側檢驗（two.sided）還是單側檢驗（less或greater）。默認雙側

 

 

方差分析

pwr.anova.test(k = , n = , f = , sig.level = , power = )

 

 

 

相關性

pwr.r.test(n = , r = , sig.level = , power = , alternative = )

　　

其中，n是觀測數目，r是效應值（通過線性相關係數衡量），sig.level是顯著性水平，power是功效水平，alternative指定顯著性檢驗是雙邊檢驗（two.sided）還是單邊檢驗（less或greater） 

 

 

線性模型

pwr.f2.test(u = , v = , f2 = , sig.level = , power = )

　　

 

當要評價一組預測變量對結果的影響程度時，適宜用第一個公式計算f2；當要評價一組預測變量對結果的影響超過第二組變量（協變量）多少時，適宜用第二個公式。

 

 

比例檢驗

pwr.2p.test(h = , n = , sig.level = , power = )

　　

其中，h是效應值，n是各組相同的樣本量。效應值h定義如下： 

 

可用ES.h(p1, p2)函數進行計算

當各組n不相同時：

pwr.2p2n.test(h = , n1 = , n2 = , sig.level = , power = )

同樣，alternative可以設定單邊檢驗或雙邊檢驗（默認）

 

卡方檢驗

卡方檢驗常常用來評價兩個類別型變量的關係：典型的零假設是變量之間獨立，備擇假設是不獨立。

pwr.chisq.test(w = , N = , df = , sig.level = , power = )

其中，w是效應值，N是總樣本大小，df是自由度。 

 

此處從1到m進行求和，連加號上的m指的是列聯表中單元格的數目。函數ES.w2(P)可以計算雙因素列聯表中備擇假設的效應值，P是一個假設的雙因素概率表。

 

 

新情況中選擇合適的效應值：

功效分析中，預期效應值是最難決定的參數。通常需要你對主題有一定的了解，並有相應的測量經驗。當沒有經驗時，可以用一個基準進行參考，該基準由Cohen(1988)提出，可為各種統計檢驗劃分小、中、大三種效應值：

統計方法	效應測量值	建議的效應值基準
		小	中	大
t檢驗	d	0.20	0.50	0.80
方差分析	f	0.10	0.25	0.40
線性模型	f2	0.02	0.15	0.35
比例檢驗	h	0.20	0.50	0.80
卡方檢驗	w	0.10	0.30	0.50

 

但是要注意，這個參考僅僅是參考，是一般性建議，特殊的研究領域內可能不適用。

 

 

 

繪製功效分析圖形：使用for循環將樣本量與相關係數之間的關係可視化出來，用圖來判斷需要的樣本量。

 

其他專業化的功效分析軟體包：

軟體包	目的
asypow	通過漸進似然比方法計算功效
PwrGSD	組序列設計的功效分析
pamm	混合模型中隨機效應的功效分析
powerSurvEpi	流行病研究的生存分析中功效和樣本量的計算
powerpkg	患病同胞配對法和TDT（Transmission Disequilibrium Test, 傳送不均衡檢驗）設計的功效分析
powerGWASinteraction	GWAS交互作用的功效計算
pedantics	一些有助於種群基因研究功效分析的函數
gap	一些病例隊列研究設計中計算功效和樣本量的函數
ssize.fdr	微陣列實驗中樣本量的計算