方差分析指的是不同變量之間互相影響從而導致結果的變化
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不敏感,證明了在不同貨架的谷物營養成分不同的結論