方差分析的基本思想及應(yīng)用條件

方差分析的基本思想

在進行科學(xué)研究時比规，有時要按實驗設(shè)計將所研究的對象分為多個處理組進行不同的處理，其中處理因素(treatment)至少有兩個水平（level）兵迅。這類科研資料的統(tǒng)計分析含友，是通過所獲得的樣本信息來推斷各處理組均數(shù)間的差別是否有統(tǒng)計學(xué)意義，即處理是否有影響髓棋。常用采用的分析方法就是方差分析（ANOVA坑雅，analysis of variance）辈挂，這是由英國統(tǒng)計學(xué)家R.A.Fisher首創(chuàng)，以F命名裹粤，故方差分析又稱為F檢驗终蒂。

設(shè)處理因素有個不同水平，實驗對象隨機分為組遥诉，分別接受不同水平的干預(yù)拇泣，第組的樣本含量為，第處理組的第個觀測值用來表示突那，其計算結(jié)果可能可以整理成以下面的形式挫酿，如下所示：

image

方差分析的目的就是在 成立的條件下构眯，通過分析各處理組均數(shù)之間的差別大小愕难，推薦個總體均數(shù)間有無差別，從面說明處理因素的效應(yīng)是否存在惫霸。

記總均數(shù)為猫缭，各處理組均數(shù)為,總例數(shù)為，g為處理組數(shù)壹店。

實驗數(shù)據(jù)有三個不同的變異：
1. 總變異猜丹。全部觀測值大小不同，這種變異稱為總變異硅卢∩渲希總變異的大小可能用離均差平方和(sum of squares of deviations from eman,SS)來表示，即各觀測值與總均數(shù)差值的平方和将塑，記為脉顿。公式略。

2. 組間變異点寥。各處理組由于接受處理的水平不同艾疟，各組的樣本均數(shù)也大小不等，這種變異稱為組間變異，其大小用各組均數(shù)與總均數(shù)的離均差平方和表示蔽莱，記為弟疆，計算公式略。

各組均數(shù)之間相關(guān)越懸殊盗冷，它們與總均數(shù)的差值越在在怠苔，就越大，反之就越小正塌。反應(yīng)了各組均數(shù)的變異嘀略，存在這種變異的原因有：①隨機誤差；②處理的不同水平可能對實驗結(jié)果的影響乓诽。

3. 組內(nèi)變異帜羊。在同一處理組中，雖然每個實驗對象接受的處理相同鸠天，但觀測值仍各不相同讼育，這種變異稱為組內(nèi)變異（誤差）。組內(nèi)變異用組內(nèi)各觀測值與其所在組的均數(shù)的差值的平方和表示稠集，記為奶段，表示隨機誤差的影響。公式略剥纷。

總離均差平方和分解為組間離均差平方和與組內(nèi)離均差平方和痹籍，就有了以下公式：

各離均差平方和的自由度為：

變異程序除與離均差平方和的大小騰外晦鞋，還與其自由度有關(guān)蹲缠，由于各部分自由度不相等，因此積分離均差平方和不能直接比較悠垛，須將各部分離均差平方和除以相應(yīng)的自由度线定，其比值稱為均方差，簡稱為均方（mean square, MS)确买。公式為：

如果各組樣本的總體均數(shù)相等()斤讥，即各處理組的樣本來自相同的總體，無處理因素的作用湾趾，則組間變異同組內(nèi)變異一樣芭商，只反映隨機誤差作用的大小，組間均方與組內(nèi)無方的比值稱為F統(tǒng)計量搀缠，如下所示：

如果F值接近于1铛楣，就沒有理由拒絕;反之，F(xiàn)值越大胡嘿，拒絕的理由越充分蛉艾，數(shù)理統(tǒng)計的理論證明，當(dāng)成立時，F(xiàn)統(tǒng)計量服從F分布勿侯，方差分析是單側(cè)F檢驗拓瞪。

變異是方差分析的基本思想

上面的話可能不太好理解。現(xiàn)在用大白話來理解一下助琐，例如我們要研究某個化合物是否有改善肥胖的效果祭埂，我們最初肯定是要做動物實驗，動物實驗的話兵钮，例如采用C57的小鼠蛆橡，分為5組，第1組掘譬，給生理鹽水泰演，第2組，給減肥藥（相當(dāng)于陽性對照）葱轩，第3組睦焕，給高劑量的化合物，第4組靴拱，給中劑量的化合物垃喊，第5組，給低劑量的化合物袜炕。分別喂一段時間后本谜，我們發(fā)現(xiàn)小鼠的體重有所變化，這個變化由兩部分構(gòu)成偎窘，第一個就是外界的刺激因素導(dǎo)致的乌助，第二種就是小鼠自身導(dǎo)致的。這種變化我們可以稱為變異（variance）评架，方差分析研究的本質(zhì)就是這種變異（體重的變化眷茁，不是基因變異那種變異）炕泳，方差分析的英語就是Analysis of variance纵诞，如果外界的刺激的因素導(dǎo)致的變異程度遠遠大于小鼠自身因素導(dǎo)致的變異，那么我們就可以認為培遵，導(dǎo)致小鼠體重這種變異是由外界刺激引起的浙芙。

不過這樣還有一個問題，因為數(shù)據(jù)越多籽腕，變異程度就越大嗡呼，為了解決這個問題，就需要用變異除以自由度（例數(shù)-1）皇耗，這樣比較的就是平均的變異南窗，因此方差分析中就出現(xiàn)了均方（MS）和組內(nèi)均方的概念。組間均方/組內(nèi)均方就是通常所說的F值，實際上代表了這樣一個含義：如果組間變異遠遠大于組內(nèi)變異万伤，那么組間均方除以組內(nèi)均方的值肯定很大窒悔，反之，這一值就會很小敌买。但是简珠，到底大到什么程度才認為有統(tǒng)計學(xué)意義呢，那就得根據(jù)F分布來判斷虹钮。

方差分析的應(yīng)用條件

多個樣本均數(shù)比較的方差分析其應(yīng)用條件為：

①各樣本是相互獨立的隨機樣本聋庵，均來自正態(tài)分布總體；

②相互比較的各樣本的總體方差相等芙粱，即具有等方差齊性祭玉；

③每個樣本之間都是獨立的。

R中的方差分析函數(shù)

所用的函數(shù)為aov()：
語法為：aov（formula,data=dataframe）
其中formula可使用的特殊符號如下春畔，其中y為因變量攘宙，A、B拐迁、C為自變量：

單因素方差分析

單因素方差分析（one-way ANOVA）是指對單因素試驗結(jié)果進行分析，檢驗因素對試驗結(jié)果有無顯著性影響的方法讯壶。單因素方差分析是用來檢驗多個平均數(shù)之間的差異料仗，從而確定因素對試驗結(jié)果有無顯著性影響的一種統(tǒng)計方法。對于完全隨機設(shè)計試驗且處理數(shù)大于2時可以用單因素方差分析（等于2 時用t檢驗）伏蚊。離差平方和的分解公式為：SST(總和)=SSR(組間)+SSE(組內(nèi))立轧，F(xiàn)統(tǒng)計量為MSR/MSE，MSR=SSR/k-1,MSE=SSE/n-k躏吊。其中SST為總離差氛改、SSR為組間平方和、SSE為組內(nèi)平方和或殘差平方和比伏、MSR為組間均方差胜卤、MSE為組內(nèi)均方差。

ANOVA-案例A

