R Commander 起動

データファイル（Box3_R.txt）のインポート

> Dataset <- 
+   
+   read.table("/Users/minaka/TeachingMaterials/R-statistics/Win/R/data/Box3_R.txt",
+   
+    header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)

データ表示

> Dataset

   REP  V  N  DATA
1    1 V1 N0 3.852
2    1 V1 N1 4.788
3    1 V1 N2 4.576
4    1 V1 N3 6.034
5    1 V1 N4 5.874
6    1 V2 N0 2.846
7    1 V2 N1 4.956
8    1 V2 N2 5.928
9    1 V2 N3 5.664
10   1 V2 N4 5.458
11   1 V3 N0 4.192
12   1 V3 N1 5.250
13   1 V3 N2 5.822
14   1 V3 N3 5.888
15   1 V3 N4 5.864
16   2 V1 N0 2.606
17   2 V1 N1 4.936
18   2 V1 N2 4.454
19   2 V1 N3 5.276
20   2 V1 N4 5.916
21   2 V2 N0 3.794
22   2 V2 N1 5.128
23   2 V2 N2 5.698
24   2 V2 N3 5.362
25   2 V2 N4 5.546
26   2 V3 N0 3.754
27   2 V3 N1 4.582
28   2 V3 N2 4.848
29   2 V3 N3 5.524
30   2 V3 N4 6.264
31   3 V1 N0 3.144
32   3 V1 N1 4.562
33   3 V1 N2 4.884
34   3 V1 N3 5.906
35   3 V1 N4 5.984
36   3 V2 N0 4.108
37   3 V2 N1 4.150
38   3 V2 N2 5.810
39   3 V2 N3 6.458
40   3 V2 N4 5.786
41   3 V3 N0 3.738
42   3 V3 N1 4.896
43   3 V3 N2 5.678
44   3 V3 N3 6.042
45   3 V3 N4 6.056
46   4 V1 N0 2.894
47   4 V1 N1 4.608
48   4 V1 N2 3.924
49   4 V1 N3 5.652
50   4 V1 N4 5.518
51   4 V2 N0 3.444
52   4 V2 N1 4.990
53   4 V2 N2 4.308
54   4 V2 N3 5.474
55   4 V2 N4 5.932
56   4 V3 N0 3.428
57   4 V3 N1 4.286
58   4 V3 N2 4.932
59   4 V3 N3 4.756
60   4 V3 N4 5.362

ブロック要因の因子化

> Dataset <- within(Dataset, {
+   REP <- as.factor(REP)
+ })

箱ひげ図描画

窒素要因（N）

> Boxplot(DATA~N, data=Dataset, id=list(method="y"))

品種要因（V）

> Boxplot(DATA~V, data=Dataset, id=list(method="y"))

ブロック要因（REP）

> Boxplot(DATA~REP, data=Dataset, id=list(method="y"))

[1] "6"  "16"

線形統計モデル構築と当てはめ

> LinearModel.1 <- lm(DATA ~ REP + N + V + N:V, data=Dataset)

> summary(LinearModel.1)


Call:
lm(formula = DATA ~ REP + N + V + N:V, data = Dataset)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.87830 -0.17743  0.00708  0.24239  0.57413 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.3003     0.2130  15.493  < 2e-16 ***
REP2         -0.2203     0.1420  -1.551  0.12840    
REP3          0.0140     0.1420   0.099  0.92194    
REP4         -0.4989     0.1420  -3.513  0.00107 ** 
NN1           1.5995     0.2750   5.816 7.30e-07 ***
NN2           1.3355     0.2750   4.856 1.70e-05 ***
NN3           2.5930     0.2750   9.429 6.30e-12 ***
NN4           2.6990     0.2750   9.814 1.96e-12 ***
VV2           0.4240     0.2750   1.542  0.13063    
VV3           0.6540     0.2750   2.378  0.02203 *  
NN1:VV2      -0.3415     0.3889  -0.878  0.38490    
NN2:VV2       0.5525     0.3889   1.421  0.16281    
NN3:VV2      -0.4015     0.3889  -1.032  0.30782    
NN4:VV2      -0.5665     0.3889  -1.457  0.15266    
NN1:VV3      -0.6240     0.3889  -1.604  0.11611    
NN2:VV3       0.2065     0.3889   0.531  0.59824    
NN3:VV3      -0.8185     0.3889  -2.105  0.04135 *  
NN4:VV3      -0.5905     0.3889  -1.518  0.13643    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3889 on 42 degrees of freedom
Multiple R-squared:  0.8813,    Adjusted R-squared:  0.8333 
F-statistic: 18.35 on 17 and 42 DF,  p-value: 3.375e-14

分散分析表

> Anova(LinearModel.1, type="II")

Anova Table (Type II tests)

Response: DATA
          Sum Sq Df F value  Pr(>F)    
REP        2.600  3  5.7294 0.00222 ** 
N         41.235  4 68.1534 < 2e-16 ***
V          1.053  2  3.4801 0.03995 *  
N:V        2.291  8  1.8931 0.08671 .  
Residuals  6.353 42                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

処理平均のプロット

> with(Dataset, plotMeans(DATA, N, V, error.bars="se", connect=TRUE, 
+   legend.pos="farright"))

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)

