R commander 起動
データファイル（Box1_R.txt）インポート
> Dataset <-
+ read.table("/Users/minaka/TeachingMaterials/R-statistics/Win/R/data/Box1_R.txt",
+ header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
ドットチャート描画
> stripchart(DATA ~ TRT, vertical=TRUE, method="stack", ylab="DATA",
+ data=Dataset)
箱ひげ図描画
> Boxplot(DATA~TRT, data=Dataset, id=list(method="y"))
線形モデル構築と当てはめ
> LinearModel.1 <- lm(DATA ~ TRT, data=Dataset)
> summary(LinearModel.1)
Call:
lm(formula = DATA ~ TRT, data = Dataset)
Residuals:
Min 1Q Median 3Q Max
-572.00 -137.25 -19.25 200.25 688.00
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2128.00 153.93 13.825 5.13e-12 ***
TRT[T.Con] -812.00 217.69 -3.730 0.00124 **
TRT[T.DB] -332.00 217.69 -1.525 0.14214
TRT[T.DDT] 423.75 217.69 1.947 0.06508 .
TRT[T.DK] -447.00 217.69 -2.053 0.05269 .
TRT[T.DM1] -1.25 217.69 -0.006 0.99547
TRT[T.DM2] 550.00 217.69 2.527 0.01962 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 307.9 on 21 degrees of freedom
Multiple R-squared: 0.7373, Adjusted R-squared: 0.6623
F-statistic: 9.826 on 6 and 21 DF, p-value: 0.00003329
分散分析表の作表
> Anova(LinearModel.1, type="II")
Anova Table (Type II tests)
Response: DATA
Sum Sq Df F value Pr(>F)
TRT 5587175 6 9.8255 0.00003329 ***
Residuals 1990237 21
---
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)
> library(multcomp, pos=17)
> library(abind, pos=22)
> AnovaModel.2 <- aov(DATA ~ TRT, data=Dataset)
> summary(AnovaModel.2)
Df Sum Sq Mean Sq F value Pr(>F)
TRT 6 5587175 931196 9.826 0.0000333 ***
Residuals 21 1990237 94773
---
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
AZO 2128.00 408.2622 4
Con 1316.00 188.3808 4
DB 1796.00 162.9601 4
DDT 2551.75 193.5792 4
DK 1681.00 254.1679 4
DM1 2126.75 305.9917 4
DM2 2678.00 488.8619 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)
+ })
Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = DATA ~ TRT, data = Dataset)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
Con - AZO == 0 -812.00 217.69 -3.730 0.01806 *
DB - AZO == 0 -332.00 217.69 -1.525 0.72754
DDT - AZO == 0 423.75 217.69 1.947 0.47408
DK - AZO == 0 -447.00 217.69 -2.053 0.41331
DM1 - AZO == 0 -1.25 217.69 -0.006 1.00000
DM2 - AZO == 0 550.00 217.69 2.527 0.20000
DB - Con == 0 480.00 217.69 2.205 0.33385
DDT - Con == 0 1235.75 217.69 5.677 < 0.001 ***
DK - Con == 0 365.00 217.69 1.677 0.63757
DM1 - Con == 0 810.75 217.69 3.724 0.01815 *
DM2 - Con == 0 1362.00 217.69 6.257 < 0.001 ***
DDT - DB == 0 755.75 217.69 3.472 0.03155 *
DK - DB == 0 -115.00 217.69 -0.528 0.99807
DM1 - DB == 0 330.75 217.69 1.519 0.73090
DM2 - DB == 0 882.00 217.69 4.052 0.00859 **
DK - DDT == 0 -870.75 217.69 -4.000 0.00989 **
DM1 - DDT == 0 -425.00 217.69 -1.952 0.47073
DM2 - DDT == 0 126.25 217.69 0.580 0.99677
DM1 - DK == 0 445.75 217.69 2.048 0.41658
DM2 - DK == 0 997.00 217.69 4.580 0.00259 **
DM2 - DM1 == 0 551.25 217.69 2.532 0.19785
---
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 ~ TRT, data = Dataset)
Quantile = 3.2495
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
Con - AZO == 0 -812.0000 -1519.3681 -104.6319
DB - AZO == 0 -332.0000 -1039.3681 375.3681
DDT - AZO == 0 423.7500 -283.6181 1131.1181
DK - AZO == 0 -447.0000 -1154.3681 260.3681
DM1 - AZO == 0 -1.2500 -708.6181 706.1181
DM2 - AZO == 0 550.0000 -157.3681 1257.3681
DB - Con == 0 480.0000 -227.3681 1187.3681
DDT - Con == 0 1235.7500 528.3819 1943.1181
DK - Con == 0 365.0000 -342.3681 1072.3681
DM1 - Con == 0 810.7500 103.3819 1518.1181
DM2 - Con == 0 1362.0000 654.6319 2069.3681
DDT - DB == 0 755.7500 48.3819 1463.1181
DK - DB == 0 -115.0000 -822.3681 592.3681
DM1 - DB == 0 330.7500 -376.6181 1038.1181
DM2 - DB == 0 882.0000 174.6319 1589.3681
DK - DDT == 0 -870.7500 -1578.1181 -163.3819
DM1 - DDT == 0 -425.0000 -1132.3681 282.3681
DM2 - DDT == 0 126.2500 -581.1181 833.6181
DM1 - DK == 0 445.7500 -261.6181 1153.1181
DM2 - DK == 0 997.0000 289.6319 1704.3681
DM2 - DM1 == 0 551.2500 -156.1181 1258.6181
Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
AZO Con DB DDT DK DM1 DM2
"bc" "a" "ab" "c" "ab" "bc" "c"