單因素方差分析例4-2：某醫(yī)生為了研究一種降血脂新藥的臨床療效赁项，按統(tǒng)一納入標準選擇120名高血脂患者葛躏，采用完全隨機設(shè)計方法將患者等分為4組（具體分組方法見例4-1）澈段，進行雙盲試驗。6周后測得低密度脂蛋白作為試驗結(jié)果舰攒，見表4-3均蜜。問4個處理組患者的低密度脂蛋白含量總體均數(shù)有無差別?
數(shù)據(jù)如下所示：

image

計算過程如下所示：

1. 導(dǎo)入數(shù)據(jù)

anova1 <- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_szq_402.csv",sep=",")
library(reshape2)
anova1 <- melt(anova1)

2.正態(tài)檢驗

方差分析需要一定的假設(shè)，即數(shù)據(jù)集應(yīng)該符合正態(tài)和各組的方差相等芒率，可以分別用shapiro.test和bartlett.test檢驗從P值觀察到這兩個假設(shè)是符合的囤耳。對于不符合假設(shè)的情況，我們就要用到非參數(shù)方法偶芍，例如Kruskal-Wallis秩和檢驗

group <- c("control","high","middle","control")

for(i in group){
  print(shapiro.test(anova1[anova1$variable==i,]$value))
}

結(jié)果如下所示：

    Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.95229, p-value = 0.1947


    Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.95833, p-value = 0.2806


    Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.9445, p-value = 0.1202


    Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.95229, p-value = 0.1947

P值大于0.05說明數(shù)據(jù)正態(tài)P值大于0.05說明數(shù)據(jù)正態(tài)

3. 方差齊性檢驗:

方差齊性檢驗就是檢驗各組樣本所代表的總體方差是否一致的檢驗充择，兩樣本方差齊性檢驗使用Bartlett法，同樣匪蟀，它也適用于多樣本的方差齊性檢驗椎麦，它是它要求所檢驗的樣本總體符合正態(tài)分頁，當(dāng)不符合正態(tài)分布的時候材彪，就不能使用观挎，則要用Levene檢驗。Levene檢驗不受數(shù)據(jù)頒販限制段化，是一種穩(wěn)健的檢驗嘁捷，因而被廣泛地認為是一種標準的檢驗方差齊性的檢驗。

方差齊性通常用bartlett檢驗

bartlett.test(anova1$value~anova1$variable)

結(jié)果如下所示：

image

結(jié)果顯示p值大于0.05显熏，可認為方差齊性雄嚣。

或者用levene檢驗：

Levene檢驗既可以用于正態(tài)分布的數(shù)據(jù)，也可以用于非正態(tài)分布的數(shù)據(jù)或分布不明的數(shù)據(jù)喘蟆，具有比較穩(wěn)健的特點缓升，檢驗效果也比較理理想。

library(car)
leveneTest(anova1$value~anova1$variable)

結(jié)果如下所示：

image

結(jié)果發(fā)現(xiàn)sig值大于0.05,表明符合方差齊性假設(shè),可以進行進一步的參數(shù)檢驗蕴轨。

4. 檢驗整體均值是否有差異

方差例如技術(shù)的目的就是要比較因素A的r水平上港谊，試驗結(jié)果是否有顯著差異，以樣本作為檢驗的標準橙弱，寫出檢驗假設(shè)歧寺，如下所示：

，

如果拒絕原假設(shè)膘螟，則說明樣本來自不同的正態(tài)總體成福，則由因素A的各個水平所造成的均值差異有統(tǒng)計意義碾局；若不能拒絕原假設(shè)荆残，說明樣本來自相同的正態(tài)總體，因素的 不同水平之間無差異净当。在這個案例中内斯，我們研究的就是不同劑量的藥物對低密度脂蛋白的水平是否產(chǎn)生影響蕴潦，我們可以使用3個函數(shù)來進行計算。

aov()

aov函數(shù)的使用如下所示：

aov(formula,data=NULL,projections=FALSE,qr=TRUE,contrasts=NULL,...)

其中formula表示主方差分析的公式俘闯，在單因素方差分析中潭苞，即為x~A；data表示方差分析的數(shù)據(jù)框真朗；projections為邏輯值此疹，表示是否返回預(yù)測結(jié)果；qr是邏輯值遮婶，表示是否返回QR分解結(jié)果蝗碎，默認為TRUE；contrasts是公式中的一些因子的對比列表旗扑，計算過程如下所示：

result <- aov(value~variable,data=anova1)
summary(result)

結(jié)果如下所示：

image

其中Residuals是殘差蹦骑，variable代表不同的組，p值小于0.001臀防，因此各組之間存在顯著性差異眠菇。另外，R給出的F值是24.88袱衷，而書中的例子是24.93捎废，書中的值是由查F表得來的，是個范圍致燥，R中的是具體值缕坎。

oneway.test()

也可以采用oneway.test()函數(shù)進行檢驗，

oneway.test()函數(shù)的用法如下所示：

oneway.test(formual,data,subset,na.action,var.equal=FALSE)

若各水平下的總體方差相等篡悟，且參數(shù)var.equal=TRUE谜叹，那么等于同函數(shù)aov()所作的方差分析結(jié)果。在默認情況下搬葬，各水平的總體方差不相等(var.equal=FALSE)荷腊，則使用Welch的近似方法，在任何多個樣本的情況下作兩樣本W(wǎng)elch檢驗急凰，在滿足方差齊發(fā)的條件下女仰，使用oneway.test()的計算過程如下所示：

> result2 <- oneway.test(value~variable,data=anova1,var.equal=TRUE)
> result2

    One-way analysis of means

data:  value and variable
F = 24.884, num df = 3, denom df = 116, p-value = 1.674e-12

lm()

單因素方差分析的數(shù)學(xué)模型可以看作是一種特殊的線性模型，因此方差分析還可以通過線性模型函數(shù)lm()計算得到抡锈，再利用anova()提取其中的方差分析表疾忍，因此aov(formula)等價于anova(lm(formula))，如下所示：

> anova(lm(value~variable,data=anova1))
Analysis of Variance Table

Response: value
           Df Sum Sq Mean Sq F value    Pr(>F)    
variable    3 32.156 10.7187  24.884 1.674e-12 ***
Residuals 116 49.967  0.4308                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5. 均數(shù)間的兩兩比較

方差分析得出總體之間有差異床三，要進一步知道哪兩組之間有差異一罩，就要使用均數(shù)間的多重比較，常用的比較方法有SNK檢驗（q檢驗）撇簿，LSD檢驗聂渊，Bonferroni檢驗差购，Dunnett檢驗，TurkeyHSD檢驗汉嗽，下面是對一些常用的兩兩比較方法的匯總：

5.1 常見的兩兩比較方法

LSD：相當(dāng)于t檢驗欲逃，只不過它需要在方差分析一定要有統(tǒng)計學(xué)差異的情況下才用。所以LSD法并沒有控制假陽性錯誤饼暑。一般情況下稳析，如果你在設(shè)計初期就有很明確的目的，可以考慮這種方法弓叛，因為每一對比較都是有特定意義的迈着，不用非得控制假陽性錯誤。

SNK法是先按多組均值大小排序邪码，然后按一個有點類似于t檢驗的公式分別比較（不過誤差計算不同）裕菠。比如a、b闭专、c三組均值分別是a最小奴潘，c最大，b居中影钉。那么比較時画髓，很顯然a和c差別最大，所以在最后的界值上做一些調(diào)整平委。例如奈虾，原來可能t檢驗的界值是2.83，只要t值大于2.83就算有差異廉赔，現(xiàn)在把界值調(diào)大肉微，比如調(diào)到3.40，a和c的比較計算的t值蜡塌，只有大于3.4才算有差異碉纳。這種調(diào)整的界值叫做q界值。這種方法對假陽性錯誤的控制依然不到位馏艾。所以如果你想控制假陽性劳曹，不要選這種方法。

