多重比较中LSD/HSD的计算方法

agricolae包中的相关函数可以实现大多数多重比较的分组标记。

LSD

LSD方法采用了t检验的方式。

#
model<-aov(Sepal.Length~Species, data=iris)
summary(model) # 方差分析模型

#LSD.test结果
out <- LSD.test(model,"Species", p.adj="none",console=T)
> out$statistics$LSD
0.2034688 

#self calculation
#MSE=sum of residuals /Df
#LSD=t(1-alpha/2,DF)*sqrt(MSE*2/n)
> LSD=qt(1-0.05/2,model$df.residual)*
  sqrt(sum(model$residual^2)/
         model$df.residual*2/(
           length(model$residuals)/length(model$coefficients)))
> LSD
0.2034688 

LSD校正

严格校正的如bonferroni，即控制alpha’=alpha/k, k为比较次数。

out <- LSD.test(model,"Species", p.adj="bonferroni",console=F)
#alpha=alpha/k $k is comparison times
> out$statistics
    MSerror  Df     Mean       CV  t.value       MSD
  0.2650082 147 5.843333 8.809859 2.421686 0.2493317

# 手动计算
MSD=qt(1-0.05/2/choose(model$rank,2),model$df.residual)*
  sqrt(sum(model$residual^2)/
         model$df.residual*2/(
           length(model$residuals)/length(model$coefficients)))
> MSD
[1] 0.2493317

Turkey

Turkey又称HSD,其服从q分布，即students range。
HSD=q(1-alpha,k,N-k)*sqrt(MSE/n)
k为组数，N为样本总数,n为每组样本数

out <- HSD.test(model,"Species",console=F)
> out$statistics
   MSerror  Df     Mean       CV       MSD
  0.2650082 147 5.843333 8.809859 0.2437727

# in HSD, the difference is q distribution
# HSD=q(1-alpha,k,Df)$sqrt(MSE/n)
MSD=qtukey(1-0.05,model$rank,model$df.residual)*
  sqrt(sum(model$residual^2)/model$df.residual/
         (length(model$residuals)/length(model$coefficients)))
> MSD
[1] 0.2437727

从上面结果的LSD（0.20）<HSD（0.24）<HSD（0.25）看, 不校正的LSD的条件最宽松，即最敏感，假阳性概率较大，而bonferroni校正的LSD又最保守,HSD较为适宜。