R語言與統(tǒng)計(jì)-1：t檢驗(yàn)與秩和檢驗(yàn)
1. 基礎(chǔ)知識(shí)
方差分析適用于多組均數(shù)的比較（在完全隨機(jī)設(shè)計(jì)的實(shí)驗(yàn)中东帅，兩組均數(shù)的t檢驗(yàn)和方差分析是完全等價(jià)的城看。但t檢驗(yàn)只能用于兩組的均數(shù)比較创葡，對(duì)于三組和三組以上的均數(shù)比較骗村，就需要用到方差分析梳毙。）
-
方差分析的目的：比較多組均數(shù)的差異-->分析一個(gè)數(shù)值型變量與一個(gè)(單因素)或多個(gè)(多因素)分類變量是否相關(guān) -
方差分析的結(jié)果：不同組別的均數(shù)差異顯著-->該數(shù)值型變量與分類變量是相關(guān)的 -
方差分析需要滿足的條件：
1. 在控制變量的不同水平下如暖，觀測(cè)變量服從正態(tài)分布
2. 在控制變量的不同水平下橄教，觀測(cè)變量方差齊性
3. 各組數(shù)據(jù)相互獨(dú)立 -
方差分析基本原理：
1. 實(shí)驗(yàn)條件刻伊，即不同的處理造成的差異，稱為組間差異描验。用變量在各組的均值與總均值之偏差平方和的總和表示白嘁，記做SSb，組間自由度dfb膘流。
2. 隨機(jī)誤差絮缅，如測(cè)量誤差造成的差異或個(gè)體間的差異，稱為組內(nèi)差異呼股，用變量在各組的均值與該組內(nèi)變量值之偏差平方和的總和表示耕魄，記做SSw，組間自由度dfw彭谁。
總偏差平方和SSt=SSb+SSw
整個(gè)方差分析的基本步驟如下：
1吸奴、建立檢驗(yàn)假設(shè)；
H0：多個(gè)樣本總體均值相等；
H1：多個(gè)樣本總體均值不相等或不全等则奥。
檢驗(yàn)水準(zhǔn)為0.05考润。
2、計(jì)算檢驗(yàn)統(tǒng)計(jì)量F值读处；
3糊治、確定P值并作出推斷結(jié)果。
2. 正態(tài)性檢驗(yàn)和方差齊性檢驗(yàn)
# 載入演示數(shù)據(jù)集
# 以multcomp包中的cholesterol數(shù)據(jù)集為例
library('multcomp')
data('cholesterol')
head(cholesterol)
# trt response
# 1 1time 3.8612
# 2 1time 10.3868
# 3 1time 5.9059
# 4 1time 3.0609
# 5 1time 7.7204
# 6 1time 2.7139
str(cholesterol)
# 'data.frame': 50 obs. of 2 variables:
# $ trt : Factor w/ 5 levels "1time","2times",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ response: num 3.86 10.39 5.91 3.06 7.72 ...
可以看到這個(gè)數(shù)據(jù)集只有兩個(gè)變量罚舱，其中治療是分類變量（因子型）井辜，有5個(gè)水平。response是數(shù)值型變量馆匿。要對(duì)每種治療所對(duì)應(yīng)的response的均值進(jìn)行比較抑胎，就只能用方差分析而不能用t檢驗(yàn)。
2.1 進(jìn)行正態(tài)性檢驗(yàn)
shapiro.test(cholesterol$response)
# Shapiro-Wilk normality test
# data: cholesterol$response
# W = 0.97722, p-value = 0.4417
符合正態(tài)分布
2.2 進(jìn)行方差齊性檢驗(yàn)（如下四個(gè)函數(shù)都可進(jìn)行多組數(shù)據(jù)的方差齊性檢驗(yàn)渐北，t檢驗(yàn)中提到的var.test函數(shù)只用于兩組數(shù)據(jù)的方差齊性檢驗(yàn)）
- a.
bartlett.test函數(shù)
bartlett.test(response~trt,data = cholesterol)
# Bartlett test of homogeneity of variances
# data: response by trt
# Bartlett's K-squared = 0.57975, df = 4, p-value =
# 0.9653
要比較均值的數(shù)據(jù)寫～左邊阿逃，分組變量寫右邊。p=0.9653赃蛛，方差齊恃锉。
- b.
flinger.test函數(shù)
fligner.test(response~trt,data = cholesterol)
# Fligner-Killeen test of homogeneity of variances
# data: response by trt
# Fligner-Killeen:med chi-squared = 0.74277, df = 4,
# p-value = 0.946
寫法同上，方差齊呕臂。
- c. car包中的
ncvTest函數(shù)
library(car)
ncvTest(lm(response~trt,data = cholesterol)) #傳入的是線性模型
# Non-constant Variance Score Test
# Variance formula: ~ fitted.values
# Chisquare = 0.1338923, Df = 1, p = 0.71443
- d. car包中的
leveneTest函數(shù)
leveneTest(response~trt,data = cholesterol)
# Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
# group 4 0.0755 0.9893
# 45
需要注意的是破托，如果檢驗(yàn)出方差不齊，我們第一步不是立馬選擇進(jìn)行非參數(shù)檢驗(yàn)歧蒋，而是首先要判斷有無異常值存在土砂，因?yàn)楫惓Ｖ祵?duì)方差的影響很大。當(dāng)然谜洽，到這一步才來檢驗(yàn)有無異常值是不符合數(shù)據(jù)分析的流程的萝映，異常值在進(jìn)行數(shù)據(jù)初步處理的時(shí)候就因該被發(fā)現(xiàn)和處理掉。
3 進(jìn)行方差分析
方差分析包括單因素方差分析阐虚，多因素方差分析序臂，協(xié)方差分析，多元方差分析实束，重復(fù)測(cè)量數(shù)據(jù)方差分析奥秆。
3.1單因素方差分析
-
aov函數(shù)
library(multcomp)
data("cholesterol")
head(cholesterol)
fit <- aov(response~trt,data = cholesterol)
summary(fit)
# Df Sum Sq Mean Sq F value Pr(>F)
# trt 4 1351.4 337.8 32.43 9.82e-13 ***
# Residuals 45 468.8 10.4
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##Df是自由度，Sum Sq是總的離差平方和咸灿，Mean Sq是平均的離差平方和构订，F(xiàn) value是檢驗(yàn)統(tǒng)計(jì)量F值
##p值=0.01364，小于0.05避矢，拒絕原假設(shè)鲫咽，接受備擇假設(shè)签赃，response均值有差異
gplots包的plotmeans函數(shù) 對(duì)上述結(jié)果進(jìn)行可視化
plotmeans(response~trt,data = cholesterol)
-
oneway.test函數(shù)
oneway.test(response~trt,data = cholesterol)
# One-way analysis of means (not assuming equal variances)
# data: response and trt
# F = 30.709, num df = 4.00, denom df = 22.46, p-value = 8.227e-09
oneway.test(response~trt,data = cholesterol,var.equal = T) #??var.equal參數(shù)是對(duì)方差是否齊進(jìn)行設(shè)置
# One-way analysis of means
# data: response and trt
# F = 32.433, num df = 4, denom df = 45, p-value = 9.819e-13
fit <- aov(response~trt,data = cholesterol,var.equal = T) #和aov函數(shù)得到的結(jié)果是一樣的
summary(fit)
# Df Sum Sq Mean Sq F value Pr(>F)
# trt 4 1351.4 337.8 32.43 9.82e-13 ***
# Residuals 45 468.8 10.4
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
3.2 兩因素方差分析
使用ToothGrowth數(shù)據(jù)集進(jìn)行演示
data('ToothGrowth')
head(ToothGrowth)
# len supp dose
# 1 4.2 VC 0.5
# 2 11.5 VC 0.5
# 3 7.3 VC 0.5
# 4 5.8 VC 0.5
# 5 6.4 VC 0.5
# 6 10.0 VC 0.5
levels(ToothGrowth$supp)
# [1] "OJ" "VC"
levels(as.factor(ToothGrowth$dose))
# [1] "0.5" "1" "2"
table(as.factor(ToothGrowth$dose))
# 0.5 1 2
# 20 20 20
aov函數(shù)
fangcha <- aov(len~supp+dose,data = ToothGrowth)
summary(fangcha)
# Df Sum Sq Mean Sq F value Pr(>F)
# supp 1 205.4 205.4 11.45 0.0013 **
# dose 1 2224.3 2224.3 123.99 6.31e-16 ***
# Residuals 57 1022.6 17.9
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
不考慮supp和dose之間的交互作用的情況谷异。結(jié)果顯示兩個(gè)因素都對(duì)小鼠牙齒生長(zhǎng)影響顯著分尸。
fangcha <- aov(len~supp*dose,data = ToothGrowth)
summary(fangcha)
# Df Sum Sq Mean Sq F value Pr(>F)
# supp 1 205.4 205.4 12.317 0.000894 ***
# dose 1 2224.3 2224.3 133.415 < 2e-16 ***
# supp:dose 1 88.9 88.9 5.333 0.024631 *
# Residuals 56 933.6 16.7
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
考慮兩個(gè)因素之間的交互作用：將上面的+換成*。結(jié)果顯示兩個(gè)因素都對(duì)小鼠牙齒生長(zhǎng)影響顯著而且兩者間的互相影響也不容忽視歹嘹。
可視化
interaction.plot(ToothGrowth$dose,ToothGrowth$supp,ToothGrowth$len,type = 'b',pch = c(1,10),col=c('red','green'))
3.3 兩兩比較（多組比較）的方差分析
上述結(jié)果已經(jīng)知道了再五組數(shù)據(jù)中的均值不全相等箩绍，下一步想知道哪些相等哪些不等，就要對(duì)這五組進(jìn)行兩兩比較尺上。
思考??：這里的兩兩比較為什么不能用t檢驗(yàn)材蛛？
背景假設(shè)：共有4組，進(jìn)行兩兩比較的話怎抛，則一共需要進(jìn)行6次卑吭，設(shè)α=0.05，6次比較都不發(fā)生假陽性錯(cuò)誤的概率是(1-α)6=0.735
-
TukeyHSD函數(shù)
TukeyHSD(fit) #直接使用前面得到的fit
# Tukey multiple comparisons of means
# 95% family-wise confidence level
# Fit: aov(formula = response ~ trt, data = cholesterol, var.equal = T)
# $trt
# diff lwr upr p adj
# 2times-1time 3.44300 -0.6582817 7.544282 0.1380949
# 4times-1time 6.59281 2.4915283 10.694092 0.0003542
# drugD-1time 9.57920 5.4779183 13.680482 0.0000003
# drugE-1time 15.16555 11.0642683 19.266832 0.0000000
# 4times-2times 3.14981 -0.9514717 7.251092 0.2050382
# drugD-2times 6.13620 2.0349183 10.237482 0.0009611
# drugE-2times 11.72255 7.6212683 15.823832 0.0000000
# drugD-4times 2.98639 -1.1148917 7.087672 0.2512446
# drugE-4times 8.57274 4.4714583 12.674022 0.0000037
# drugE-drugD 5.58635 1.4850683 9.687632 0.0030633
輸出的結(jié)果：從左往右依次是：兩兩比較马绝、兩兩間的差值豆赏、lwr是95%可信線的下限，upr是上限富稻。最后是p值掷邦。
將結(jié)果可視化：
plot(TukeyHSD(fit))
線段中點(diǎn)是均值，兩端是95%置信區(qū)間椭赋，跨過0說明沒有顯著差異抚岗。
3.4 單因素協(xié)方差分析
在進(jìn)行方差分析時(shí)，所有混雜因素統(tǒng)稱為協(xié)變量哪怔。
協(xié)方差分析的前提條件：
1. 協(xié)變量是連續(xù)型變量
2. 協(xié)變量與Y存在線性關(guān)系宣蔚。且在X的不同水平下，斜率大致相同认境，僅有截距不同胚委。
# 使用litter數(shù)據(jù)集進(jìn)行演示
head(litter)
# dose weight gesttime number
# 1 0 28.05 22.5 15
# 2 0 33.33 22.5 14
# 3 0 36.37 22.0 14
# 4 0 35.52 22.0 13
# 5 0 36.77 21.5 15
# 6 0 29.60 23.0 5
檢驗(yàn)dose對(duì)weight的影響。出生時(shí)間gesttime是協(xié)變量元暴。
- Step 1: 協(xié)變量與因變量之間的線性關(guān)系是否斜率相同篷扩？
# 使用aov函數(shù)查看
facn <- aov(weight~gesttime + dose + gesttime:dose,data = litter)
summary(facn)
# Df Sum Sq Mean Sq F value Pr(>F)
# gesttime 1 134.3 134.30 8.289 0.00537 **
# dose 3 137.1 45.71 2.821 0.04556 *
# gesttime:dose 3 81.9 27.29 1.684 0.17889
# Residuals 66 1069.4 16.20
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
aov后面小括號(hào)里寫的順序：結(jié)果變量～協(xié)變量+自變量。如果要看協(xié)變量和自變量之間是否存在交互茉盏，在后面寫+協(xié)變量：自變量鉴未。最后是data=數(shù)據(jù)集。
結(jié)果顯示兩個(gè)變量之間不存在交互效應(yīng)(p=0.17889, >0.05)鸠姨，可以認(rèn)為它們的斜率是相同的铜秆。
- Step2: 因變量的正態(tài)性和方差齊性
library(car)
qqPlot(lm(weight~gesttime+dose, data = litter), simulate=TRUE,
main='QQ Plot',labels=FALSE)
- Step3: 方差分析
fangcha=aov(weight~gesttime+dose,data = litter)
summary(fangcha)
# Df Sum Sq Mean Sq F value Pr(>F)
# gesttime 1 134.3 134.30 8.049 0.00597 **
# dose 3 137.1 45.71 2.739 0.04988 *
# Residuals 69 1151.3 16.69
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- Step4: 使用effects包里的effect函數(shù)去除協(xié)變量的影響
library(effects)
effect('dose',facn)
# NOTE: dose is not a high-order term in the model
# dose effect
# dose
# 0 5 50 500
# 32.30722 28.50625 30.61078 29.29005
3.5 多元方差分析MANOVA
因變量不止一個(gè)，但是需要將它們作為一個(gè)整體同時(shí)進(jìn)行分析讶迁。例如：某種藥物對(duì)患者血紅蛋白濃度连茧，紅細(xì)胞計(jì)數(shù)，外周血細(xì)胞因子水平等多種因素的影響。
一元方差分析不能滿足啸驯，必須使用多元方差分析的情況：
1. 因變量之間存在曖昧關(guān)系（性質(zhì)類似或存在相關(guān)性）
2. 有時(shí)候單個(gè)變量難以反應(yīng)出組間真實(shí)差異
3. 自變量為分類變量
使用manovs()函數(shù)進(jìn)行性多元方差分析
library(MASS)
attach(UScereal)
y <- cbind(calories,fat,sugars)
fit_m <- manova(y~shelf)
summary(fit_m)
# Df Pillai approx F num Df den Df Pr(>F)
# shelf 1 0.19594 4.955 3 61 0.00383
# Residuals 63
# shelf **
# Residuals
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
3.6 重復(fù)測(cè)量方差分析 REPEAT ANOVA
基本特點(diǎn):
1. 同一受試者在不同時(shí)間點(diǎn)上進(jìn)行多次測(cè)量
2. 所測(cè)指標(biāo)大多為連續(xù)型變量
3. 處理因素：g個(gè)水平（g>=1）客扎，每個(gè)水平下有n個(gè)受試者，共計(jì)g*n個(gè)實(shí)驗(yàn)對(duì)象
people <- rep(1:20,each=3)
time <- rep(LETTERS[1:3],20)
value <- sample(60:180,60)
dat <- data.frame(peo=people,time=time,value=value)
DT::datatable(dat)
fit1 <- aov(value~time,data=dat)
summary(fit1)
# Df Sum Sq Mean Sq F value Pr(>F)
# time 2 3810 1905 1.865 0.164
# Residuals 57 58213 1021
#其本質(zhì)就是隨機(jī)區(qū)組方差分析罚斗，只不過時(shí)間是控制變量
4. 拓展：
參考：https://blog.csdn.net/dingming001/article/details/72822270