Bonferroni法大概是最流行的一種方法琅摩，因為簡單铁孵。它的思想是調(diào)整檢驗水準，根據(jù)比較的次數(shù)重新設(shè)定檢驗水準房资，然后根據(jù)P值做出結(jié)論蜕劝。比如常規(guī)的檢驗水準是0.05，只要P小于0.05就認為有統(tǒng)計學(xué)差異志膀。但是如果用Bonferroni法調(diào)整熙宇，則需要0.05除以比較次數(shù)鳖擒，如比較6次溉浙，這時調(diào)整后的檢驗水準是0.05/6=0.0083宏胯，也就是說甥厦，P值小于0.0083才算有差異。這種方法其實更像是一種檢驗水準調(diào)整方法，簡單易用吃靠，不僅可以用于均值的兩兩比較，也可用于率的兩兩比較烙样。比如你要做三組的卡方檢驗开仰，想進一步兩兩比較，就可以考慮用這種方式颂郎，即分別做兩組的卡方吼渡，但是一共比較3次，所以P值小于0.0167才算有差異乓序。這種方法的缺點是寺酪，當(dāng)比較次數(shù)多時，結(jié)果過于保守替劈，比如比較10次寄雀，那就需要P值小于0.005才算有差異，有時很難達到陨献。

再說點其它國內(nèi)課本上不大介紹但是國外介紹比較多的幾種方法盒犹。

Tukey法，有時也叫Tukey HSD法（Honestly Significant Difference test）眨业，這種方法大概是統(tǒng)計學(xué)家最認可一致的方法了急膀，絕大多數(shù)統(tǒng)計學(xué)家都相信，Tukey法的檢驗效率可能是最高的龄捡。Tukey法也是基于q檢驗脖阵，大概意思是先確定一個最大差異的臨界值，然后分別對其中兩組比較墅茉，看看哪兩組差值大于這個界值命黔，就算有差異。Tukey法是大多數(shù)統(tǒng)計學(xué)家首先推薦的兩兩比較方法就斤，不過這種方法只適用于組間例數(shù)相等的情況悍募。對于組間例數(shù)不等的時候，可用修正的Tukey法洋机，也叫作Tukey-Kramer法坠宴。

現(xiàn)在計算各組之間的均數(shù)與SD

aggregate(anova1$value,by=list(anova1$variable),FUN=mean)
aggregate(anova1$value,by=list(anova1$variable),FUN=sd)

結(jié)果如下所示：

image

繪制箱形圖可能觀察到不同因素對于因變量的影響

plot.design(value ~ variable,data =anova1, main = "Group means") 
# plot.design was inclued in packages "graphics"

image

繪制有置信區(qū)間的組均值圖

library(gplots)
plotmeans(value ~ variable,data =anova1 )

圖片如下所示：

image

單因素方差分析是從總體的角度上說明各效應(yīng)的均值之間存在著顯著差異，但具體在哪些水平下的均值存在較大差異绷旗，還不清楚喜鼓，因此我們要對每一對樣本均值進行一一比較副砍，既要進行組間均值的兩兩比較。最常用庄岖，也是最直接的方法就是多重t檢驗方法豁翎，其對因子A每個水平下的數(shù)據(jù)兩兩比較均值進行t檢驗，但估計方差時用的是全部數(shù)據(jù)的誤差均方和隅忿，檢驗的假設(shè)為：

構(gòu)建t檢驗統(tǒng)計量心剥，如下所示：

image

在原假設(shè)成立的情況下，該統(tǒng)計量服從自由度為n-r的t分布背桐，所以檢驗拒絕為优烧。

多重t檢驗的思路簡單，但是這樣的方法不能直接使用链峭，這是因為在樣本之間多次重復(fù)使用t檢驗會大大增加犯第一類錯誤的概率畦娄。例如因子A有三個水平，需要比較3次(n=k(k-1)/2)弊仪，顯著性水平為0.05熙卡，回此每次比較犯第一類錯誤的概率就是0.05，那么截取次比較同時進行撼短，犯第一次錯誤的總概率為再膳，這個數(shù)字就比較大了，因此直接使用多重t檢驗的結(jié)果不一定是可靠的曲横，統(tǒng)計學(xué)家們提供了很多方法對P值進行修正喂柒，R中的P值修正函數(shù)是p.ajdust()，其調(diào)用格式如下所示：

p.adjust(p,method=p.ajdust.methods,n=length(p))

在R中禾嫉，輸入如下指令可以看到調(diào)整方法的列表灾杰，如下所示：

> p.adjust.methods
[1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"         "BY"         "fdr"       
[8] "none" 

當(dāng)多重檢驗次數(shù)較多時，bonferroni修正方法的效果比較好熙参，思路也很簡單：如果在同一個數(shù)據(jù)集上同時進行n個獨立的假設(shè)檢驗艳吠，那么用于每一個假設(shè)的統(tǒng)計顯著水平，應(yīng)為僅檢驗一個假設(shè)時顯著水平的1/n孽椰，即昭娩，其中n為多重t檢驗的次數(shù)。

在R中計算均值的多重t檢驗時黍匾，用函數(shù)pairwise.t.test()返回多重比較的修正P值栏渺，其格式如下所示：

pairwise.t.test(x, g, p.adjust.method = p.adjust.methods,
                pool.sd = !paired, paired = FALSE,
                alternative = c("two.sided", "less", "greater"),
                ...)

其中參數(shù)x是響應(yīng)向量，g是因子向量锐涯，p.adjust.methods給出了P值的修正方法磕诊，默認按Holm方法修正，若不作任何調(diào)整則設(shè)為none；paired為邏輯值霎终，指示是否要配對t檢驗滞磺；alternative指明檢驗的方向問題。現(xiàn)在對前面的結(jié)果進行多重t檢驗莱褒，使用bonferroni修正方法击困，如下所示：

> pairwise.t.test(anova1$value,anova1$variable,p.adjust.method = "bonferroni")

    Pairwise comparisons using t tests with pooled SD 

data:  anova1$value and anova1$variable 

       control low     middle 
low    0.00029 -       -      
middle 0.00020 1.00000 -      
high   2.1e-13 0.00013 0.00020

P value adjustment method: bonferroni 

經(jīng)過修正后的p值會比原理增大很多，這在一定程度上克服了多重t檢驗增加犯第一類錯誤概率的缺點保礼，從上面結(jié)果來看沛励，除了middle與low組外责语，其余組之間的p值都很小炮障，說明這幾個樣本之間有差異。

再看其它的兩兩比較方法坤候，這些方法有LSD胁赢，TukeyHSD，Scheffe檢驗白筹，

LSD檢驗

library(agricolae) #此包中有LSD檢驗
result <- LSD.test(result,"variable",p.adj="bonferroni")
result

結(jié)果如下所示：

image

注：R給出的F值是24.88智末，而書中的例子是24.93，書中的值是由查F表得來的徒河，是個范圍系馆，R中的是具體值。
Bonferroni校正("bonferroni")顽照，用于多重比較的p值校正由蘑，次數(shù)不多時適用。如果在同一數(shù)據(jù)集上同時檢驗n個獨立的假設(shè)代兵，那么用于每一假設(shè)的統(tǒng)計顯著水平尼酿，應(yīng)為僅檢驗一個假設(shè)時的顯著水平的1/n。舉個例子：如要在同一數(shù)據(jù)集上檢驗兩個獨立的假設(shè)植影，顯著水平設(shè)為常見的0.05裳擎。此時用于檢驗該兩個假設(shè)應(yīng)使用更嚴格的0.025。即0.05* (1/2)思币。該方法是由Carlo Emilio Bonferroni發(fā)展的鹿响，因此稱Bonferroni校正。在藥理實驗中谷饿，動物數(shù)通常不會太多惶我，用Bonferroni校正的居多。

TukeyHSD檢驗

result2 <- aov(value~variable,data=anova1)
result.tukey <- TukeyHSD(result2)
result.tukey
plot(result.tukey)

結(jié)果如下所示：

image

圖片結(jié)果：

image

glht()函數(shù)做Tukey檢驗

此外各墨，multcomp包中glht()函數(shù)提供了多重均值更全面的方法指孤，適用于線性模型和廣義線性模型，下面的代碼重現(xiàn)Tukey HSD檢驗，如下所示：

library(multcomp)
result4 <- aov(value ~ variable, data = anova1)
tukey4 <- glht(result4, linfct=mcp(variable="Tukey"))
summary(tukey4)

計算結(jié)果如下所示：

> summary(tukey4)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = value ~ variable, data = anova1)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)    
low - control == 0    -0.71500    0.16946  -4.219 0.000291 ***
middle - control == 0 -0.73233    0.16946  -4.322 0.000171 ***
high - control == 0   -1.46400    0.16946  -8.639  < 1e-04 ***
middle - low == 0     -0.01733    0.16946  -0.102 0.999615    
high - low == 0       -0.74900    0.16946  -4.420 0.000108 ***
high - middle == 0    -0.73167    0.16946  -4.318 0.000170 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

