# ————————————————————————————
# Tukey HSD 法による多重比較（Rcmdrでのまちがった計算）
# ————————————————————————————
library(Rcmdr)
## Loading required package: splines
## Loading required package: RcmdrMisc
## Loading required package: car
## Loading required package: sandwich
## Loading required package: effects
## Loading required package: carData
##
## Attaching package: 'carData'
## The following objects are masked from 'package:car':
##
## Guyer, UN, Vocab
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## コマンダーの GUI はインタラクティブなセッションでしか起動しません．
# 二要因乱塊法データ「Box3_R.txt」を読み込む
Dataset <- read.table("/Users/minaka/TeachingMaterials/R-statistics/Mac/R/data/Box3_R.txt", header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
Dataset <- within(Dataset, { REP <- as.factor(REP)
})
# 二要因乱塊法モデルのフィッティング
LinearModel.1 <- lm(DATA ~ REP + N + V +N:V, data=Dataset)
# ——————————————————————————
# Tukey HSD検定（モデルのまちがいに注意）
# ——————————————————————————
library(mvtnorm, pos=17)
library(survival, pos=17)
library(MASS, pos=17)
library(TH.data, pos=17)
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
library(multcomp, pos=17)
library(abind, pos=22)
AnovaModel.2 <- aov(DATA ~ N, data=Dataset)
# ↑このモデル「DATA ~ N」はまちがっている！
with(Dataset, numSummary(DATA, groups=N, statistics=c("mean", "sd")))
## mean sd data:n
## N0 3.483333 0.5139808 12
## N1 4.761000 0.3313680 12
## N2 5.071833 0.6897728 12
## N3 5.669667 0.4411381 12
## N4 5.796667 0.2710335 12
local({
.Pairs <- glht(AnovaModel.2, linfct = mcp(N = "Tukey"))
print(summary(.Pairs)) # pairwise tests
print(confint(.Pairs)) # confidence intervals
print(cld(.Pairs)) # compact letter display
old.oma <- par(oma=c(0,5,0,0))
plot(confint(.Pairs))
par(old.oma)
})
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DATA ~ N, data = Dataset)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## N1 - N0 == 0 1.2777 0.1930 6.619 < 1e-04 ***
## N2 - N0 == 0 1.5885 0.1930 8.229 < 1e-04 ***
## N3 - N0 == 0 2.1863 0.1930 11.326 < 1e-04 ***
## N4 - N0 == 0 2.3133 0.1930 11.984 < 1e-04 ***
## N2 - N1 == 0 0.3108 0.1930 1.610 0.497395
## N3 - N1 == 0 0.9087 0.1930 4.707 0.000157 ***
## N4 - N1 == 0 1.0357 0.1930 5.365 < 1e-04 ***
## N3 - N2 == 0 0.5978 0.1930 3.097 0.024600 *
## N4 - N2 == 0 0.7248 0.1930 3.755 0.003744 **
## N4 - N3 == 0 0.1270 0.1930 0.658 0.964415
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DATA ~ N, data = Dataset)
##
## Quantile = 2.8202
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## N1 - N0 == 0 1.27767 0.73328 1.82206
## N2 - N0 == 0 1.58850 1.04411 2.13289
## N3 - N0 == 0 2.18633 1.64194 2.73072
## N4 - N0 == 0 2.31333 1.76894 2.85772
## N2 - N1 == 0 0.31083 -0.23356 0.85522
## N3 - N1 == 0 0.90867 0.36428 1.45306
## N4 - N1 == 0 1.03567 0.49128 1.58006
## N3 - N2 == 0 0.59783 0.05344 1.14222
## N4 - N2 == 0 0.72483 0.18044 1.26922
## N4 - N3 == 0 0.12700 -0.41739 0.67139
##
## N0 N1 N2 N3 N4
## "a" "b" "b" "c" "c"
AnovaModel.3 <- aov(DATA ~ V, data=Dataset)
# ↑このモデル「DATA ~ V」はまちがっている！
with(Dataset, numSummary(DATA, groups=V, statistics=c("mean", "sd")))
## mean sd data:n
## V1 4.7694 1.0548145 20
## V2 5.0420 0.9574259 20
## V3 5.0581 0.8559827 20
local({
.Pairs <- glht(AnovaModel.3, linfct = mcp(V = "Tukey"))
print(summary(.Pairs)) # pairwise tests
print(confint(.Pairs)) # confidence intervals
print(cld(.Pairs)) # compact letter display
old.oma <- par(oma=c(0,5,0,0))
plot(confint(.Pairs))
par(old.oma)
})
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DATA ~ V, data = Dataset)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## V2 - V1 == 0 0.2726 0.3034 0.898 0.644
## V3 - V1 == 0 0.2887 0.3034 0.951 0.610
## V3 - V2 == 0 0.0161 0.3034 0.053 0.998
## (Adjusted p values reported -- single-step method)
##
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DATA ~ V, data = Dataset)
##
## Quantile = 2.4067
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## V2 - V1 == 0 0.2726 -0.4576 1.0028
## V3 - V1 == 0 0.2887 -0.4415 1.0189
## V3 - V2 == 0 0.0161 -0.7141 0.7463
##
## V1 V2 V3
## "a" "a" "a"
AnovaModel.4 <- aov(DATA ~ REP, data=Dataset)
# ↑このモデル「DATA ~ REP」はまちがっている！
with(Dataset, numSummary(DATA, groups=REP, statistics=c("mean", "sd")))
## mean sd data:n
## 1 5.132800 0.9331476 15
## 2 4.912533 0.9578146 15
## 3 5.146800 1.0105777 15
## 4 4.633867 0.9102152 15
local({
.Pairs <- glht(AnovaModel.4, linfct = mcp(REP = "Tukey"))
print(summary(.Pairs)) # pairwise tests
print(confint(.Pairs)) # confidence intervals
print(cld(.Pairs)) # compact letter display
old.oma <- par(oma=c(0,5,0,0))
plot(confint(.Pairs))
par(old.oma)
})
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DATA ~ REP, data = Dataset)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 -0.2203 0.3482 -0.633 0.921
## 3 - 1 == 0 0.0140 0.3482 0.040 1.000
## 4 - 1 == 0 -0.4989 0.3482 -1.433 0.485
## 3 - 2 == 0 0.2343 0.3482 0.673 0.907
## 4 - 2 == 0 -0.2787 0.3482 -0.800 0.854
## 4 - 3 == 0 -0.5129 0.3482 -1.473 0.460
## (Adjusted p values reported -- single-step method)
##
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DATA ~ REP, data = Dataset)
##
## Quantile = 2.6478
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 2 - 1 == 0 -0.2203 -1.1423 0.7018
## 3 - 1 == 0 0.0140 -0.9080 0.9360
## 4 - 1 == 0 -0.4989 -1.4210 0.4231
## 3 - 2 == 0 0.2343 -0.6878 1.1563
## 4 - 2 == 0 -0.2787 -1.2007 0.6434
## 4 - 3 == 0 -0.5129 -1.4350 0.4091
##
## 1 2 3 4
## "a" "a" "a" "a"
# 以上のまちがいは「一元配置分散分析」の枠内でのみ
# TukeyHSD検定が利用できるというRcmdrの制約による．
# 一要因完全無作為化法ならばそのまま使ってよい．
# 上のように複数の要因を含む場合は，まちがったモデルを
# 補正した上で再計算する必要がある（スクリプトを書く）．