Last updated: December 20120 -January 2010. Kajiyama                    [ 目次に戻る ]

datasets"iris"による,あやめのがく・花びらの形状分析

多変量分散分析 (多変量・対応あり・標本数は同じ)


  1. あやめ 「多変量・対応あり・標本数は同じ」場合の課題

    表1に示すように,測定値(特性値)が2種類以上(多変量)の実験・調査データがある場合,

    分散分析では測定値ごとに分割して一元配置として検討することになる.

    これに対し,測定値が多変量のまま各水準の平均値とその変動を分析するのが「多変量分散分析 MANOVA」である.



    あやめの花の形状が品種(因子A 3水準)により異なるか検討する.

    花は一箇所でなく,4部位の形状(測定値),がく片(Sepal)の長さと幅,および花びら(Petal)の長さと幅を測定した.(対応あり)

    表1の測定値から3品種のあやめの花の形状に差があるか,4種類の花の形状部位に差があるか.有意水準5%で検討する.


    表1 あやめのがく・花びらの形状「アンスタック形式」

     
    Number of cases: 150
    Variable Names:
    Sepal.Length: がく片の長さ
    Sepal.Width : がく片の幅
    Petal.Length: 花びらの長さ
    Petal.Width : 花びらの幅
    Species: あやめ3品種 
    [setosa・versicolor・virginica] 
      
     
    測定値 4部位
    因子A 品種
    Sepal.Length Sepal.Width Petal.Length Petal.Width
    A1
    setosa
    50例
    A2
    versicolor
    50例
    A2
    virginica
    50例


    表2 あやめのがく・花びらの形状 「アンスタック形式」

    因子 A 品種
    水準 A1 setosa 水準 A2 versicolor 水準 A3virginica
    測定値
    4部位
    Sepal.
    Length
    Sepal.
    Width
    Petal.
    Length
    Petal.
    Width
    Sepal.
    Length
    Sepal.
    Width
    Petal.
    Length
    Petal.
    Width
    Sepal.
    Length
    Sepal.
    Width
    Petal.
    Length
    Petal.
    Width
    サンプル01                        
    サンプル02                        
    ・・・                        
    サンプル50                        
    平均値 5.006 3.428 1.462 0.246 5.936 2.770 4.26 1.326 6.588 2.974 5.552 2.026



  2. あやめ Rに「スタック形式」のデータを読み込み,データの確認を行う


    表3 あやめのがく・花びらの形状「スタック形式」 データは「表シート」をマウスでコピーし入力する.

        Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
    1            5.1         3.5          1.4         0.2     setosa
    2            4.9         3.0          1.4         0.2     setosa
    3            4.7         3.2          1.3         0.2     setosa
    4            4.6         3.1          1.5         0.2     setosa
    5            5.0         3.6          1.4         0.2     setosa
        ・・・・・・・・
        ・・・・・・・・
    145          6.7         3.3          5.7         2.5  virginica
    146          6.7         3.0          5.2         2.3  virginica
    147          6.3         2.5          5.0         1.9  virginica
    148          6.5         3.0          5.2         2.0  virginica
    149          6.2         3.4          5.4         2.3  virginica
    150          5.9         3.0          5.1         1.8  virginica
    


    (1) 「スタック形式」での操作の手順

    > x <- read.table ("clipboard", header=TRUE )
    > x 
    1            5.1         3.5          1.4         0.2     setosa
    2            4.9         3.0          1.4         0.2     setosa
    3            4.7         3.2          1.3         0.2     setosa
    4            4.6         3.1          1.5         0.2     setosa
    5            5.0         3.6          1.4         0.2     setosa
        ・・・・・・・・
        ・・・・・・・・
    145          6.7         3.3          5.7         2.5  virginica
    146          6.7         3.0          5.2         2.3  virginica
    147          6.3         2.5          5.0         1.9  virginica
    148          6.5         3.0          5.2         2.0  virginica
    149          6.2         3.4          5.4         2.3  virginica
    150          5.9         3.0          5.1         1.8  virginica
    
    > str (x) 
    'data.frame':   150 obs. of  5 variables:
     $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
     $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
     $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
     $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
     $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
    


    (2) データフレームの読み込み

    ここでは,表シートからクリップボード経由で読み込む操作を行う.

    表シートの 表3 あやめのがく・花びらの形状「スタック形式」の緑の部分をマウスで選択し,これをコピーする.

    次に,Rの「コンソール」画面に,

    x <- read.table ("clipboard", header=TRUE

    と記述し,コピーした[clipboard]データファイルを,Rの内部のデータフレーム,ファイル名x に直接読み込む.

    header=TRUE は第1行が列の変数名になっていることを指示している.

    ) #    データをコピーした後で,カッコを入力する.ENTERーキーを押す.


    (3) 読み込んだデータフレームの確認

    x      と記述し,Rの内部に作成したデータフレーム x を,表示し確認する.


    (4) データフレームの構造確認

    str (x)   と記述し,Rの内部のデータフレーム x の3変数の内容を表示確認する.

    Rが理解した,因子名(変数名)のリストと因子ごとのカテゴリの数とカテゴリ名が示される.

    ・分析で使用する因子名と水準(カテゴリ)名はRが理解した名前を使用する.


    注意: Speciesのデータ形式が数値で [Factor] 形式でないときは,

    x$Species <- factor( x$Species ) と記述し,データフレームの水準の変数を[Facter]形式に変更する.



  3. 1因子Aの基礎統計


    a. 因子Aの基礎統計量  あやめの花の品種別の部位ごとの基礎統計量

    Rの「コンソール」画面に,

    by ( x$Sepal.Length, x$Species, summary)
    by ( x$Sepal.Width, x$Species, summary)
    by ( x$Petal.Length, x$Species, summary)
    by ( x$Petal.Width, x$Species, summary)
       と記述し,「スタック形式」のデータで,品種別(因子A水準別)の基礎統計利用を計算する.

    b. 因子Aの平均値の推移図   # plotMeans関数を使用するにはパッケージ 「Rcmdr」を読み込んでおくこと

    因子A 品種Speciesごとに4つの部位の測定値の「平均値の推移図」を作成する.

    plotMeans( x$Sepal.Length, x$Species, error.bars="se")

    plotMeans( x$Sepal.Width, x$Species, error.bars="se")

    plotMeans( x$Petal.Length, x$Species, error.bars="se")

    plotMeans( x$Petal.Width, x$Species, error.bars="se")



    > by ( x$Sepal.Length, x$Species, summary)
    x$Species: setosa
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      4.300   4.800   5.000   5.006   5.200   5.800 
    ---------------------------------------------------------------- 
    x$Species: versicolor
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      4.900   5.600   5.900   5.936   6.300   7.000 
    ---------------------------------------------------------------- 
    x$Species: virginica
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      4.900   6.225   6.500   6.588   6.900   7.900 
    
    >
    > by ( x$Sepal.Width, x$Species, summary)
    x$Species: setosa
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      2.300   3.200   3.400   3.428   3.675   4.400 
    ---------------------------------------------------------------- 
    x$Species: versicolor
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      2.000   2.525   2.800   2.770   3.000   3.400 
    ---------------------------------------------------------------- 
    x$Species: virginica
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      2.200   2.800   3.000   2.974   3.175   3.800 
    
    >
    > by ( x$Petal.Length, x$Species, summary)
    x$Species: setosa
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1.000   1.400   1.500   1.462   1.575   1.900 
    ---------------------------------------------------------------- 
    x$Species: versicolor
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
       3.00    4.00    4.35    4.26    4.60    5.10 
    ---------------------------------------------------------------- 
    x$Species: virginica
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      4.500   5.100   5.550   5.552   5.875   6.900
    
    >
    > by ( x$Petal.Width, x$Species, summary) 
    x$Species: setosa
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0.100   0.200   0.200   0.246   0.300   0.600 
    ---------------------------------------------------------------- 
    x$Species: versicolor
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1.000   1.200   1.300   1.326   1.500   1.800 
    ---------------------------------------------------------------- 
    x$Species: virginica
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1.400   1.800   2.000   2.026   2.300   2.500 
    
    > plotMeans( x$Sepal.Length, x$Species, error.bars="se") 
    > 
    > plotMeans( x$Sepal.Width, x$Species, error.bars="se") 
    > 
    > plotMeans( x$Petal.Length, x$Species, error.bars="se") 
    > 
    > plotMeans( x$Petal.Width, x$Species, error.bars="se") 
    
    
    
    
    
    
    


    c. グラフの観察

    品種Specieにより4つの部位の寸法の平均値が異なっていることが観察できる.

    Sepal.Width部位の平均値のパタ−ンが異なる.



  4. あやめ manova 関数による多変量分散分析


    (1) あやめの4つの部位の測定値をまとめて検討する.

    帰無仮説は「各水準の母集団の平均値は等しい」.

    対立仮説は,「全ての水準の組み合わせに平均値の差がある」のでなく,「各水準の少なくとも一組に平均値の差がある」.


    Rの「コンソール」画面に,

    result <- manova(cbind( x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species) #   と記述し,多変量分散分析を行う.

    多変量検定として、次の4 つを行う.

     ・Wilks のラムダ(Wilks’ Lambda)
     ・Pillai のトレース(Pillai’s Trace)
     ・Hotelling-Lawley のトレース(Hotelling-Lawley Trace)
     ・Roy の最大根(Roy’s Greatest Root)

    summary(result, test="Pillai") #   と記述し,4測定値まとめて多変量分散分析の結果で「Pillai のトレース検定」する.

    # 何も指定しないとPillai の検定になる.

    summary(result, test="Wilks") #   と記述し,4測定値まとめて多変量分散分析の結果で「Wilk's Lambdaの検定」する.

    summary(result, test="Hotelling-Lawley") #   と記述し,4測定値まとめて多変量分散分析の結果で「Hotelling-Lawleyのトレース検定」する.

    summary(result, test="Roy") #   と記述し,4測定値まとめて多変量分散分析の結果で「Roy の最大根検定」する.


    > result <- manova(cbind( x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species) 
    > result
    Call:
       manova(cbind(x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species)
    
    Terms:
                    x$Species Residuals
    resp 1            63.2121   38.9562
    resp 2            11.3449   16.9620
    resp 3           437.1028   27.2226
    resp 4            80.4133    6.1566
    Deg. of Freedom         2       147
    
    Residual standard error: 0.5147894 0.3396877 0.4303345 0.2046500 
    Estimated effects may be unbalanced
    
    > summary(result, test="Pillai") 
    
               Df Pillai approx F num Df den Df    Pr(>F)    
    x$Species   2  1.192   53.466      8    290 < 2.2e-16 ***
    Residuals 147                                            
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    > summary(result, test="Wilks")  
               Df   Wilks approx F num Df den Df    Pr(>F)    
    x$Species   2   0.023  199.145      8    288 < 2.2e-16 ***
    Residuals 147                                             
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    > summary(result, test="Hotelling-Lawley")  
               Df Hotelling-Lawley approx F num Df den Df    Pr(>F)    
    x$Species   2            32.48   580.53      8    286 < 2.2e-16 ***
    Residuals 147                                                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    > summary(result, test="Roy")  
               Df     Roy approx F num Df den Df    Pr(>F)    
    x$Species   2   32.19  1166.96      4    145 < 2.2e-16 ***
    Residuals 147                                             
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    

    分析結果

    有意水準5%でF検定を行ったところ,(4つの検定法で)P値<.05であるから,帰無仮説は棄却する.

    品種によってあやめの花の4つの部位からなる形状の平均値は異なるといえる.


    (2) あやめの品種間の平均値に差があるかを,花の部位ごとに検討する.これは部位の測定値ごとの一元配置の分散分析である.

    多変量分散分析の結果を元に,部位ごとに測定値に差があるかを検討する.

    summary.aov(result)

    > summary.aov(result) 
    
     Response 1 :測定値 Sepal.Length
                 Df Sum Sq Mean Sq F value    Pr(>F)    
    x$Species     2 63.212  31.606  119.26 < 2.2e-16 ***
    Residuals   147 38.956   0.265                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
     Response 2 :測定値 Sepal.Width 
                 Df  Sum Sq Mean Sq F value    Pr(>F)    
    x$Species     2 11.3449  5.6725   49.16 < 2.2e-16 ***
    Residuals   147 16.9620  0.1154                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
     Response 3 :測定値 Petal.Length
                 Df Sum Sq Mean Sq F value    Pr(>F)    
    x$Species     2 437.10  218.55  1180.2 < 2.2e-16 ***
    Residuals   147  27.22    0.19                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
     Response 4 :測定値 Petal.Width 
                 Df Sum Sq Mean Sq F value    Pr(>F)    
    x$Species     2 80.413  40.207     960 < 2.2e-16 ***
    Residuals   147  6.157   0.042                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    

    分析結果

    品種と測定値Sepal.Lengthは,P値<.05であるから,帰無仮説は棄却する.品種間でSepal.Lengthに違いがある

    品種と測定値Sepal.Width は,P値<.05であるから,帰無仮説は棄却する.品種間でSepal.Widthに違いがある

    品種と測定値Petal.Lengthは,P値<.05であるから,帰無仮説は棄却する.品種間でPetal.Lengthに違いがある

    品種と測定値Petal.Widは,P値<.05であるから,帰無仮説は棄却する.品種間でPetal.Widに違いがある



  5. あやめの2品種間ごとに差があるか比較検討する


    # 品種"setosa","versicolor"間の平均値の差の検定
    summary(manova(cbind(x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species,
    subset=x$Species %in% c("setosa","versicolor")))


    # 品種"setosa","virginica"間の平均値の差の検定
    summary(manova(cbind(x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species,
    subset=x$Species %in% c("setosa","virginica")))


    # 品種"versicolor","virginica"間の平均値の差の検定
    summary(manova(cbind(x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species,
    subset=x$Species %in% c("versicolor","virginica")))



    > summary(manova(cbind(x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species, 
    subset=x$Species %in% c("setosa","versicolor")))
              Df Pillai approx F num Df den Df    Pr(>F)    
    x$Species  1   0.96   625.46      4     95 < 2.2e-16 ***
    Residuals 98                                            
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    > 
    > summary(manova(cbind(x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species, 
    subset=x$Species %in% c("setosa","virginica")))
              Df  Pillai approx F num Df den Df    Pr(>F)    
    x$Species  1    0.98  1182.57      4     95 < 2.2e-16 ***
    Residuals 98                                             
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    > 
    > summary(manova(cbind(x$Sepal.Length, x$Sepal.Width, x$Petal.Length, x$Petal.Width) ~ x$Species, 
    subset=x$Species %in% c("versicolor","virginica")))
              Df Pillai approx F num Df den Df    Pr(>F)    
    x$Species  1  0.784   86.148      4     95 < 2.2e-16 ***
    Residuals 98                                            
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    

    分析結果

    帰無仮説は2つの品種間には差が無いである.有意水準5%で検定を行う.

    ["setosaとversicolor]には違いがある.P値<.05であるから,帰無仮説は棄却する.

    ["setosaとvirginica]には違いがある.P値<.05であるから,帰無仮説は棄却する.

    ["versicolorとvirginica]には違いがある.P値<.05であるから,帰無仮説は棄却する.

    まとめると,あやめの3つの品種間にそれぞれ花の形状の違いがある.

多変量分散分析 (多変量・対応あり・標本数は同じ)     [ 目次に戻る ]