繪制不同組的箱線圖恃轩，如下所示：

plot(cld(tukey4,level = .05),col="lightgrey")
#cld()函數(shù)中l(wèi)evel選項設(shè)置了使用顯著水平0.05结洼，即95%的置信區(qū)間

image

7.殘差分析

這一步做殘差分析就是為了再次驗證原始數(shù)據(jù)是否服從正態(tài)分布，如下所示：

residual <- rstudent(result4)
qqnorm(residual, pch=20, cex=2)
qqline(residual, col="gray60", lwd=2)

殘差的QQ圖如下所示：

image

shapiro檢驗叉跛，如下所示：

> shapiro.test(residual)

    Shapiro-Wilk normality test

data:  residual
W = 0.9884, p-value = 0.4026

p值的為0.4026松忍，可以認為殘差滿足正態(tài)性。

繪制殘差圖筷厘，如下所示：

plot(residual ~ anova1$variable, main="各組的殘差圖")

image

對殘差進行方差齊性檢驗鸣峭，如下所示：

> library(car)
> leveneTest(result4)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   3   1.493 0.2201
      116

P值大于0.05，可以認為殘差滿足方差齊性酥艳。

ANOVA-案例B

再看一個案例摊溶，進行一次完整的ANOVA分析，如下所示：

某單位研究不同藥物對小白鼠的鎮(zhèn)咳作用充石，抽取40只小白鼠隨機分配到各個藥物組中莫换，實驗時先用0.2mL NH4OH對小白鼠噴霧，測定其發(fā)生咳嗽時間骤铃。以給藥前后發(fā)生咳嗽時間的差值衡量不同藥物的鎮(zhèn)咳作用拉岁，結(jié)果見表5.1.試比較3種藥物的平均推遲咳嗽時間的差異有無差異。（王炳順. 醫(yī)學(xué)統(tǒng)計學(xué)及SAS應(yīng)用[M]. 上海交通大學(xué)出版社, 2007.）

計算過程如下所示：

anova2 <- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_anova_01.csv",sep=",")
# input raw data 
# 王炳順. 醫(yī)學(xué)統(tǒng)計學(xué)及SAS應(yīng)用[M]. 上海交通大學(xué)出版社, 2007.例5.1

anova2$group <- as.factor(anova2$group)
library(car)
leveneTest(anova2$time~anova2$group)
result <- aov(anova2$time~anova2$group)
summary(result)
compareresult <- TukeyHSD(result)
compareresult
plot(compareresult)

結(jié)果如下所示：

> for(i in 1:3){
+   result <- shapiro.test(filter(anova2,group==i)$time)
+   print(result)
+ }

    Shapiro-Wilk normality test

data:  filter(anova2, group == i)$time
W = 0.9778, p-value = 0.9523


    Shapiro-Wilk normality test

data:  filter(anova2, group == i)$time
W = 0.98349, p-value = 0.9878


    Shapiro-Wilk normality test

data:  filter(anova2, group == i)$time
W = 0.94005, p-value = 0.5536

> 
> # test for homogeneity var
> leveneTest(anova2$time~anova2$group)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  1.4002 0.2593
      37               
> 
> ## test for mean of all samples
> result <- aov(anova2$time~anova2$group)
> summary(result)
             Df Sum Sq Mean Sq F value Pr(>F)  
anova2$group  2   5062  2531.2   3.485 0.0411 *
Residuals    37  26877   726.4                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> ## GLM for mean of all samples
> anova(lm(anova2$time~anova2$group))
Analysis of Variance Table

Response: anova2$time
             Df  Sum Sq Mean Sq F value  Pr(>F)  
anova2$group  2  5062.5 2531.23  3.4845 0.04107 *
Residuals    37 26877.4  726.42                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> ## test for each pair of samples
> ## TukeyHSD
> TukeyHSD(result)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = anova2$time ~ anova2$group)

$`anova2$group`
        diff        lwr      upr     p adj
2-1 12.33333 -11.694604 36.36127 0.4302024
3-1 29.03333   2.169283 55.89738 0.0317315
3-2 16.70000 -10.164050 43.56405 0.2944100

> pairwise.t.test(anova2$time,anova2$group,p.adjust.method = "bonferroni")

    Pairwise comparisons using t tests with pooled SD 

data:  anova2$time and anova2$group 

  1     2    
2 0.654 -    
3 0.036 0.413

P value adjustment method: bonferroni 

重復(fù)測量方差分析

未設(shè)立平等對照的前后測量

image

... ...

[圖片上傳失敗...(image-cd09b0-1561781157663)]

雖然這個表和配對設(shè)計的t檢驗結(jié)果表達完全相同，但卻是屬于兩種不同類型的設(shè)計撕瞧，其區(qū)別在于

1. 配對設(shè)計中同一對子的兩個實驗單位可以隨機分配處理陵叽，兩個實驗單位同期觀察試驗結(jié)果，哦可以比較處理組間差別风范。前后測量設(shè)計不能同期觀察試驗結(jié)果咨跌，雖然可以在前后測量之間安排處理，但本質(zhì)上比較的是前后差別硼婿，推論處理是否有效是有條件的锌半，即假定測量時間對觀察結(jié)果沒有影響。
2. 配對t檢驗要求同一對子的兩個實驗單位的觀察結(jié)果分別與差值相互獨立寇漫，差值服從正態(tài)分布刊殉。而膠后測量設(shè)計前后兩次觀察結(jié)果通常與差值不獨立，大多數(shù)情況下第一次觀察結(jié)果與差值存在負相關(guān)關(guān)系州胳，例如上表记焊，治療前舒張壓與差值的相關(guān)系數(shù)為-0.602。
3. 配對設(shè)計用平均差值推論處理的作用栓撞，前后測量測序除了分析平均差值外遍膜，還可以進行相關(guān)回歸分析碗硬，由上表計算，治療前后舒張壓的相關(guān)系數(shù)為0.963瓢颅，P<0.01恩尾，用治療前舒張壓(X)推論治療后舒張壓(Y)的回歸方向為，截矩檢驗P=0.014挽懦，回歸系數(shù)檢驗P<0.01翰意。