窒素要因（N）の多重比較

> AnovaModel.2 <- aov(DATA ~ REP + N + V + N:V, data=Dataset)

> summary(AnovaModel.2)

            Df Sum Sq Mean Sq F value  Pr(>F)    
REP          3   2.60   0.867   5.729 0.00222 ** 
N            4  41.23  10.309  68.153 < 2e-16 ***
V            2   1.05   0.526   3.480 0.03995 *  
N:V          8   2.29   0.286   1.893 0.08671 .  
Residuals   42   6.35   0.151                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> 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)
+ })

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: aov(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.000175 ***
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.871311    
N3 - N1 == 0   0.9935     0.2750   3.613 0.006836 ** 
N4 - N1 == 0   1.0995     0.2750   3.998 0.002217 ** 
N3 - N2 == 0   1.2575     0.2750   4.573 0.000431 ***
N4 - N2 == 0   1.3635     0.2750   4.958 0.000115 ***
N4 - N3 == 0   0.1060     0.2750   0.385 0.995140    
---
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 + N + V + N:V, data = Dataset)

Quantile = 2.8496
95% family-wise confidence level
 

Linear Hypotheses:
             Estimate lwr     upr    
N1 - N0 == 0  1.5995   0.8158  2.3832
N2 - N0 == 0  1.3355   0.5518  2.1192
N3 - N0 == 0  2.5930   1.8093  3.3767
N4 - N0 == 0  2.6990   1.9153  3.4827
N2 - N1 == 0 -0.2640  -1.0477  0.5197
N3 - N1 == 0  0.9935   0.2098  1.7772
N4 - N1 == 0  1.0995   0.3158  1.8832
N3 - N2 == 0  1.2575   0.4738  2.0412
N4 - N2 == 0  1.3635   0.5798  2.1472
N4 - N3 == 0  0.1060  -0.6777  0.8897

 N0  N1  N2  N3  N4 
"a" "b" "b" "c" "c"

品種要因（V）の多重比較

> AnovaModel.3 <- aov(DATA ~ REP + N + V + N:V, data=Dataset)

> summary(AnovaModel.3)

            Df Sum Sq Mean Sq F value  Pr(>F)    
REP          3   2.60   0.867   5.729 0.00222 ** 
N            4  41.23  10.309  68.153 < 2e-16 ***
V            2   1.05   0.526   3.480 0.03995 *  
N:V          8   2.29   0.286   1.893 0.08671 .  
Residuals   42   6.35   0.151                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> 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)
+ })

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: aov(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.0561 .
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: aov(formula = DATA ~ REP + N + V + N:V, data = Dataset)

Quantile = 2.4298
95% family-wise confidence level
 

Linear Hypotheses:
             Estimate lwr     upr    
V2 - V1 == 0  0.4240  -0.2442  1.0922
V3 - V1 == 0  0.6540  -0.0142  1.3222
V3 - V2 == 0  0.2300  -0.4382  0.8982

 V1  V2  V3 
"a" "a" "a"

ブロック要因（REP）の多重比較

> AnovaModel.4 <- aov(DATA ~ REP + N + V + N:V, data=Dataset)

> summary(AnovaModel.4)

            Df Sum Sq Mean Sq F value  Pr(>F)    
REP          3   2.60   0.867   5.729 0.00222 ** 
N            4  41.23  10.309  68.153 < 2e-16 ***
V            2   1.05   0.526   3.480 0.03995 *  
N:V          8   2.29   0.286   1.893 0.08671 .  
Residuals   42   6.35   0.151                    
---
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 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 + 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.41712   
3 - 1 == 0   0.0140     0.1420   0.099  0.99965   
4 - 1 == 0  -0.4989     0.1420  -3.513  0.00579 **
3 - 2 == 0   0.2343     0.1420   1.650  0.36275   
4 - 2 == 0  -0.2787     0.1420  -1.962  0.21847   
4 - 3 == 0  -0.5129     0.1420  -3.612  0.00415 **
---
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 + N + V + N:V, data = Dataset)

Quantile = 2.6761
95% family-wise confidence level
 

Linear Hypotheses:
           Estimate lwr     upr    
2 - 1 == 0 -0.2203  -0.6003  0.1598
3 - 1 == 0  0.0140  -0.3660  0.3940
4 - 1 == 0 -0.4989  -0.8790 -0.1189
3 - 2 == 0  0.2343  -0.1458  0.6143
4 - 2 == 0 -0.2787  -0.6587  0.1014
4 - 3 == 0 -0.5129  -0.8930 -0.1329

   1    2    3    4 
 "b" "ab"  "b"  "a"

三中信宏「生物統計学」2018年6月7日実習記録

MINAKA Nobuhiro

AUTOMATIC

R Commander 起動

データファイル（Box3_R.txt）のインポート

データ表示

ブロック要因の因子化

箱ひげ図描画

窒素要因（N）

品種要因（V）

ブロック要因（REP）

線形統計モデル構築と当てはめ

分散分析表

処理平均のプロット

Tukey HSD 検定

窒素要因（N）の多重比較

品種要因（V）の多重比較

ブロック要因（REP）の多重比較