多重比較中LSD HSD的計算方法

agricolae包中的相關函式可以實現大多數多重比較的分組標記。

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

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

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