設(shè)立平等對照的前后測量

第1張表中的高血壓患者治療后的舒張壓平均下降了16mmHg，雖然經(jīng)配對t檢驗信柿，t=16.18冀偶，P<0.01，這也未必說明治療有效渔嚷，因為住院休息进鸠，環(huán)境和情緒的改變同樣可以使血壓恢復(fù)平衡。因此圃伶，確定療效的前后測量設(shè)計必須增加平行對照堤如，例如將20位輕度高血壓患者隨機分配到處理組和對照組蒲列，試驗結(jié)果見表12-2.如下所示：

image

比較這張表中處理組與對照組之間的差別窒朋，似乎可以用前后測量的差值(d)做兩組均數(shù)比較的t檢驗，但上表中處理組與對照組的差值d方差不齊（F值為6.58蝗岖，P<0.01）侥猩，不符合兩無數(shù)比較t檢驗的前提條件，更為重要的是抵赢，如果處理組與對照組基線值（治療前的舒張壓）不同欺劳，那么差值不能作為判別組間差別的依據(jù)。

重復(fù)測試設(shè)計

當(dāng)前后測量設(shè)計的重復(fù)測量次數(shù)m>=3時铅鲤，稱為重復(fù)測試設(shè)計或重復(fù)測量數(shù)據(jù)划提，重復(fù)測試數(shù)據(jù)類似于隨機區(qū)間設(shè)計數(shù)據(jù)，如下表所示邢享，而且同樣可以計算出隨機區(qū)組設(shè)計的方差分析表：

案例分析B

image

下表是隨機區(qū)組數(shù)據(jù)的方差分析表：

image

但重復(fù)測量設(shè)計與隨機區(qū)級設(shè)計有區(qū)別鹏往，它們的區(qū)別在于：

（一）重復(fù)測量設(shè)計中的“處理”是在區(qū)組（受試者）間隨機分配，區(qū)組內(nèi)的各時間點是固定的骇塘，不能隨機分配伊履，如下表所示：

image

A、B兩種處理隨機分配給各個患者后款违，每個患者測量的時間是相同的唐瀑，隨機區(qū)組設(shè)計要求每個區(qū)組內(nèi)實驗單位彼此獨立，處理只能在區(qū)組內(nèi)隨機分配插爹，每個實驗單位接受的處理是不相同的哄辣，如下所示（隨機區(qū)間設(shè)計）：

image

（二）重復(fù)測量設(shè)計區(qū)組內(nèi)實驗單位彼此不獨立请梢，如在表12-3中，每個受試者血糖濃度的個體特征是用4個時間點的測量值刻畫的力穗，如下圖所示：

[圖片上傳失敗...(image-c425d8-1561781157663)]

同一受試者的血樣重復(fù)測量結(jié)果是高度相關(guān)的溢陪，相關(guān)系數(shù)如下所示：

image

更準確地說，表12-3重復(fù)測量數(shù)據(jù)用隨機區(qū)組方差分析比較處理組間差異睛廊，前提條件是要滿足“球?qū)ΨQ(sphericity)”假設(shè)形真，即重復(fù)測量誤差的協(xié)方差矩陣經(jīng)正交對比變換后，與單位矩陣成比較超全，表12-3數(shù)據(jù)Mauchly“球?qū)ΨQ”檢驗結(jié)果見表12-7咆霜，其中，拒絕“球?qū)ΨQ”假設(shè)嘶朱，即表12-3的數(shù)據(jù)不滿足隨機區(qū)組方差分析（表12-4）“球?qū)ΨQ”的條件蛾坯，如下所示：

image

此時，就需要對組內(nèi)效應(yīng)的F界值進行校正疏遏，校正的方法是用“球?qū)ΨQ”系數(shù)乘以組內(nèi)效應(yīng)F界值的自由度和脉课，得到，用作為檢驗界值财异，“球?qū)ΨQ”系數(shù)的常用估計方法有Greenhouse-Geisser倘零、Huynh-Feldt和Lower-bound三種方法，例如戳寸，表12-4組內(nèi)效應(yīng)的F界值為呈驶，Greenhouse-Geisser的校正系數(shù)0.536校正后的F界值為，大于未校正的界值3.07疫鹊，也就是說當(dāng)重復(fù)測量數(shù)據(jù)不滿足”球?qū)ΨQ“假設(shè)時袖瞻，采用隨機區(qū)組設(shè)計方差分析時，增大了I型錯誤（無差別判斷為有差別）的概率拆吆，當(dāng)樣本含量較小時聋迎，Mauchly”球?qū)ΨQ“檢驗交通較低，易發(fā)生II型錯誤（不滿足”球?qū)ΨQ“判斷為滿足）枣耀，此時霉晕，對處理組內(nèi)效應(yīng)F界值的自由度和仍要進行校正，以得到可信的統(tǒng)計結(jié)論奕枢。

重復(fù)測量數(shù)據(jù)的兩因素兩水平分析

兩因素離均差平方和分解

將重復(fù)測量測序的干預(yù)因素作為A因素娄昆，共兩個水平，1水平為”對照“缝彬，2水平為”干預(yù)“萌焰；前后兩次測量時間作為B因素，共兩個水平谷浅，1水平為”第1次測量時間“扒俯，2水平為”第2次測量時間“奶卓，這樣共有(a1b1)、(a1b2)撼玄、(a2b1)夺姑、(a2b2)4個”處理“組，各組觀察值的合計分別為T1掌猛、T2盏浙、T3、T4表示荔茬，A因素兩水平的小計分別用A1废膘、A2表示，B因素兩水平的小計分別用B1慕蔚、B2表示丐黄，處理組和對照組的例數(shù)相等，即n=n1=n2孔飒，SS處理可分解為A因素的離均差平方和灌闺，B因素的離均差平方和，AB交互作用的離均差平方和三部分坏瞄，見表12-8桂对，如下所示：

image

將表12-2的數(shù)據(jù)按照表12-8進行兩因素離均差平方和進行分解，如下所示：

image

[圖片上傳失敗...(image-44e062-1561781157663)]

主效應(yīng)與交互作用的離均差平方和分解：

image

組間惦积、組內(nèi)離均差平方和分解與方差分析

重復(fù)測量數(shù)據(jù)的差異由兩部分組成接校。一是觀察對象的個體差異，其離均差平方和記作SS組間狮崩，二是每個觀察對象前后兩次觀察之間的差異，其離均差平方和記作SS組內(nèi)鹿寻，設(shè)兩組觀察對象都為n睦柴，每個觀察對象前后兩次觀察的合計值為Mj，SS組間和SS組內(nèi)的計算方法見表12-9毡熏，如下所示：

image

SS組間坦敌、SS組內(nèi)分解完畢后，就可以用方差分析對A因素主效應(yīng)痢法、B因素主效應(yīng)和AB因素交互作用進行假設(shè)檢驗了狱窘，由于A因素（干預(yù)分組）的作用于觀察以剛，作用效應(yīng)在SS組間财搁，方差分析按表12-10進行蘸炸，B因素主效應(yīng)（測量前后）和AB交互作用的效應(yīng)在SS組內(nèi),方差分析按表12-11進行，如下所示：

