Last updated: October 2009. Kajiyama                    [ 目次に戻る ]

datasets"chickwts"による,六種類の飼料サプリメントと鶏の体重

くり返し数が同じでない一元配置分散分析と多重比較 (対応なし・標本数が異なる)


  1. 「対応なし・標本数が異なる」場合の課題

    三以上の群の平均値を比較するのであるが,各群の標本数が同じでない場合がある.

    例えば実験に失敗してデータを取れない場合や,予定した被験者が当日いない場合など.実験装置をゴジラが壊したこともありました.

    これらのデータを大切に扱うのが「各群の標本数が異なる場合の分散分析」である.



    鶏の餌にサプリメント混ぜ与え,体重を増やすのに効果的なサプリメントを探している.

    六種類のサプリメントを準備し,種類ごとに14羽の鶏を飼育した.しかし,途中で消えたものもあり,標本数が異なってしまった.

    六種類のサプリメント(水準A1〜水準A6)と鶏の体重(weight 測定値)の関係を表したのが表1である.

    サプリメントごとの鶏の体重の母集団の平均値に差があるといえるか,有意水準5%で検討する.


    表1 六種類の飼料サプリメントと鶏の体重 「対応なし・標本数が異なる アンスタック形式」

    このデータセット"chickwts"はすでにRに組み込まれている.
    > ? chickwts # データセット解説


    Description: An experiment was conducted to measure and compare the effectiveness of various feed supplements on the growth rate of chickens.
    A data frame with 71 observations on 2 variables.
    feed: a factor giving the feed type.
    weight: a numeric variable giving the chick weight.



    サンプル
    01
    02
    03
    ・・・
    13
    14
    平均値
    水準A1
    casein
    牛乳カス?
    水準A2
    horsebean
    ソラマメ
    水準A3
    linseed
    亜麻仁
    水準A4
    meatmeal
    水準A5
    soybean
    大豆
    水準A6
    sunflower
    ヒマワリの種
               
               
               
               
               
               
    323.6 160.2 218.8 276.9 246.4 328.9


  2. 分散分析の関数は「スタック形式」のデータで処理する.


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


    表2 六種類の飼料サプリメントと鶏の体重 「対応なし・標本数が異なる スタック形式」

       weight      feed
    1     179 horsebean
    2     160 horsebean
    3     136 horsebean
    4     227 horsebean
    5     217 horsebean
    6     168 horsebean
    7     108 horsebean
    8     124 horsebean
    9     143 horsebean
    10    140 horsebean
    11    309   linseed
    12    229   linseed
    13    181   linseed
    14    141   linseed
    15    260   linseed
    16    203   linseed
    17    148   linseed
    18    169   linseed
    19    213   linseed
    20    257   linseed
    21    244   linseed
    22    271   linseed
    23    243   soybean
    24    230   soybean
    25    248   soybean
    26    327   soybean
    27    329   soybean
    28    250   soybean
    29    193   soybean
    30    271   soybean
    31    316   soybean
    32    267   soybean
    33    199   soybean
    34    171   soybean
    35    158   soybean
    36    248   soybean
    37    423 sunflower
    38    340 sunflower
    39    392 sunflower
    40    339 sunflower
    41    341 sunflower
    42    226 sunflower
    43    320 sunflower
    44    295 sunflower
    45    334 sunflower
    46    322 sunflower
    47    297 sunflower
    48    318 sunflower
    49    325  meatmeal
    50    257  meatmeal
    51    303  meatmeal
    52    315  meatmeal
    53    380  meatmeal
    54    153  meatmeal
    55    263  meatmeal
    56    242  meatmeal
    57    206  meatmeal
    58    344  meatmeal
    59    258  meatmeal
    60    368    casein
    61    390    casein
    62    379    casein
    63    260    casein
    64    404    casein
    65    318    casein
    66    352    casein
    67    359    casein
    68    216    casein
    69    222    casein
    70    283    casein
    71    332    casein
    


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


    > y <- read.table ("clipboard", header=TRUE )
    > y
       weight      feed
    1     179 horsebean
    2     160 horsebean
    3     136 horsebean
    4     227 horsebean
      ・・・・・
    68    216    casein
    69    222    casein
    70    283    casein
    71    332    casein
    
    > str (y)
    'data.frame':   71 obs. of  2 variables:
     $ weight: num  179 160 136 227 217 168 108 124 143 140 ...
     $ feed  : Factor w/ 6 levels "casein","horsebean",..: 2 2 2 2 2 2 2 2 2 2 ...
    


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

    ここでは,表シートあるいはエクセルの表計算シートからクリップボード経由で読み込む操作を行う.

    表2 六種類の飼料サプリメントと鶏の体重 「スタック形式」の緑の部分をマウスで選択し,これをコピーする.

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

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

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

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

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


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

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


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

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

    Rが理解した,鶏の体重(weight 測定値)と飼料サプリメントの2変数名のリストと変数ごとのカテゴリの数とカテゴリ名が示される.

    ・分析で使用する変数名とカテゴリ名はRが理解した名前を使用する.


    注意: 水準A1〜A6を含む因子(飼料サプリメント)のデータ形式が [Factor] 形式でないときは,正しい統計処理ができなくなる.

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


  4. 1因子Aの基礎統計


    a. 因子Aの基礎統計量

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

    by ( y$weight, y$feed, summary ) #   と記述し,「スタック形式」のデータで,群別の基礎統計利用を計算する.

    > by ( y$weight, y$feed, summary )
    y$feed: casein
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      216.0   277.2   342.0   323.6   370.8   404.0 
    --------------------------------------------------------------- 
    y$feed: horsebean
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      108.0   137.0   151.5   160.2   176.2   227.0 
    --------------------------------------------------------------- 
    y$feed: linseed
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      141.0   178.0   221.0   218.8   257.8   309.0 
    --------------------------------------------------------------- 
    y$feed: meatmeal
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      153.0   249.5   263.0   276.9   320.0   380.0 
    --------------------------------------------------------------- 
    y$feed: soybean
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      158.0   206.8   248.0   246.4   270.0   329.0 
    --------------------------------------------------------------- 
    y$feed: sunflower
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      226.0   312.8   328.0   328.9   340.2   423.0 
    


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

    boxplot ( y$weight ~ y$feed) #  と記述し,箱ひげグラフを表示する.

    plotMeans ( y$weight, y$feed, error.bars="se") #  と記述し,平均値の推移図を表示する.

    c. グラフの観察




  5.  バートレット(Bartlett)の検定--分散の等質性(等分散性)の検定


    飼料サプリメント6水準(6群)の分散が等しいか,あるいは等しくないかを検定する.

    6水準のデータの分散は等質であるとの帰無仮説を立てる.

    分散分析( R関数 aov )は,分散が等しいという仮定の上で構築されているので標本の測定値にもとづいて分散が一様であるか検討する.

    ただし,分散分析のF検定はその前提条件が崩れても,検定結果は信頼できるといわれている.


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

    bartlett.test ( y$weight ~ y$feed ) # と記述し,「スタック形式」のデータで,分散の等質性の検定を行う.

    バートレットの検定は水準ごとの標本の数が同じでない場合も使用できる.

    p値が0.05以下のときに水準の分散は等しくないと結論する.


    > bartlett.test ( y$weight ~ y$feed )
    
            Bartlett test of homogeneity of variances
    
    data:  y$weight by y$feed 
    Bartlett's K-squared = 3.2597, df = 5, p-value = 0.66
    

    分析結果

    p値が0.66であり,p値>0.05であるので各水準の分散は等しいと結論し,分散分析に進む.



  6. 「対応なし・標本数が異なる」場合の一元配置・分散分析


    「6種類の飼料サプリメント(A1〜A6)で飼育した鶏の体重(weight 測定値)の母集団の平均値に差があるといえるか」,有意水準5%で検討する.

    帰無仮説は「各群の母集団の鶏の体重の平均値は等しい」.

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



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

    summary ( aov ( y$weight ~ y$feed )) # と記述し,「スタック形式」のデータで,分散分析を行う.

    summary ( aov ( weight ~ feed, y )) # の書式も同じ結果となる.


    > summary ( aov ( y$weight ~ y$feed ))
    
                Df Sum Sq Mean Sq F value    Pr(>F)    
    y$feed       5 231129   46226  15.365 5.936e-10 ***
    Residuals   65 195556    3009                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    有意性の記号 :    ***  (0.1%で有意)  **  (1%で有意)   *  (5%で有意)   . (有意の傾向がある)
    
    
    分散分析表が表示される.これを整形すると,以下のようになる.
    -----------------------------------------------------------          (変動)  (不偏分散) (分散比) 変動要因 自由度 平方和  平均平方和 F値 P値 ----------------------------------------------------------- 因子A 5 231129 46226 15.365 5.936e-10 *** 誤差E 65 195556 3009 -----------------------------------------------------------
    分析結果

    因子(feed)AのA1〜A6の水準間の変動について,P値>0.05であることから,5%水準で帰無仮説は棄却される.

    因子Aは有意差がある.飼料サプリメントの種類間で体重の平均値に違いがある.



  7. 「対応なし・標本数が異なる」場合の多重比較


    有意差があれば,どの水準の飼料サプリメントの種類が効率よく体重を増やすことができるかを知りたくなる.

    しかし分散分析は,どの水準とどの水準に差があるかは,明らかにできない.

    このため多重比較という事後検定(下位検定)に進み,水準(群)間の比較を行う.

    Tukey(テューキー)の方法は各水準の標本数が数が同じ場合に用いる.

    標本数が異なる場合はテューキー・クレーマーの方法 Tukey-Kramer's method を用いる.



    1. 「対応なし・標本数が異なる」場合の青木の Tukey-Kramer(テューキー・クレーマー)の方法

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

      source("http://aoki2.si.gunma-u.ac.jp/R/src/tukey.R", encoding="euc-jp") #  を貼り付け青木先生のプログラムとリンクさせる.

      # 心配性の受講生への青木先生のコメント

      tukey ( y$weight, y$feed ) #   と記述し,青木のTukey-Kramer(テューキー・クレーマー)の方法で多重比較を行う.

      > source("http://aoki2.si.gunma-u.ac.jp/R/src/tukey.R", encoding="euc-jp")
      > tukey ( y$weight, y$feed )
      $result1
              n     Mean Variance
      Group1 12 323.5833 4151.720
      Group2 10 160.2000 1491.956
      Group3 12 218.7500 2728.568
      Group4 11 276.9091 4212.091
      Group5 14 246.4286 2929.956
      Group6 12 328.9167 2384.992
      
      $Tukey
                  t            p
      1:2 6.9567776 3.070197e-08
      1:3 4.6816194 2.100151e-04
      1:4 2.0385502 3.324584e-01
      1:5 3.5756235 8.365309e-03
      1:6 0.2381746 9.998902e-01
      2:3 2.4930286 1.413329e-01
      2:4 4.8698150 1.062092e-04
      2:5 3.7969132 4.216654e-03
      2:6 7.1838681 1.219886e-08
      3:4 2.5401639 1.276965e-01
      3:5 1.2827225 7.932853e-01
      3:6 4.9197940 8.843233e-05
      4:5 1.3792208 7.391356e-01
      4:6 2.2714895 2.206962e-01
      5:6 3.8227890 3.884521e-03
      
      $phi
      [1] 65
      
      $v
      [1] 3008.554
      


      分析結果

      多重比較の帰無仮説は二群間の平均値には差が無いである.

      [$Tukey]の欄で有意水準5%以下のp値を探すと,[1:2][1:3][1:5][2:4][2:5][2:6][3:6][5:6]の8個の水準間に有意差がある.

      caseinとsunflowerが体重増加に効果があった.



    2. 「対応なし・標本数が異なる」場合の テューキーの方法(Tukey's honestly significant difference test)

      標本数が異なる場合はテューキー・クレーマーの方法 Tukey-Kramer's method を用いる.

      しかし,Rの関数「TukeyHSD」(テューキーの方法) は標本数が異なった場合も使えるようにTukey-Kramerのへ拡張が行われているかもしれない?.

      > help(TukeyHSD)

      "Technically the intervals constructed in this way would only apply to balanced designs where there are the same number of observations made at each level of the factor.

      This function incorporates an adjustment for sample size that produces sensible intervals for mildly unbalanced designs."



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

      TukeyHSD ( aov ( y$weight ~ y$feed) ) #   と記述し,Tukey(テューキー)の方法で多重比較を行う.

      require(graphics) # グラフパッケージの読み込み

      plot (TukeyHSD ( aov ( y$weight ~ y$feed) ) ) # 結果をグラフで表す


      > TukeyHSD ( aov ( y$weight ~ y$feed) ) 
        Tukey multiple comparisons of means
          95% family-wise confidence level
      
      Fit: aov(formula = y$weight ~ y$feed)
      
      $`y$feed`
                                 diff         lwr       upr     p adj
      horsebean-casein    -163.383333 -232.346876 -94.41979 0.0000000
      linseed-casein      -104.833333 -170.587491 -39.07918 0.0002100
      meatmeal-casein      -46.674242 -113.906207  20.55772 0.3324584
      soybean-casein       -77.154762 -140.517054 -13.79247 0.0083653
      sunflower-casein       5.333333  -60.420825  71.08749 0.9998902
      linseed-horsebean     58.550000  -10.413543 127.51354 0.1413329
      meatmeal-horsebean   116.709091   46.335105 187.08308 0.0001062
      soybean-horsebean     86.228571   19.541684 152.91546 0.0042167
      sunflower-horsebean  168.716667   99.753124 237.68021 0.0000000
      meatmeal-linseed      58.159091   -9.072873 125.39106 0.1276965
      soybean-linseed       27.678571  -35.683721  91.04086 0.7932853
      sunflower-linseed    110.166667   44.412509 175.92082 0.0000884
      soybean-meatmeal     -30.480519  -95.375109  34.41407 0.7391356
      sunflower-meatmeal    52.007576  -15.224388 119.23954 0.2206962
      sunflower-soybean     82.488095   19.125803 145.85039 0.0038845
      



      分析結果

      多重比較の帰無仮説は二群間の平均値には差が無いである.

      グラフの見方は横棒グラフが,ゼロの値の位置をまたがなければ,有意であると読む.

      [p adj]の欄で有意水準5%以下を探すと,8個の水準間に有意差があった.caseinとsunflowerが体重増加に効果があった.


    3. 「対応なし・標本数が異なる」場合の ボンフェローニ(Bonferroni)法

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

      pairwise.t.test( y$weight, y$feed, p.adjust.method="bonferroni" ) # と記述し,Bonferroni法で多重比較を行う.

      >  pairwise.t.test( y$weight, y$feed, p.adjust.method="bonferroni" )
      
              Pairwise comparisons using t tests with pooled SD 
      
      data:  y$weight and y$feed 
      
                casein  horsebean linseed meatmeal soybean
      horsebean 3.1e-08 -         -       -        -      
      linseed   0.00022 0.22833   -       -        -      
      meatmeal  0.68350 0.00011   0.20218 -        -      
      soybean   0.00998 0.00487   1.00000 1.00000  -      
      sunflower 1.00000 1.2e-08   9.3e-05 0.39653  0.00447
      
      P value adjustment method: bonferroni 
      

      分析結果

      多重比較の帰無仮説は二群間の平均値には差が無いである.有意水準5%以下を捜すと,.

      8個の水準間に有意な差がある.


    4. 「対応なし・標本数が異なる」場合の Holm法

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

      pairwise.t.test( y$weight, y$feed, p.adjust.method="holm" ) # と記述し,Holm法で多重比較を行う.

      >  pairwise.t.test( y$weight, y$feed, p.adjust.method="holm" )
      
              Pairwise comparisons using t tests with pooled SD 
      
      data:  y$weight and y$feed 
      
                casein  horsebean linseed meatmeal soybean
      horsebean 2.9e-08 -         -       -        -      
      linseed   0.00016 0.09435   -       -        -      
      meatmeal  0.18227 9.0e-05   0.09435 -        -      
      soybean   0.00532 0.00298   0.51766 0.51766  -      
      sunflower 0.81249 1.2e-08   8.1e-05 0.13218  0.00298
      
      P value adjustment method: holm 
      

      分析結果

      多重比較の帰無仮説は二群間の平均値には差が無いである.有意水準5%以下を捜すと,.

      8個の水準間に有意な差がある.


    5. まとめ

      体重を増やすのに効果的な飼料サプリメントはcasein(牛乳カス?)とsunflower(ヒマワリの種)であり,逆にhorsebean(そら豆)の効果は低かった.


一元配置分散分析と多重比較 (対応なし・標本数が異なる)     [ 目次に戻る ]