方差分析指的是不同變量之間互相影響從而導致結果的變化

1.單因素方差分析:

　　案例:50名患者接受降低膽固醇治療的藥物,其中三種治療條件使用藥物相同(20mg一天一次,10mg一天兩次,5mg一天四次),剩下的兩種方式是(drugE和drugD),代表候選藥物

　　　　 哪種藥物治療降低膽固醇的最多?

 1 library(multcomp)
 2 attach(cholesterol)
 3 # 1.各組樣本大小
 4 table(trt)
 5 # 2.各組均值
 6 aggregate(response,by=list(trt),FUN=mean)
 7 # 3.各組標準差
 8 aggregate(response,by=list(trt),FUN=sd)
 
 9 # 4.檢驗組間差異
10 fit <- aov(response ~ trt)
11 summary(fit)
12 library(gplots)
13 # 5.繪制各組均值和置信區間
14 plotmeans(response ~ trt,xlab = ‘Treatment‘,ylab = ‘Response‘,main=‘MeanPlot\nwith 95% CI‘)
15 detach(cholesterol)

　　結論:

　　　　1.均值顯示drugE降低膽固醇最多,1time降低膽固醇最少.

　　　　2.說明不同療法之間的差異很大

　　多重比較藥品和服藥次數

1 library(multcomp)
 
2 par(mar=c(5,4,6,2))
3 tuk <- glht(fit,linfct=mcp(trt=‘Tukey‘))
4 plot(cld(tuk,level=.05),col=‘lightgrey‘)

　　　　結論:每天復用4次和使用drugE的時候治療膽固醇效果最好

　　評估檢驗的假設條件

1 library(car)
2 qqPlot(lm(response ~ trt,data=cholesterol),simulate=T,main=‘Q-Q Plot‘,labels=F)
3 bartlett.test(response ~ trt,data=cholesterol)
 
4 # 檢測離群點
5 outlierTest(fit)

　　結論:數據落在95%置信區間的範圍內,說明數據點滿足正態性假設

　2.單因素協方差分析

　　案例:懷孕的小鼠被分為4各小組,每個小組接受不同劑量的藥物劑量(0.5,50,500)產下小鼠體重為因變量,懷孕時間為協變量

 1 data(litter,package = ‘multcomp‘)
 2 attach(litter)
 3 table(dose)
 4 aggregate(weight,by=list(dose),FUN=mean)
 5 fit2 <- aov(weight ~ gesttime + dose)
 6 summary(fit2)
 7 library(effects)
 8 # 取出協變量計算調整的均值
 9 effect(‘dose‘,fit2)
10 contrast <- rbind(‘no drug vs drug‘ = c(3,-1,-1,-1))
11 summary(glht(fit2,linfct=mcp(dose=contrast)))
12 library(HH)
13 ancovaplot(weight ~ gesttime + dose,data=litter)

　結論:0劑量產仔20個,500劑量產仔17個

　　　 0劑量的體重在32左右,500劑量在30左右

　　 懷孕時間和體重相關

　　 用藥劑量和體重相關

　　結論:小鼠的體重和懷孕時間成正比和劑量成反比

3.雙因素方差分析

　　案例:隨機分配60只豚鼠,分別采用兩種餵食方法(橙汁或者維C),各種餵食方法中含有抗壞血酸3鐘含量(0.5,1,2)

　　　　 每種處理組合都分配10只豚鼠,牙齒長度為因變量

1 attach(ToothGrowth)
2 table(supp,dose)
3 aggregate(len,by=list(supp,dose),FUN=mean)
4 aggregate(len,by=list(supp,dose),FUN=sd)
5 # 將dose轉換為因子變量,這樣就不是一個協變量
6 dose <- factor(dose)
7 fit3 <- aov(len ~ supp*dose)
8 summary(fit3)
9 detach(ToothGrowth)

　　結論:主效應的對豚鼠牙齒影響很大

　　結論:在0.5~1mg的區間中維C的豚鼠的牙齒長度超過使用橙汁的小鼠,在1~2的區間內同理,當超過2mg時,兩者對豚鼠牙齒的影響相同

4.重復測量方差

　　案例:在一定濃度的CO2的環境中比較寒帶植物和非寒帶植物的光合作用率進行比較

 1 CO2$conc <- factor(CO2$conc)
 2 w1b1 <- subset(CO2,Treatment == ‘chilled‘)
 3 fit4 <- aov(uptake ~ conc*Type + Error(Plant/(conc)),w1b1)
 4 summary(fit4)
 5 par(las=2)
 6 par(mar=c(10,4,4,2))
 7 with(w1b1,interaction.plot(conc,Type,uptake,type=‘b‘,col=c(‘red‘,‘blue‘),pch=c(16,18),
 8                            main=‘Interaction plot for plant type and concentration‘))
 9 boxplot(uptake~Type*conc,data=w1b1,col=c(‘gold‘,‘green‘),
10         main = ‘Chilled Quebec and Mississippi Plants‘,
11         ylab="Carbon dioxide uptake rate (umol/m^2 sec)")

　　結論:魁北克的植物比密西西比州的二氧化碳的吸收率高,隨著CO2的濃度體高,效果越明顯

5.多元方差分析

　　案例:研究美國食物中的卡路裏,脂肪,糖分是否會因貨架的不同而不同

1 library(MASS)
2 attach(UScereal)
3 shelf <- factor(shelf)
4  y <- cbind(calories,fat,sugars)
5  aggregate(y,by=list(shelf),FUN=mean)
6  cov(y)
7  fit5 <- manova(y ~ shelf)
8  summary(fit5)
9  summary.aov(fit5)

　　找出離群點

 1 center <- colMeans(y)
 2 n <- nrow(y)
 3 p <- ncol(y)
 4 cov <- cov(y)
 5 d <- mahalanobis(y,center,cov)
 6 coord <- qqplot(qchisq(ppoints(n),df=p),d,
 7                 main="QQ Plot Assessing Multivariate Normality",
 8                 ylab="Mahalanobis D2")
 9 abline(a=0,b=1)
10 identify(coord$x,coord$y,labels = row.names(UScereal))

　　結論:在不同的貨架上的谷物營養成分不同,有兩個產品不符合多元正態分布

1 library(rrcov)
2 # 穩健多元方差分析
3 Wilks.test(y,shelf,method=‘mcd‘)

　　結論:穩健檢測對離群點和違反MANOVA不敏感,證明了在不同貨架的谷物營養成分不同的結論

R語言-方差分析