image

現(xiàn)在根據(jù)表12-2中的數(shù)據(jù)尖奔，對處理組與對照組搭儒、治療前后舒張壓的差別進行統(tǒng)計分析穷当。

1.計算、

表12-2共有20名患者淹禾，每個患者舒張壓的合計為M1=130+114=244馁菜，M2=124+110=234，...铃岔，M19=120+124=244汪疮，M20=134+128=262，由于前面已經(jīng)計算出C=580328.1毁习，按表12-9中的公式計算铲咨，如下所示：

image

2.測量前后比較與交互作用分析

由于已經(jīng)計算出了，現(xiàn)在SS組內(nèi)=1702蜓洪，代入表12-11纤勒，得到方差分析表，如下所示：

image

3.處理組與對照組比較

前面已經(jīng)計算出隆檀，現(xiàn)在SS組間=2517.9摇天，代入表12-10，得到方差分析表恐仑，如下所示：

image

4.結(jié)論

①由表12-12的方差分析表可知泉坐，治療前后舒張壓的交效應(yīng)有差別(P<0.01)，由于表12-2計算可知裳仆，治療前的平均舒張壓為(126.2+124.8)/2=125.5mmHg腕让，治療后的平均舒張壓為(110.2+120.6)/2=115.4mmHg，治療前后舒張壓主效應(yīng)的估計值為115.4-125.5=-9.6mmHg歧斟。①由表12-13方差分析表纯丸，不考慮測量時間，處理組與對照組舒張壓的主效應(yīng)未見差別(P>0.05)静袖，由原始數(shù)據(jù)過計算觉鼻，處理組的平均舒張壓為(126.2+110.2)/2=118.2mmHg，對照組的平均舒張壓為(124.8+120.6)/2=122.7mmHg队橙，處理組與對照組舒張壓主效應(yīng)的估計值為118.2-122.7=-4.5mmHg坠陈。③由表12-12方差分析表可知，測量前后與處理存在交互作用(P<0.01)捐康，即處理組和對照組治療前后的舒張壓的變化幅度不同仇矾。由單獨效應(yīng)、主效應(yīng)和交互作用的概念可知解总，若存在交互作用贮匕，單獨分析主效應(yīng)的意義不大，需要愛一分析各因素的單位效應(yīng)，按表12-2計算處理組和對照組的單獨效應(yīng)振愿，處理組和對照組治療前的基線數(shù)據(jù)沒有統(tǒng)計學(xué)差異，如下所示：

image

但處理組和舒張壓平均下降16.0mmmHg腰池，對照組平均僅下降4.2mmHg隙疚，AB交互作用F=18.8壤追，P<0.01，說明處理組的降壓效應(yīng)優(yōu)于對照組供屉。

隨機區(qū)組設(shè)計

隨機區(qū)組設(shè)計(randomized block design)又稱為配伍組設(shè)計行冰，是配對設(shè)計的擴展。具體做法是：先按影響實驗結(jié)果的非處理因素（如性別伶丐、體重悼做、年齡、職業(yè)哗魂、病情肛走、病程等）將實驗對象配成區(qū)組(block)，再分別將各區(qū)組內(nèi)的實驗對象隨機分配到各處理組或?qū)φ战M录别。與完全隨機設(shè)計相比朽色，隨機區(qū)組設(shè)計的特點是隨機同分配的次數(shù)要重復(fù)多次，每次隨機分配都對同一個區(qū)組內(nèi)的實驗對象進行组题，且各個處理組的實驗對象數(shù)量相同葫男，區(qū)組內(nèi)均衡。在進行統(tǒng)計分析時崔列，將區(qū)組變異離均差平方和從完全隨機設(shè)計的組內(nèi)離均差平方和中分離出來梢褐，從而減小組內(nèi)離均差平方和（誤差平方和），提高統(tǒng)計檢驗效率赵讯。若將區(qū)組作為另一處理因素的不同水平盈咳，隨機區(qū)組設(shè)計等同于無重復(fù)觀察的兩因素設(shè)計。

現(xiàn)在看一下隨機區(qū)組設(shè)計的案例瘦癌。

如何按隨機區(qū)組設(shè)計猪贪，分配5個區(qū)組的15只小白鼠接受甲、乙讯私、丙三種抗癌藥物？

方法為先將小白鼠的體重從輕到重編號西傀，體重相近的3只小白鼠配成一個區(qū)組（下表中第一行和第二行）斤寇，然后對這些小鼠分配隨機數(shù)，在每個區(qū)組內(nèi)將隨機數(shù)按大小排序（下表第四行）拥褂，各區(qū)組中序號為1的接受甲藥娘锁、序號為2的接受乙藥、序號為3的接受丙藥饺鹃，如下所示：

image

隨機區(qū)組設(shè)計資料在進行統(tǒng)計分析時莫秆，需要根據(jù)數(shù)據(jù)特征選擇方法间雀，對于正態(tài)分布且方差齊同的資料，采用雙向分類的方差分析(two-way classification ANOVA)或配對t檢驗(g=2)镊屎，當(dāng)不滿足方差分析或t檢驗條件時惹挟，可以進行變量變換后采用雙向分析的方差分析或采用Fredman M檢驗。

重復(fù)測量數(shù)據(jù)的兩因素多水平分析

大多數(shù)醫(yī)學(xué)實驗都有重復(fù)測量記錄缝驳，如果統(tǒng)計分析時只分析最后一次的測量結(jié)果连锯，那就會喪失很多“過程“信息，如果通過記錄不同時間點的重復(fù)觀察數(shù)據(jù)了解測量指標的時間變化趨勢用狱、藥物在血中濃度變化的時間分析特征运怖、不同時間的治療效果等。在臨床上夏伊，各個單個時間的測量值或許都在正常參考值范圍以內(nèi)（或以外）摇展，但如果發(fā)現(xiàn)該測量值持續(xù)上升（或下降），也說明患者的狀態(tài)處理變化中溺忧，在統(tǒng)計分析時咏连，保留”處理“前的基線(baseline)，如前面所述的患者治療前收縮壓砸狞，可以評價隨機分級的均衡性捻勉，也能夠提高統(tǒng)計分析的效率。

實驗設(shè)計方法

一般來說刀森，任何實驗設(shè)計都可以采用重復(fù)測量計算踱启，即在實驗過得中定期記錄觀察結(jié)果。所謂重復(fù)測量數(shù)據(jù)的兩因素多水平設(shè)計研底，兩因素指干預(yù)（A因素）和測量時間（B因素）埠偿，多水平指干預(yù)（A因素）有g(shù)（大于等于2）個水平，測量時間（B因素）有m（大于等于2）個水平（時間點）榜晦，即每個觀察對象有m個重復(fù)測量數(shù)據(jù)冠蒋，隨機化分組是將作用于觀察對象的g個干預(yù)隨機分配，前面提到的重復(fù)測量數(shù)據(jù)的兩因素兩水平是兩因素多水平的設(shè)計的特征乾胶，即g=m=2.

試驗結(jié)果的方差分析

