置換檢驗(yàn)是區(qū)別于參數(shù)檢驗(yàn)進(jìn)行t檢驗(yàn)撞反、卡方檢驗(yàn)歌殃、方差分析撼港，回歸分析（參看前幾節(jié)）的另一種思路方法窖维；相比以前學(xué)過的參數(shù)法盅称，置換檢驗(yàn)更適合處理非正態(tài)數(shù)據(jù)虎忌，存在離群點(diǎn)佑惠，樣本很小嘹履，或者無法做參數(shù)檢驗(yàn)等情況宙址，主要用于生成零假設(shè)的p值轴脐，即回答“效應(yīng)是否存在”這樣的問題；
自助法則主要針對(duì)潛在分布未知抡砂、出現(xiàn)離群點(diǎn)大咱、樣本量過小，或者沒有可供選擇的參數(shù)方法下注益，用于生成置信區(qū)間碴巾。比如像估計(jì)樣本中位數(shù)的置信區(qū)間或者是兩樣本的中位數(shù)之差，而這些在正態(tài)分布理論沒有簡單公式理論套用丑搔。
置換檢驗(yàn)
原理參考文章厦瓢，主要思想我認(rèn)為是求出所有分布的可能（中間的一般為零假設(shè)），出現(xiàn)這種分布的概率啤月。
1煮仇、安裝R包
（1）coin包為針對(duì)獨(dú)立性問題的置換法檢驗(yàn)包
install.packages("coin")
（2）lmPerm包為針對(duì)方差分析與回歸分析的置換法檢驗(yàn)包
- 因?yàn)樽髡咭讶ナ溃创a被CRAN歸檔谎仲，安裝比較麻煩浙垫。
首先在 https://cran.r-project.org/src/contrib/Archive/lmPerm/ 下載lmPerm_1.1-2.tar.gz包文件；然后在https://cran.r-project.org/bin/windows/Rtools/ 下載并安裝RTOOLs（Mac與Linux系統(tǒng)可以跳過此步）郑诺；最后打開Rstudio夹姥，執(zhí)行下述代碼，選擇lmPerm_1.1-2.tar.gz進(jìn)行安裝即可间景。
install.packages(file.choose(),repos = NULL,type="source")
2佃声、coin 包的置換檢驗(yàn)
2.1、t 檢驗(yàn)--兩組間差異
score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)
treatment <- factor(c(rep("A",5), rep("B",5))) #一定要轉(zhuǎn)化為因子
mydata <- data.frame(treatment, score) #要求是數(shù)據(jù)框格式
#比較A,B兩組數(shù)據(jù)是否存在差異
- 一般參數(shù)法
t.test(score~treatment, data=mydata, var.equal=TRUE)
#var.equal=TRUE 參數(shù)假定方差相等倘要，并使用合并方差估計(jì)
#返回結(jié)果p值有顯著意義
- 置換法
library(coin)
oneway_test(score~treatment, data=mydata, distribution="exact")
#返回結(jié)果p值無顯著意義圾亏，教材作者認(rèn)為由于觀測(cè)數(shù)據(jù)少十拣，還是更傾向于oneway_test置換檢驗(yàn)法
distribution=參數(shù)可為exact（精確模式，即依據(jù)所有可能的排列組合志鹃，僅適用于兩樣本問題）夭问、approxiamate(nresample=#)（蒙特卡洛抽樣，#指需要重復(fù)的次數(shù)）曹铃、asymptotic（漸進(jìn)分布抽樣）
2.2缰趋、單因素方差分析--多組間差異
lmPerm包更擅長方差分析。示例實(shí)驗(yàn)設(shè)計(jì)仍為5組接受不同治療方法的多組結(jié)果比較陕见。
library(multcomp)
- 參數(shù)法
attach(cholesterol)
fit <- aov(response ~ trt)
summary(fit)
#返回結(jié)果p值有顯著意義
- 置換法
set.seed(1234)
#因?yàn)橛?0個(gè)樣本秘血，不適合采用精確模式，本例使用了蒙特卡洛法抽樣评甜，為保證分享他人結(jié)果重現(xiàn)性灰粮，所以設(shè)置了隨機(jī)數(shù)種子。
oneway_test(response~trt, data=cholesterol,
distribution=approximate(nresample=9999))
#返回結(jié)果p值有顯著意義忍坷，但尷尬的是與教材結(jié)果不一致粘舟，與參數(shù)法的p值也有一定的差距。
值得注意的是定義重復(fù)次數(shù)的參數(shù)佩研，教材中給的是B柑肴，可能由于版本更新，現(xiàn)在由nresample代替旬薯。
2.3晰骑、卡方檢驗(yàn)--類別變量獨(dú)立性檢驗(yàn)
實(shí)驗(yàn)示例仍為關(guān)節(jié)炎的治療（兩種）與效果（無、部分袍暴、顯著）間的關(guān)系
- 參數(shù)法
library(vcd)
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
chisq.test(mytable)
#返回結(jié)果p值有顯著意義
- 置換法
Arthritis <- transform(Arthritis,
Improved = as.factor(as.numeric(Improved)))
#Improved本身是一個(gè)有序因子些侍，上述操作將其變成一個(gè)分類因子隶症。
#若是有序因子政模，coin() 將會(huì)生成一個(gè)線性與線性趨勢(shì)檢驗(yàn)，而不是卡方檢驗(yàn)
set.seed(1234)
chisq_test(Treatment~Improved, data=Arthritis,
distribution=approximate(nresample=9999))
#返回結(jié)果p值有顯著意義蚂会，且與參數(shù)法較接近淋样。
-
若不轉(zhuǎn)換因子的趨勢(shì)檢驗(yàn)結(jié)果如下：
置換趨勢(shì)檢驗(yàn)
2.4、相關(guān)性檢驗(yàn)
實(shí)驗(yàn)示例為研究文盲率與謀殺率是否相關(guān)
states <- as.data.frame(state.x77)
- 參數(shù)法
attach(states)
cor.test(Illiteracy,Murder)
detach(states)
#返回結(jié)果p值有顯著意義胁住，并且計(jì)算了相關(guān)系數(shù)
- 置換法
set.seed(1234)
spearman_test(Illiteracy ~ Murder, data=states,
distribution=approximate(nresample=9999))
#僅計(jì)算p值趁猴，返回結(jié)果p值有顯著意義(不計(jì)算相關(guān)系數(shù))
3、lmPerm 包的置換檢驗(yàn)
主要為lmp()彪见、aovp()兩個(gè)函數(shù)分別對(duì)應(yīng)參數(shù)法的lm()線性回歸儡司、aov()方差分析。主要格式上的區(qū)別是添加了perm= 參數(shù)余指〔度可以為Exact（精確模式）跷坝、Prob（不斷從可能的序列中抽樣，直至估計(jì)的標(biāo)準(zhǔn)差在估計(jì)的p值0.1之下）碉碉、SPR（使用貫序概率比檢驗(yàn)來判斷何時(shí)停止抽樣）柴钻。值得注意的是當(dāng)樣本觀測(cè)大于10，perm="Exact"自動(dòng)默認(rèn)轉(zhuǎn)為"Prob"垢粮，因此精確檢驗(yàn)只適用于小樣本問題贴届。
3.1、回歸分析 lmp
（1）簡單線性回歸
實(shí)驗(yàn)示例仍為以身高預(yù)測(cè)體重的設(shè)計(jì)
- 參數(shù)法
fit <- lm(weight ~ height, data=women)
summary(fit)
#返回回歸系數(shù)結(jié)果p值有顯著意義
- 置換法
library(lmPerm)
set.seed(1234)
fit1 <- lmp(weight ~ height, data=women, perm="Prob")
summary(fit1)
#返回回歸系數(shù)結(jié)果p值有顯著意義蜡吧，好像不顯示截距項(xiàng)值
（2）多項(xiàng)式回歸
高精度擬合身高體重回歸關(guān)系
- 參數(shù)法
fit0 <- lm(weight ~ height + I(height^2), data=women)
summary(fit0)
- 置換法
library(lmPerm)
set.seed(1234)
fit <- lmp(weight ~ height + I(height^2), data=women, perm="Prob")
summary(fit)
（3）多元回歸
探究謀殺率與多因素的回歸關(guān)系
- 參數(shù)法
states <- as.data.frame(state.x77)
fit0 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
summary(fit0)
- 置換法
set.seed(1234)
fit <- lmp(Murder ~ Population + Illiteracy+Income+Frost,data=states, perm="Prob")
summary(fit)
3.2毫蚓、方差分析 aovp
（1）單因素方差分析
- 參數(shù)法及coin包oneway_test()函數(shù)方法見上 2.2
- aovp()置換法
library(lmPerm)
library(multcomp)
set.seed(1234)
fit <- aovp(response ~ trt, data=cholesterol, perm="Prob")
anova(fit)
#返回p值具有顯著意義，表明各組的療法效果不全相同
有意思的是教材隨送的參考代碼如下昔善，即aovp被替換成lmp绍些，但做出來的結(jié)果與上面一模一樣，發(fā)現(xiàn)下面兩點(diǎn)也有相同的問題........存疑耀鸦！
fit <- lmp(response ~ trt, data=cholesterol, perm="Prob")
（2）單因素協(xié)方差分析
實(shí)驗(yàn)示例仍為藥物對(duì)剛出生小鼠體重影響柬批，協(xié)變量為懷孕時(shí)間
- 參數(shù)法
library(multcomp)
fit <- aov(weight ~ gesttime + dose, data=litter)
summary(fit)
#結(jié)果表明控制懷孕時(shí)間相同時(shí)，不同藥物對(duì)小鼠體重的影響不同
- 置換法
set.seed(1234)
fit1 <- aovp(weight ~ gesttime + dose, data=litter, perm="Prob")
anova(fit1)
#結(jié)果同上
（3）雙因素方差分析（交互效應(yīng)）
實(shí)驗(yàn)示例：兩種藥物分別在不同劑量下對(duì)小鼠牙齒長度的影響袖订。
- 參數(shù)法
fit <- aov(len ~ supp*dose, data=ToothGrowth)
summary(fit)
- 置換法
set.seed(1234)
fit1 <-aovp(len ~ supp*dose, data=ToothGrowth, perm="Prob")
anova(fit1)
自助法求置信區(qū)間
核心思想是有放回的抽樣多次（1000次）
1氮帐、一般步驟
（1）寫一個(gè)能返回帶研究統(tǒng)計(jì)量的函數(shù)；
（2）確定重復(fù)數(shù)洛姑，使用boot()函數(shù)處理上沐；（一般重復(fù)1000次即可；此外有人認(rèn)為初始樣本大小為20-30即可得到足夠好的結(jié)果）楞艾；
（3）boot.ci()函數(shù)計(jì)算統(tǒng)計(jì)量置信區(qū)間参咙。
在回歸分析中，confint()函數(shù)用于提供回歸系數(shù)的置信區(qū)間（95%）
2硫眯、單個(gè)統(tǒng)計(jì)量自助法
實(shí)驗(yàn)示例：使用mtcar數(shù)據(jù)框蕴侧，采用多元回歸，根據(jù)車重和發(fā)動(dòng)機(jī)排量來預(yù)測(cè)汽車的每加侖行駛的英里數(shù)两入。想獲得95%的R平方值（預(yù)測(cè)變量對(duì)響應(yīng)變量可解釋的方差比）的置信區(qū)間
（1）首先寫函數(shù)
rsq <- function(formula, data, indices) { #固定格式
d <- data[indices,] #固定格式净宵；必須聲明，因?yàn)閎oot() 要用其來選擇樣本
fit <- lm(formula, data=d) #計(jì)算公式
return(summary(fit)$r.square) #返回特定值
}
（2）然后使用boot()函數(shù)
library(boot)
set.seed(1234)
results <- boot(data=mtcars, statistic=rsq,
R=1000, formula=mpg~wt+disp)
#分別為數(shù)據(jù)框裹纳、函數(shù)择葡、重復(fù)次數(shù)，供函數(shù)使用的公式
#可以進(jìn)一步查看results結(jié)果
results$t0 #中心值
results$t #所有枚舉值
print(results)
plot(results) #繪制結(jié)果剃氧，看出自助R方值不呈正態(tài)分布
（3）最后boot.ci()函數(shù)求置信區(qū)間
boot.ci(results, type=c("perc", "bca"))
- type參數(shù)設(shè)定了獲取置信區(qū)間的方法敏储，perc方法展示的是樣本均值；bca將根據(jù)偏差對(duì)區(qū)間做簡單調(diào)整朋鞍。作者認(rèn)為一般bca的結(jié)果更可取已添。
3迫横、多個(gè)統(tǒng)計(jì)量自助法
實(shí)驗(yàn)示例：使用mtcar數(shù)據(jù)框，采用多元回歸酝碳，根據(jù)車總和發(fā)動(dòng)機(jī)排量來預(yù)測(cè)汽車的每加侖行駛的英里數(shù)矾踱。想獲取一個(gè)統(tǒng)計(jì)量向量--三個(gè)回歸系數(shù)（截距項(xiàng)、車總疏哗、發(fā)動(dòng)機(jī)排量）95%的置信區(qū)間呛讲。
#首先還是寫函數(shù)，返回回歸系數(shù)向量（三個(gè)統(tǒng)計(jì)量）
bs <- function(formula, data, indices) {
d <- data[indices,] #注意格式
fit <- lm(formula, data=d)
return(coef(fit))
}
library(boot)
set.seed(1234)
results <- boot(data=mtcars, statistic=bs,
R=1000, formula=mpg~wt+disp)
results$t0 #中心值
results$t #所有枚舉值
print(results)
plot(results, index=2)
#索引index1,2,3分別指截距項(xiàng)返奉、車重贝搁，發(fā)動(dòng)機(jī)排量
boot.ci(results, type="bca", index=2)
boot.ci(results, type="bca", index=3)
本次主要學(xué)習(xí)了區(qū)別于參數(shù)法的實(shí)驗(yàn)分析思路--置換檢驗(yàn)；以及一個(gè)求置信區(qū)間的利器--自助法芽偏。
參考教材《R語言實(shí)戰(zhàn)（第2版）》
