ANOVA检验实例：不同药物药效的检测

在本文中将通过R自带的ANOVA检验方法对三种药物A、B、C的镇痛效果进行检验，分析不同药物的药效。

首先检测小鼠用药前平均热痛反应时间。然后分别注射定量的药物给同种小鼠，并在15min、30min、60min检测其用药后平均热痛反应时间。每种药物检测重复10次，所得数据如下：

RAWDATA.png

其中%代表痛阈提高百分率，计算公式为：

痛阈提高百分率计算公式

然后我们将数据整理为以下格式并保存为csv文件：

整理后的数据

接着，我们首先将数据导入R：

drug<-read.csv('drugdata.csv',header=TRUE)#导入R
attach(drug)
time<-factor(time)#对数据格式进行处理将其转化为因子

然后我们初步先了解一下数据情况，计算不同药物在不同时间点的用药后平均热痛反应时间的平均值及方差：

#calculate mean and sd of each group (drugA,B,C)
#DrugA
A_mean<-tapply(A,time,mean)
A_mean
A_sd<-tapply(A,time,sd)
A_sd
#DrugB
B_mean<-tapply(B,time,mean)
B_mean
B_sd<-tapply(B,time,sd)
B_sd
#DrugC
C_mean<-tapply(C,time,mean)
C_mean
C_sd<-tapply(C,time,sd)
C_sd

结果如下：

平均值和方差计算结果

我们可以用绘图的方式更加直观地了解一下不同药物的药效，有一个初步的判断：

plot(c(0,15,30,60),A_mean,type = 'b',xlab='Time',ylab='Efficacy',main = 'The efficacy of Drugs',pch=16,col='red')
lines(c(0,15,30,60),B_mean,type = 'b',pch=16,col='blue')
lines(c(0,15,30,60),C_mean,type = 'b',pch=16,col='green')
legend('topright',c('A','B','C'),lty=c(1,1,1),pch=c(16,16,16),col=c('red','blue','green'))

image.png

通过上图我们可以初步发现，A的药效最强但持续时间较短，B几乎没有药效，而C虽然药效不如A但是药效更加持久。
但是，我们的判断真的可靠吗？因此我们要进行检验。
进一步，我们检验不同药物的药效，首先我们来剖析一下这个问题：
我们要探究不同药物发挥药效的过程，即在比较不同时间点用药后平均热痛反应时间是否存在差异，这是一个组间比较问题，其零假设为：

零假设

代表不同时间点之间药后平均热痛反应时间不存在差异。
我们可以先用ANOVA方法粗略了解一下：

#analysis drugA
fitA<-aov(A~time)
summary(fitA)
#analysis drugB
fitB<-aov(B~time)
summary(fitB)
#analysis drugC
fitC<-aov(C~time)
summary(fitC)

结果如下：

ANOVA的结果

可以发现药物A、C确实起了药效，其p-value小于0.05。而药物B的p-value大于0.05，说明其各时间点并无差别。因此，证实了我们的判断药物B没有作用。
现在，我们知道了药物A、C起了药效，但是不同时间点是否存在差异呢？
aov()方法可以告诉我们0min、15min、30min、60min之间是否存在差异，但是不能具体哪两者之间存在差异，此时我们进一步用TukeyHSD()的方法进行分析。

TukeyHSD(fitA)
TukeyHSD(fitC)

结果如下：

TukeyHSD of DrugA

我们可以发现15min与0min、30min与0min之间存在明显差异说明药效在15min时已经产生，30min仍存在药效。但是15min与30min之间也存在明显的差异，根据均值，我们可以判断30min药效已经有一定的衰减，而60min与0min没有差异，可以判断60min时药效已经基本消退。

TukeyHSD of DrugC

而DrugC的分析可能有点难以判断。首先，15min与0min存在差异，说明15min已经产生药效，但是30min与0min的p-value为0.0559恰好大于0.05，我们是否要认为其已经失去药效了呢？
根据之前绘制的图片，结合实际情况，应该认为此时仍有药效，毕竟0.0559与0.05的差异并不大。

进一步我们发现30min与15min之间p-value远远大于0.05，因此我们认为15min与30min之间没有差异，即药效没有发生衰退。
而进一步我们发现60min与30min、60min与0min之间的p-value都大于0.05，那么根据p-value的大小及折线图，我们认为药效发生了消退，在60min基本完全失效。

根据上述分析，我们可以得出以下3个结论：
1、A的药效最强但持续时间较短，15min后就已经发生衰退，60min已经基本失效。
2、B没有药效
3、C虽然药效不如A但是药效更加持久，15min时已经产生药效，30min时能保持较高的药效，而60min也基本失效。

上述分析看起来很顺，但是却忽略了重复检测，也就是15min、30min、60min都是针对同一只小鼠进行检测。
因此我们在使用aov的时候可以加入组内误差进行。具体地：
首先要给数据添加一列小鼠编号并命名为subject，然后使用aov

aov(A~time+Error(subject/time))

不过由于TukeyHSD的方法并不支持带有error的anova检验，所以此时我们无法使用TukeyHSD