重復(fù)測量數(shù)據(jù)的統(tǒng)計分析有許多統(tǒng)計方法供選用抖剿，可以用單因素方差分析，也可以用多變量方差分析(MANOVA)识窿，其中單因素方差分析最簡單斩郎。設(shè)有觀察對象隨機等分成g個干預(yù)組，每組倒數(shù)為n喻频，重復(fù)測量次數(shù)為m缩宜，每個觀察對象測量值合計為，g個干預(yù)組，每組的測量值合計為锻煌。多個干預(yù)組比較的方差分析用表12-14妓布，設(shè)表示第j個時間的測量值合計，表示第i個干預(yù)宋梧，第j個時間的點測量值合計匣沼，多個時間點測量前后與交互作用的方差分析見表12-15，如下所示：

image

案例分析B

將手術(shù)要求基本相同的15名患者隨機分為3組乃秀，在手術(shù)過程中分別采用A肛著、B、C三種麻醉誘導(dǎo)方法跺讯，在T0（誘導(dǎo)前）枢贿、T1、T2刀脏、T3局荚、T4五個時間點測量患者的收縮壓，數(shù)據(jù)記錄見表12-16愈污，試進行方差分析耀态，如下所示：

image

1.分解

g=3,m=5,n=5,15名患者的合計分別為：

image

2.分解

分組計算不同麻醉誘導(dǎo)首装、不同時間點患者的收縮壓的合計值()，見表12-17杭跪，如下所示：

image

3.方差分析表

按表12-14仙逻、12-15列出方差分析表。

4.F值的界值校正

按表12-15的Lower-Bound自由度校正方法涧尿，的校正自由度系奉，查F界值表，的校正界值為（校正前為）姑廉。的校正自由度缺亮，查F界值表，的校正界值為（校正前為）桥言，校正后的自由度和P值見下表：

image

5.結(jié)論

不同麻醉誘導(dǎo)方法存在組間差別萌踱，見表12-18：

image

患者的收縮壓在不同的誘導(dǎo)方法下，不同誘導(dǎo)時相變化的趨勢不同号阿，見表12-19:

image

其中A組不同誘導(dǎo)時間收縮壓相對穩(wěn)定虫蝶，見表12-20，如下所示：

image

重復(fù)測量數(shù)據(jù)的多重比較

重復(fù)測量數(shù)據(jù)的多重比較是一個比較復(fù)雜的問題倦西，對不同研究有不同的要求，例如前面提到的赁严，A扰柠、B粉铐、C三種麻醉誘導(dǎo)方法存在差別，每組患者不同誘導(dǎo)時相收縮壓的趨勢不同卤档，如下所示：

image

多重比較可能包括3種情況蝙泼，一是A、B劝枣、C三組收縮壓均數(shù)差別的兩兩比較汤踏，需要重復(fù)檢驗3次（A與B，A與C舔腾，B與C）溪胶；二是各組T0，T1稳诚，T2哗脖，T3和T4這5個不同時相收縮壓均數(shù)變化曲線是否有固定趨勢（)其中，是收縮壓均數(shù)扳还，如果拒絕H0才避，對于5次重復(fù)測量結(jié)果，曲線經(jīng)下次多項式(polynomial)變換氨距，需要檢驗下次多項式1次(linear)桑逝、2次(quadratic)、3次(cubic)俏让、4次(order 4)是否為0楞遏，上圖中的三條曲線共需要重復(fù)檢驗12次；三是各組T0舆驶，T1橱健，T2，T3沙廉，T4這5個不同時相收縮壓均數(shù)的兩兩比較拘荡，每組需要比較10次，3組共要做30次重復(fù)檢驗 撬陵。因此珊皿，前面提到的例題的多重檢比較需要做2+12+30=44次。下面我們列出了SPSS重復(fù)測量數(shù)據(jù)的多重比較的計算結(jié)果巨税。

組間差別多重比較

在表12-18基礎(chǔ)上進行兩兩比較蟋定，誤差均方為，由表12-20可知草添，A驶兜、B、C三組收縮壓均數(shù)分別為119.68、124.48抄淑、128.20屠凶，兩兩比較可以選擇LSD、SNK肆资、Bonferroni等多重比較方法矗愧，下圖為LSD兩兩比較結(jié)果，其中A與C的差別有統(tǒng)計學(xué)意義郑原，如下所示：

image

時間趨勢比較

由表12-20可知唉韭，A組5個時相的收縮壓分別為121.00，112.40犯犁，118.40属愤，125.80，120.80栖秕，在SPSS中選擇Polynomial正交多項式春塌。表12-23就列出了A組收縮壓均數(shù)5個時相的時間趨勢檢驗結(jié)果，正交多項式1次項簇捍、3次項有統(tǒng)計學(xué)意義(P<0.05)只壳，拒絕。B組和C組的時間趨勢檢驗結(jié)果見表12-24暑塑、表12-25吼句，正交多項式1次項、2次項事格、4次項均有統(tǒng)計學(xué)意義惕艳，拒絕。如下所示：

image

時間點的多重比較

一些醫(yī)學(xué)研究關(guān)心同一個處理重復(fù)測量數(shù)據(jù)時間點的兩兩差別驹愚，如例12-3中A組5個時相收縮壓均數(shù)之間的差別远搪。由于當(dāng)時間點較多時，重復(fù)檢驗次數(shù)將難以承受逢捺，如果假定5個處理組谁鳍，10個時間點，重復(fù)檢驗次數(shù)就是5x45=225次劫瞳。因此倘潜，時間點多重比較通常采用事前檢驗(prior tests)的策略。事前檢驗與事后檢驗(ost hoc tests)的不同之處在于志于，事前檢驗的檢驗次數(shù)是在試驗開始前就確定了涮因，無論方差分析結(jié)論如何都執(zhí)行重復(fù)檢驗，比較的時間點按專業(yè)意義設(shè)定伺绽，如試驗后各時間點均數(shù)與基線比較养泡，重復(fù)檢驗次數(shù)比較少嗜湃。事后檢驗則是方差分析結(jié)果拒絕才進行兩兩比較，通常所有的時間上熾都要比較瓤荔、重復(fù)檢驗次數(shù)較多净蚤。

按照事前檢驗策略，例12-3中A組5個時間收縮壓均數(shù)差別的兩兩比較输硝，只比較麻醉誘導(dǎo)后各時相T1，T2程梦，T3点把，T4與誘導(dǎo)前T0的差別，只需要檢驗4次屿附，采用Bonferroni多重比較方法郎逃，設(shè)定α=0.05，檢驗的名義水準為挺份，下表為SPSS的檢驗結(jié)果褒翰，可以發(fā)現(xiàn)，A組麻醉誘導(dǎo)后的4個時相與誘導(dǎo)前比較的重復(fù)檢驗匀泊，僅T1收縮壓降低(P<0.05)优训。B組誘導(dǎo)后T2收縮壓降低，T3各聘、T4收縮壓增高(P<0.05)揣非，C組誘導(dǎo)后T2收縮壓降低，T3收縮壓升高(P<0.05)躲因，因此A組麻醉誘導(dǎo)后收縮壓的波動小于B組和C組早敬，如下所示：

image

例11-1：2×2析因設(shè)計

將20只家兔隨機等分4組，每組5只大脉，進行神經(jīng)損傷后的縫合試驗搞监。處理由兩個因素組合而成，A因素為縫合方法镰矿，有兩水平琐驴，一水平為外膜縫合，另一水平為束膜縫合衡怀，B因素為縫合后的時間棍矛，有兩水平，一水平為縫合后1月抛杨，另一水平為縫合后2月够委。試驗結(jié)果為家兔神經(jīng)縫合后的軸突通過率(%)（注：測量指標，視為計量資料）怖现，見表11-1茁帽。欲比較不同縫合方法及縫合后時間對軸突通過率的影響玉罐，試做析因方差分析。

