#方差分析
#小写字母表示定量变量，大写字母表示组别因子
#单因素ANOVA y~A
#含单个协变量的单因素 y~x+A
#双因素ANOVA y~A*B
#含两个协变量的双因素 y~x1+X2+A*B
#随机化区组 y~B+A
#单因素组内ANOVA y~A+Error(Subject/A)
#含组内因子(w)和单个组间因子(B)的重复测量 y~B*W+Error(Subject/W)
#表达式中各项的顺序
#情况一：因子不止一个，并且是非平衡设计
#情况二：存在协变量
#范例一
y~A+B+A:B
#A对于y的影响
#控制A时，B对y的影响
#控制A和B的主效应时，A与B的交互效应
#类型1，序贯型--第一个不调整，后面的顺次依据前n-1做调整
#效应根据表达式中先出现的效应做调整，A不调整，B根据A调整，A:B交互项根据A和B调整
#类型2，分层型--同层的相互参照，交互层的参照单层
#效应根据同水平或低水平的效应做调整。A根据B调整，B根据A调整，A:B交互项同时根据A和B调整
#类型3，边界型--除自己外，与其他各效应同时调整
#每个效应应根据模型其他各效应做响应调整。A根据B和A:B做调整，A:B交互项根据A和B调整
#R默认调用类型1方法，其他软件(比如SAS和SPSS)默认调用类型为类型3方法
#单因素方差分析
#范例一
library("mvtnorm")
library("survival")
library("TH.data")
library("MASS")
library(multcomp)
attach(cholesterol)
table(trt) #统计每个不同水平下，样本数量
aggregate(response,by=list(trt),FUN=mean) #基于trt分组，计算各组Y的均值
aggregate(response,by=list(trt),FUN=sd) #基于trt分组，计算各组Y的标准差
fit <- aov(response~trt) #生成模型
#F检验，Pr(>F)=9.82e-13 ***，落入拒绝域，显著性，说明所有疗法均值不同
#原假设：mean(1time)=mean(2time)=mean(3time)=mean(drugD)=mean(drugE)
summary(fit)
library(gplots)
#plotmeans可以用来绘制带有置信区间的组均值图形
plotmeans(response~trt,xlab="Treatment",ylab="Response",main="Mean Plot\nwith 95% CI")
detach(cholesterol)
#多重比较
#范例一
#2times-1time p adj=0.1380949,未落入拒绝域，均值差异不显著
#4times-1time p adj=0.0003542，落入拒绝域，均值差异显著
TukeyHSD(fit)
p <- par(no.readonly = TRUE)
par(las=2) #全局画图都有影响，同时影响纵，横坐标
par(mar=c(5,8,4,2))
plot(TukeyHSD(fit)) #置信区间包含0的疗法说明差异不显著，跟p值列可以对应起来看
par(p)
#范例二
library(multcomp)
par(mar=c(5,4,6,2))
#glht(模型,linfct=mcp(B=设定的值))
#linfct:a specification of the linear hypotheses(假设) to be tested
#mcp(..., interaction_average = FALSE, covariate_average = FALSE)
#covariate_average: 逻辑指示比较是否在其他协变量上平均,默认值为FALSE
#interaction_average: 逻辑指示比较是否在交互项上平均,默认值为FALSE
#生成多重比较对象
#trt为类别变量，分组变量
#Tukey检验 使用学生化的范围统计量进行组间所有成对比较。Tukey的检验特点：
#所有各组的样本数相等；
#各组样本均数之间的全面比较；
#可能产生较多的假阴性结论。
m <- mcp(trt="Tukey")
class(m) #返回mcp对象
tuk <- glht(fit,linfct=mcp(trt="Tukey"))
#cld函数中的level选项设置了使用的显著水平(0.05,即95的置信区间)
plot(cld(tuk,level=0.05),col="lightgrey")
#用相同的字母(用箱线图表示)说明均值差异不显著
#1time-2times 均标识成a, 差异不显著，未落入拒绝域
#2times-4times 均标识为b, 差异不显著，未落入拒绝域
#4times-drugD 均标识为c, 差异不显著，未落入拒绝域
#评估检验的假设条件
#报错
#Error in rank(x, ties.method = "min", na.last = "keep") :
#unimplemented type 'list' in 'greater'
#注意函数大小写，否则容易引用错误的函数，导致调用失败
#正态性
library(carData)
library(car)
#数据落在95%的置信区间内说明满足正态假设
qqPlot(lm(response~trt,data=cholesterol),
simulate=TRUE,main="Q-Q Plot",labels=FALSE)
#方差齐性检验，等价于方差相等
#p-value = 0.9653,未落入拒绝域，不显著，接受原假设，方差相等
bartlett.test(response~trt,data=cholesterol)
#方差齐性分析对离群点非常敏感
#p-value=0.029422<0.05,落入拒绝域，拒绝原假设：包含离群点，接受备择假设：不包含离群点
library(car)
outlierTest(fit)
#结论：从Q-Q图，Bartlett检验，离群点检验，该数据似乎可以用ANOVA模型拟合的很好
#单因素协方差分析 ANCOVA
#单因素协方差分析扩展了单因素方差分析，包含一个或多个定量的协变量
#data 函数
#Loads specified data sets, or list the available data sets.
data() #list all available data sets
try(data(package = "rpart") ) # list the data sets in the rpart package
data(litter,package="multcomp") #加载到Environment 内存中，变量名：litter
class(litter) #返回"data.frame"
attach(litter)
table(dose) #将litter数据框，基于dose各水平值，统计样本数量
aggregate(weight,by=list(dose),FUN=mean) #基于dose列值分组，计算y列各组下的均值
fit <- aov(weight~gesttime+dose) #生成方差分析模型，gesttime为协变量，dose为分组变量
#F检验表明
#Pr(>F)=0.00597，落入拒绝域，说明怀孕时间与幼崽出生体重相关
#Pr(>F)=0.04988，落入拒绝域，说明当控制怀孕时间，药物剂量和出生体重相关
summary(fit)
detach(litter)
#范例二
#去除协变量效应后的组均值
library(effects)
effect("dose",fit)
#与aggregate(weight,by=list(dose),FUN=mean)相比
#当前值相似，但不是所有时候都这样
#范例三 自定义多重比较
#假定你对未用药条件与其他三种用药条件影响是否不同感兴趣
#-1表示该组为对照组，0或者不为-1表示该组不参与比较；
#dose下分组0,5,50,500 对应3,-1,-1,-1
library(multcomp)
contrast <- rbind("no drug vs. drug"=c(3,-1,-1,-1))
#rownames(contrast) #返回"no drug vs. drug"
g <- glht(fit,linfct=mcp(dose=contrast))
summary(g)
#Pr(>|t|) =0.012，落入拒绝域，t统计量
#得出未用药组比其他用药组条件下的出生体重高的结论
#评估假设检验的条件
#正态性
#同方差性
#回归斜率相同
#ANCOVA模型包含怀孕时间X剂量的交互项时，可对回归斜率的同质性进行检验
library(multcomp)
fit2 <- aov(weight~gesttime*dose,data=litter)
summary(fit2)
#Pr(>F)=0.17889,不显著，未落入拒绝域，接受原假设：斜率相等
#若假设不成立，可以尝试变换协变量或因变量，或使用能对每一个斜率独立解释的模型
#结果可视化
#ancova 绘制因变量，协变量，和因子之间的关系图
install.packages("HH")
library(lattice)
library(grid)
library(latticeExtra)
library(RColorBrewer)
library(gridExtra)
library(HH)
#三个维度，纵轴为因变量，下方横轴为gesttime,上方横轴为dose
#基于类别变量画成一个个长方形块，都包含一个gesttime(up or down),一个Y轴，
#用gesttime来预测weight的回归线，相互平行，仅截距项不同
#基于每一个长方形块来解释，每个分组内部，体重会随着怀孕时间的增长而增长
#回归线，截距项中0组最大，500组最大截距项最小
ancova(weight~gesttime + dose,data=litter)
#包含交互项，生成的图形将允许斜率和截距项依据组别而发生变化
#这对那些违背回归斜率同质性的实例非常有用
ancova(weight~gesttime*dose,data=litter)
#双因素方差分析
#范例一
attach(ToothGrowth)
table(supp,dose) #supp为行维度，dose为列维度，基于当前行维度条件，列维度条件，统计样本量
aggregate(len,by=list(supp,dose),FUN=mean)#基于行维度条件，列维度条件具体分组后，计算因变量的均值
aggregate(len,by=list(supp,dose),FUN=sd)#基于行维度条件，列维度条件具体分组后，计算因变量的标准差
class(dose) #返回 "numeric"
#dose 原始为 "numeric"，aov无法把它识别成分组变量而不是数值型协变量
#factor类型才可以被自动识别为分组类型
dose <- factor(dose)
fit <- aov(len~supp*dose)
#F检验，落入拒绝域，拒绝原假设：Y与每一个预测变量的均值相等，
summary(fit)
detach(ToothGrowth)
#结果的可视化
#方法一
#supp: vc,oj,两种喂食方式，橙汁或者维生素c (对应两色线)
#dose: 0.5,1,2,三种抗坏血酸含量，0.5mg,1mg,2mg
#共享x轴，y轴，supp做为第二个因素画多条线
#Y轴包含aggregate(len,by=list(supp,dose),FUN=mean)的值
interaction.plot(dose,supp,len,type="b",
col=c("red","blue"),pch=c(16,18),
main="Interaction between Dose and Supplement Type")
#方法二
#interaction函数
#gl(n, k, length = n*k, labels = 1:n, ordered = FALSE)
#n 整数的级别（因子）数
#k 一个整数，重复数
#length 一个整数，结果的长度
#labels 可选的向量因子水平的标签
#ordered 一个逻辑表明是否排序
a <- gl(2, 4, 8)
a #返回1 1 1 1 2 2 2 2，Levels: 1 2
class(a) #返回"factor"
b <- gl(2, 2, 8, labels = c("ctrl", "treat"))
b #返回ctrl ctrl treat treat ctrl ctrl treat treat，Levels: ctrl treat
s <- gl(2, 1, 8, labels = c("M", "F"))
s #返回M F M F M F M F，Levels: M F
#[1] 1.ctrl 1.ctrl 1.treat 1.treat 2.ctrl 2.ctrl 2.treat 2.treat
#连接字符串，a[1].b[1],依次类推
x <-interaction(a, b)
class(x) #返回"factor"
#[1] 1:ctrl:M 1:ctrl:F 1:treat:M 1:treat:F 2:ctrl:M 2:ctrl:F 2:treat:M 2:treat:F
#连接字符串，a[1]:b[1]:s[1],依次类推
interaction(a, b, s, sep = ":")
#interaction(supp,dose,sep=" ")
#将每一行的两列数值按照格式：supp dose，改造成分组变量，factor类型，减少实际分组类型
#Levels: OJ 0.5 VC 0.5 OJ 1 VC 1 OJ 2 VC 2 对应横坐标
#connect=list(c(1,3,5),c(2,4,6))，指定将哪些level用线串联
library(gplots)
plotmeans(len~interaction(supp,dose,sep=" "),
connect=list(c(1,3,5),c(2,4,6)),
col=c("red","darkgreen"),
main="Interaction Plot with 95% CIs",
xlab="Treatment and Dose Combination"
)
#方法三
#因变量在y轴
#|前面的表示第一个变量，对应x轴
#|后面的表示第二个变量，对应不同的线
library(HH)
interaction2wt(len~supp*dose)
#重复测量方差分析
#即受试者被测量不止一次
#范例一
#CO2数据集：包含了北方和南方植物的寒冷容忍度研究结果
#uptake 二氧化碳吸收量，因变量
#Treatment 类型：chilled 耐寒，nonchilled 非耐寒
#Type 组间因子，自变量，植物类型：Mississippi(密西西比)，Quebec(魁北克)
#conc 组内因子，二氧化碳浓度
CO2$conc <- factor(CO2$conc) #将数值列因子化，Levels: 95 175 250 350 500 675 1000
w1b1 <- subset(CO2,Treatment=='chilled') #抽取样本子集
#y ~ B*W+Error(Subject/w)
#w 单个组内因子
#B 单个组间因子
#Subject 对被试者独有的标识变量
fit <- aov(uptake~conc*Type+Error(Plant/conc),w1b1)
summary(fit)
par(las=2)
par(mar=c(10,4,4,2))
with(w1b1,
interaction.plot(conc,Type,uptake,type="b",col=c("red","blue"),pch=c(16,18),
main="Interaction Plot for Plant Type and Concentration")
)
boxplot(uptake~Type*conc,data=w1b1,col=(c("gold","green")),
main="Chilled Quebec and Mississippi Plants",
ylab="Carbon dioxide uptake rate (umol/m^2 sec)"
)
#多元方差分析
#当因变量(结果变量)不止一个时，可用多元方差分析(MANOVA)对它们同时进行分析
#范例一 单因素多元方差分析
library(MASS)
attach(UScereal)
shelf <- factor(shelf) #将数值列转换成因子，方便作为分组依据
y <- cbind(calories,fat,sugars) #将因变量提取出来，改造成为一个源数据框
aggregate(y,by=list(shelf),FUN=mean) #基于shelf分组条件，计算组内各因变量平均值
cov(y) #反馈方差，协方差
fit <- manova(y~shelf)#生成模型
#F检验，说明基于shelf分组以后，三个因变量的三组的均值显然不同
summary(fit)
#F检验，分别说明shefl分组以后，每一个因变量的三组均值显然不同
summary.aov(fit)
detach(UScereal)
#评估假设检验
#多元正态性
#范例一
center <- colMeans(y) #针对数据框y，对每一个列求均值
n <- nrow(y) #求行数，返回65
p <- ncol(y) #求列数，返回3
cov <- cov(y) #返回协方差，方差矩阵
#Mahalanobis
#Mahalanobis距离是表示数据的协方差距离。它是一种有效的计算两个未知样本集的相似度的方法
d <- mahalanobis(y,center,cov)
#ppoints(n) #用于概率绘图的纵坐标
x <- ppoints(n)
class(x) #反馈 "numeric"
#qchisq(ppoints(n),df=p) 基于自由度，向量中求取分位数q(quantile function)，做为横坐标值
#d 纵坐标值
coord <- qqplot(qchisq(ppoints(n),df=p),d,main="Q-Q Plot Assessing Multivariate Normality",
ylab = "Mahalanobis D2"
)
abline(a=0,b=1) #画直线，截距=0，斜率=1
identify(coord$x,coord$y,labels=row.names(UScereal)) #互动式，十字架时，点击原点，ESC退出后可以看到文字描述
#方差-协方差矩阵同质性
#R中暂时没有Box‘s M检验
#检验多元离群点
install.packages("mvoutlier")
library("sgeostat")
library(mvoutlier)
outliers <- aq.plot(y)
class(outliers) #返回list
v <- unlist(outliers) #强制转换成向量
grep("TRUE",v,fixed=TRUE)
#Error in C("TRUE", "FALSE", "TRUE") :
#object not interpretable as a factor
#解决方法：c小写，非大写
v <- c("TRUE","FALSE","TRUE")
class(v) #返回"character"
grep("TRUE",v,fixed=TRUE)
#稳健多元方差分析
#如果多元正态性或者方差-协方差均值假设都不满足，或者你担心多元离群点，那么可以考虑用稳健或者非参数版本的MANOVA检验
library(robustbase)
library(rrcov)
Wilks.test(y,shelf,method="mcd")
#用回归来做ANOVA
#ANOVA和线性回归都是广义线性模型的特例
#范例一
library(mvtnorm)
library(survival)
library(multcomp)
levels(cholesterol$trt)
fit.aov <- aov(response~trt,data=cholesterol)
summary(fit.aov)
#分组变量被拆分成全部水平的预测变量
#因为线性模型要求预测变量是数值型
#当lm函数碰到因子时，它会用一系列与因子水平相对应的数值型对照变量来代替因子
#如果有k个水平，将会创建k-1个对照变量
fit.lm <- lm(response~trt,data=cholesterol)
#trt2times，表示1time和2times的一个对照
#trt4times, 表示1time和4times的一个对照
#trtdrugD, 表示1time和drugD的一个对照
#trtdrugE, 表示1time和drugE的一个对照
summary(fit.lm)
#contrasts="contr.helmert",缺省调用
fit.lm <- lm(response~trt,data=cholesterol,contrasts="contr.helmert")
summary(fit.lm)
#内置对照
#contr.helmert 第二个水平对照第一个水平，第三个水平对照前两个的均值，第四个水平对照前三个的均值
#contr.poly 基于正交多项式的对照，用于趋势分析（线性，二次，三次等）和等距水平的有序因子
#contr.sum 对照变量纸盒限制为0，也称做偏差找对，对各水平的均值与所有水平的均值进行比较
#contr.treatment 各水平对照基线水平(默认第一个水平)，也称为虚拟编码
#contr.SAS 类似于contr.treatment,只是基线水平变成了最后一个水平。生成的系数类似于大部分
#SAS过程中使用的对照变量
#范例一
contrasts(cholesterol$trt)
#纵向看，当患者处于drugD,变量drugD=1，其他1time,2times,4times,drugE均为0
#横向看，当其他变量都是0的时候，表示患者目前处于1time条件，所以无需列出第一行条件
#R环境中默认的对照方法
#无序因子的默认对比方法为contr.SAS,有序因子的默认对比方法为contr.helmert
options(contrasts = c("contr.SAS","contr.helmert"))
扫一扫
3368
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
