# ————————————————————————————
# 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)
# ↑このモデル「DATA ~ REP + N + V +N:V」を以下のTukeyHSD検定で用いる．
# ———————————————————
# TukeyHSD検定（正しいモデルでの計算）
# ———————————————————
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)
#【1】要因「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(LinearModel.1, 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)
})
## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found
## -- default contrast might be inappropriate
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = DATA ~ REP + N + V + N:V, data = Dataset)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## N1 - N0 == 0 1.5995 0.2750 5.816 < 1e-04 ***
## N2 - N0 == 0 1.3355 0.2750 4.856 0.000151 ***
## N3 - N0 == 0 2.5930 0.2750 9.429 < 1e-04 ***
## N4 - N0 == 0 2.6990 0.2750 9.814 < 1e-04 ***
## N2 - N1 == 0 -0.2640 0.2750 -0.960 0.871313
## N3 - N1 == 0 0.9935 0.2750 3.613 0.006746 **
## N4 - N1 == 0 1.0995 0.2750 3.998 0.002215 **
## N3 - N2 == 0 1.2575 0.2750 4.573 0.000386 ***
## N4 - N2 == 0 1.3635 0.2750 4.958 0.000110 ***
## N4 - N3 == 0 0.1060 0.2750 0.385 0.995141
## ---
## 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: lm(formula = DATA ~ REP + N + V + N:V, data = Dataset)
##
## Quantile = 2.8503
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## N1 - N0 == 0 1.5995 0.8157 2.3833
## N2 - N0 == 0 1.3355 0.5517 2.1193
## N3 - N0 == 0 2.5930 1.8092 3.3768
## N4 - N0 == 0 2.6990 1.9152 3.4828
## N2 - N1 == 0 -0.2640 -1.0478 0.5198
## N3 - N1 == 0 0.9935 0.2097 1.7773
## N4 - N1 == 0 1.0995 0.3157 1.8833
## N3 - N2 == 0 1.2575 0.4737 2.0413
## N4 - N2 == 0 1.3635 0.5797 2.1473
## N4 - N3 == 0 0.1060 -0.6778 0.8898
##
## N0 N1 N2 N3 N4
## "a" "b" "b" "c" "c"
#【2】要因「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(LinearModel.1, 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)
})
## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found
## -- default contrast might be inappropriate
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = DATA ~ REP + N + V + N:V, data = Dataset)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## V2 - V1 == 0 0.424 0.275 1.542 0.2821
## V3 - V1 == 0 0.654 0.275 2.378 0.0562 .
## V3 - V2 == 0 0.230 0.275 0.836 0.6828
## ---
## 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: lm(formula = DATA ~ REP + N + V + N:V, data = Dataset)
##
## Quantile = 2.4297
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## V2 - V1 == 0 0.42400 -0.24419 1.09219
## V3 - V1 == 0 0.65400 -0.01419 1.32219
## V3 - V2 == 0 0.23000 -0.43819 0.89819
##
## V1 V2 V3
## "a" "a" "a"
#【3】要因「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(LinearModel.1, 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: lm(formula = DATA ~ REP + N + V + N:V, data = Dataset)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 -0.2203 0.1420 -1.551 0.41700
## 3 - 1 == 0 0.0140 0.1420 0.099 0.99965
## 4 - 1 == 0 -0.4989 0.1420 -3.513 0.00561 **
## 3 - 2 == 0 0.2343 0.1420 1.650 0.36279
## 4 - 2 == 0 -0.2787 0.1420 -1.962 0.21838
## 4 - 3 == 0 -0.5129 0.1420 -3.612 0.00436 **
## ---
## 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: lm(formula = DATA ~ REP + N + V + N:V, data = Dataset)
##
## Quantile = 2.6756
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 2 - 1 == 0 -0.2203 -0.6002 0.1597
## 3 - 1 == 0 0.0140 -0.3660 0.3940
## 4 - 1 == 0 -0.4989 -0.8789 -0.1190
## 3 - 2 == 0 0.2343 -0.1457 0.6142
## 4 - 2 == 0 -0.2787 -0.6586 0.1013
## 4 - 3 == 0 -0.5129 -0.8929 -0.1330
##
## 1 2 3 4
## "b" "ab" "b" "a"
