R Commander 起動
データファイル（Box2_R.txt）インポート
> Dataset <-
+ read.table("/Users/minaka/TeachingMaterials/R-statistics/Win/R/data/Box2_R.txt",
+ header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
ブロック要因（REP）の因子変換
> Dataset <- within(Dataset, {
+ REP <- as.factor(REP)
+ })
箱ひげ図（処理要因）
> Boxplot(DATA~TRT, data=Dataset, id=list(method="y"))
箱ひげ図（ブロック要因）
> Boxplot(DATA~REP, data=Dataset, id=list(method="y"))
[1] "5"
乱塊法の線形統計モデル構築と当てはめ
> LinearModel.1 <- lm(DATA ~ REP + TRT, data=Dataset)
> summary(LinearModel.1)
Call:
lm(formula = DATA ~ REP + TRT, data = Dataset)
Residuals:
Min 1Q Median 3Q Max
-415.7 -215.6 21.0 156.8 627.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5047.00 203.62 24.787 1.38e-13 ***
REP2 55.17 191.97 0.287 0.77776
REP3 -184.50 191.97 -0.961 0.35175
REP4 -667.67 191.97 -3.478 0.00337 **
TRTS125 -139.75 235.12 -0.594 0.56111
TRTS150 -144.50 235.12 -0.615 0.54803
TRTS25 276.25 235.12 1.175 0.25833
TRTS50 222.50 235.12 0.946 0.35897
TRTS75 456.50 235.12 1.942 0.07121 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 332.5 on 15 degrees of freedom
Multiple R-squared: 0.6546, Adjusted R-squared: 0.4704
F-statistic: 3.553 on 8 and 15 DF, p-value: 0.01651
分散分析表の作表
> Anova(LinearModel.1, type="II")
Anova Table (Type II tests)
Response: DATA
Sum Sq Df F value Pr(>F)
REP 1944361 3 5.8622 0.007416 **
TRT 1198331 5 2.1678 0.112809
Residuals 1658376 15
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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)
Tukey HSD 検定（処理要因）
> AnovaModel.2 <- aov(DATA ~ REP + TRT, data=Dataset)
> summary(AnovaModel.2)
Df Sum Sq Mean Sq F value Pr(>F)
REP 3 1944361 648120 5.862 0.00742 **
TRT 5 1198331 239666 2.168 0.11281
Residuals 15 1658376 110558
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> with(Dataset, numSummary(DATA, groups=TRT, statistics=c("mean", "sd")))
mean sd data:n
S100 4847.75 321.9900 4
S125 4708.00 188.4958 4
S150 4703.25 497.3941 4
S25 5124.00 320.2093 4
S50 5070.25 736.4185 4
S75 5304.25 411.6515 4
> local({
+ .Pairs <- glht(AnovaModel.2, linfct = mcp(TRT = "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 + TRT, data = Dataset)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
S125 - S100 == 0 -139.75 235.12 -0.594 0.990
S150 - S100 == 0 -144.50 235.12 -0.615 0.988
S25 - S100 == 0 276.25 235.12 1.175 0.842
S50 - S100 == 0 222.50 235.12 0.946 0.928
S75 - S100 == 0 456.50 235.12 1.942 0.416
S150 - S125 == 0 -4.75 235.12 -0.020 1.000
S25 - S125 == 0 416.00 235.12 1.769 0.512
S50 - S125 == 0 362.25 235.12 1.541 0.646
S75 - S125 == 0 596.25 235.12 2.536 0.175
S25 - S150 == 0 420.75 235.12 1.790 0.500
S50 - S150 == 0 367.00 235.12 1.561 0.634
S75 - S150 == 0 601.00 235.12 2.556 0.169
S50 - S25 == 0 -53.75 235.12 -0.229 1.000
S75 - S25 == 0 180.25 235.12 0.767 0.969
S75 - S50 == 0 234.00 235.12 0.995 0.912
(Adjusted p values reported -- single-step method)
Simultaneous Confidence Intervals
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = DATA ~ REP + TRT, data = Dataset)
Quantile = 3.2509
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
S125 - S100 == 0 -139.7500 -904.0961 624.5961
S150 - S100 == 0 -144.5000 -908.8461 619.8461
S25 - S100 == 0 276.2500 -488.0961 1040.5961
S50 - S100 == 0 222.5000 -541.8461 986.8461
S75 - S100 == 0 456.5000 -307.8461 1220.8461
S150 - S125 == 0 -4.7500 -769.0961 759.5961
S25 - S125 == 0 416.0000 -348.3461 1180.3461
S50 - S125 == 0 362.2500 -402.0961 1126.5961
S75 - S125 == 0 596.2500 -168.0961 1360.5961
S25 - S150 == 0 420.7500 -343.5961 1185.0961
S50 - S150 == 0 367.0000 -397.3461 1131.3461
S75 - S150 == 0 601.0000 -163.3461 1365.3461
S50 - S25 == 0 -53.7500 -818.0961 710.5961
S75 - S25 == 0 180.2500 -584.0961 944.5961
S75 - S50 == 0 234.0000 -530.3461 998.3461
S100 S125 S150 S25 S50 S75
"a" "a" "a" "a" "a" "a"
Tukey HSD 検定（ブロック要因）
> AnovaModel.3 <- aov(DATA ~ REP + TRT, data=Dataset)
> summary(AnovaModel.3)
Df Sum Sq Mean Sq F value Pr(>F)
REP 3 1944361 648120 5.862 0.00742 **
TRT 5 1198331 239666 2.168 0.11281
Residuals 15 1658376 110558
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> with(Dataset, numSummary(DATA, groups=REP, statistics=c("mean", "sd")))
mean sd data:n
1 5158.833 192.3168 6
2 5214.000 558.5213 6
3 4974.333 382.6035 6
4 4491.167 275.7248 6
> local({
+ .Pairs <- glht(AnovaModel.3, 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 + TRT, data = Dataset)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
2 - 1 == 0 55.17 191.97 0.287 0.9914
3 - 1 == 0 -184.50 191.97 -0.961 0.7729
4 - 1 == 0 -667.67 191.97 -3.478 0.0161 *
3 - 2 == 0 -239.67 191.97 -1.248 0.6074
4 - 2 == 0 -722.83 191.97 -3.765 0.0089 **
4 - 3 == 0 -483.17 191.97 -2.517 0.0977 .
---
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 ~ REP + TRT, data = Dataset)
Quantile = 2.884
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
2 - 1 == 0 55.1667 -498.4682 608.8016
3 - 1 == 0 -184.5000 -738.1349 369.1349
4 - 1 == 0 -667.6667 -1221.3016 -114.0318
3 - 2 == 0 -239.6667 -793.3016 313.9682
4 - 2 == 0 -722.8333 -1276.4682 -169.1984
4 - 3 == 0 -483.1667 -1036.8016 70.4682
1 2 3 4
"b" "b" "ab" "a"