data1101 <- c(10,10,40,50,10,30,30,70,60,30,10,20,30,50,30,50,50,70,60,30)
factor_a <- gl(2,10) # a取兩個水平潘拨，每個水平數(shù)量是10個吊输，排列順序是1...2...
factor_b <- rep(gl(2,5),2) # b取兩個水平，每個水平是10個铁追，排列順序是1..2..1..2
data1101.aov <- aov(data1101~factor_a+factor_b+factor_a*factor_b) # factor_a*factor_b表示a與b的交互作用
summary(data1101.aov)

從上述分析結(jié)果可以看出季蚂，只有B因素的主效應(yīng)達到了p<0.05，而A因素的主效應(yīng)p>0.05琅束，A與B的交互作用的p>0.05扭屁，結(jié)論為：沿不能認為兩咱縫合方法對神經(jīng)軸突通過率有影響，但縫合后2個月與1個月相比涩禀，神經(jīng)軸突通過率提高了料滥。

例11-2：I×J兩因素析因分析

觀察A,B兩種鎮(zhèn)痛藥物聯(lián)合運用在產(chǎn)婦分娩時的鎮(zhèn)痛效果。A藥取3個劑量：1.0mg艾船，2.5mg葵腹，5.0mg；B藥也取3個劑量：5 屿岂，15 践宴，30 。共9個處理組雁社。將27名產(chǎn)婦隨機等分為9組浴井，每組3名產(chǎn)婦，記錄每名產(chǎn)婦分娩時的鎮(zhèn)痛時間霉撵，結(jié)果見表11-7磺浙。試分析A,B兩藥聯(lián)合運用的鎮(zhèn)痛效果。

data1102 <- c(105,80,65,75,115,80,85,120,125,115,105,80,125,130,90,65,120,100,75,95,85,135,120,150,180,190,160)
drug_a <- factor(rep(gl(3,3),3))
drug_b <- factor(gl(3,9))
data1102.aov <- aov(data1102~drug_a+drug_b+drug_a*drug_b)
summary(data1102.aov)

例11-3：I×J×K三因素析因分析

用5×2×2析因設(shè)計研究5種類型的軍裝在兩種環(huán)境徒坡、兩種活動狀態(tài)下的散熱效果撕氧，將100名受試者隨機等分20組，觀察指標是受試者的主觀熱感覺(從“冷”到“熱”按等級評分)喇完，結(jié)果見表11-11伦泥。試進行方差分析。

data1103 <- c(0.25,0.3,0.75,0.2,-0.1,-0.25,0.1,-0.5,-1,0,1.25,0.5,0.6,0.85,2.5,-0.75,-0.35,0.4,-0.5,0.1,0.4,0.05,-0.2,0.9,-0.1,4.75,4.6,4.55,4.25,4.72,3.45,4.8,3.5,3.1,4.3,4,4,4.25,4,4.1,4.85,5.2,4.1,5,4.8,4.55,4.3,4.4,4.2,3.6,0.5,1.5,0.75,-0.75,1.75,2.1,1.5,2.65,0.9,2.4,2.75,1.25,3,0.95,1.75,1,1.37,0.05,0.62,3.05,2.35,2.55,1.17,1.05,2.75,3.75,4,4.1,3.27,4.8,4,4.05,5,4.25,4.02,4,4.15,4.2,4,4.15,4.25,4.1,4.15,4.25,4.75,4.6,4.25,4.17,4.25,4.8)

a <- factor(rep(1:5,20))
b <- factor(rep(1:2,each=50))
c <- factor(rep(c(1,2,1,2),each=25))

data1103.aov <- aov(data1103~a+b+c+a*b+a*c+b*c+a*b*c)
summary(data1103.aov)

例11-4

研究雌螺產(chǎn)卵的最優(yōu)條件锦溪，在20cm2的泥盒里飼養(yǎng)同齡雌螺10只不脯，試驗條件有4個因素（見表11-15），每個因素2個水平刻诊。試在考慮溫度與含氧量對雌螺產(chǎn)卵有交互作用的情況下安排正交試驗防楷。

data1104 <- c(86, 95,91,94,91,96,83,88)
A <- factor(rep(1:2,each=4))
B <- factor(gl(2,2,8))
AB <- factor(c(1,1,2,2,2,2,1,1))
C <- factor(rep(1:2,4))
D <- factor(c(1,2,2,1,2,1,1,2))
data1104.aov <- aov(data1104~A+B+AB+C+D)
summary(data1104.aov)

例11-5

研究高頻呼吸機A，B则涯，C复局，D冲簿，E五個參數(shù)對通氣量的影響，每個參數(shù)有高亿昏、低兩個水平峦剔，試按正交設(shè)計安排試驗。

data1105 <- c(16.26, 19.38,23.6,28.43,20.28,34.88,49.1,47.44,18.32,24.85,39.45,32.08,45.5,50.3,55.26,66.64)
A <- factor(rep(1:2,each=8)) # 產(chǎn)生2個水平的因子角钩，每個因子8個
B <- factor(gl(2,4,16)) # 產(chǎn)生2個水平的因子吝沫，每個因子4個，并且把向量B補充到16個元素
C <- factor(gl(2,2,16))
D <- factor(rep(1:2,8)) # 產(chǎn)生2個水平的因子彤断，2個水平的因子作為整體野舶，復(fù)制8遍
E <- factor(c(1,2,2,1,2,1,1,2,2,1, 1,2,1,2,2,1))

data1105.aov <- aov(data1105~A+B+C+D+E)
summary(data1103.aov)

例11-6

試驗甲、乙宰衙、丙三種催化劑在不同溫度下對某化合物的轉(zhuǎn)化作用。由于各催化劑所要求的溫度范圍不同睹欲，將催化劑作為一級實驗因素（I=3）供炼，溫度作為二級實驗因素（J=3），采用嵌套設(shè)計窘疮，每個處理重復(fù)2次（n=2）袋哼，試驗結(jié)果見表11-25。試做方差分析闸衫。

catalyst <- factor(gl(3,6,18))
temperature <- factor(rep(c(70,80,90,55,65,75,90,95,100),each=2))
data1106 <- c(82, 84,91,88,85,83,65,61,62,59,56,60,71,67,75,78,85,89)

data1106.aov <- aov(data1106~catalyst+temperature)
summary(data1106.aov)

例11-7

試驗一種全身注射抗毒素對皮膚損傷的保護作用試驗涛贯，將10只家兔隨機等分兩組，一組注射抗毒素蔚出，一組注射生理鹽水作對照弟翘。分組后，每只家兔取甲骄酗、乙兩部位稀余，分別隨機分配注射低濃度毒素和高濃度毒素，觀察指標為皮膚受損直徑(mm)趋翻，試驗結(jié)果見表11-31睛琳。試做方差分析。

A <- factor(gl(2,5,20))
B <- factor(rep(1:2,each=10))
data1107 <- c(15.75,15.5,15.5,17,16.5,18.25,18.5,19.75,21.5,20.75,19,20.75,18.5,20.5,20,22.25,21.5,23.5,24.75,23.75)
data1107.aov <- aov(data1107~A+B)
summary(data1107.aov)

參考資料

1. 白話統(tǒng)計.馮國雙
2. 醫(yī)學(xué)統(tǒng)計學(xué).孫振球.第四版
3. 數(shù)據(jù)分析：R語言實戰(zhàn).李詩羽.電子工業(yè)出版社