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

datasets"foster: HSAUR"による,ラット親子の遺伝子型と子供の体重の差の分析: 二元配置の分散分析と多重比較

(因子Aに対応なしX因子Bに対応なし・標本数が異なる) (多重比較・交互作用なし)


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

    くり返しがある二元配置のデータを分散分析することは,

    総変動から,因子A,因子Bによる変動と交互作用AXBによる変動と誤差変動を分離して,測定値への影響を検討することである.



    ラットの子供と母親の遺伝子型の組み合わせの違いにより子供の体重(g)の違いがあるかを検討した.

    子供の遺伝子型はA,B,I,Jとし,母親の遺伝子型はA,B,I,Jとした.

    標本数が異なる場合の分析になる.

    親子の遺伝子型の違いで体重の差があるか.2因子の交互作用はあるかを検討する.



    表1 ラット親子の遺伝子型と子供の体重(g)を検討する実験計画

    基本インストールではこのデータセットは準備されていないので,データセットのパッケージ「HSAUR」を新規にインストール.

    library("HSAUR") # でRに読み込む.
    data() # で確認.
    ? foster # でデータの解説が読める.
    litgen: genotype of the litter, a factor with levels A, B, I, and J.
    motgen: genotype of the mother, a factor with levels A, B, I, and J.
    weight: the weight of the litter after a feeding period.

    foster # でデータを表示する.61サンプル


    ラットの子供の遺伝子型(A,B,I,J)(因子A・4水準),母親の遺伝子型(A,B,I,J)(因子B・水準)の組み合わせ,4x4=16のケースを考える.

    それぞれの枠の実験に5匹のラットを使用する.4x4x5=80 サンプルで計画したが5例集まらなかった実験枠が出てきた.

     
    因子A litgen(子供の遺伝子型)
    因子B motgen(親の遺伝子型) 水準A1(A) 水準A2(B) 水準A3(I) 水準A4(J)
    水準B1(A) 実験01  5例 実験05  4例 実験09  3例 実験13  4例
    水準B2(B) 実験02  3例 実験06  5例 実験10  3例 実験14  3例
    水準B2(I) 実験03  4例 実験07  4例 実験11  5例 実験15  3例
    水準B2(J) 実験04  5例 実験08  2例 実験12  3例 実験16  5例


    表2 ラット親子の遺伝子型と子供の体重(g)を検討する実験計画 「対応なし・くり返し数が同じではない」

    因子A
    子供の遺伝子型
    水準A1(A) 水準A2(B) 水準A3(I) 水準A4(J)
    因子B
    親の遺伝子型
    水準B1(A) 水準B2(B) 水準B3(I) 水準B4(J) 水準B1(A) 水準B2(B) 水準B3(I) 水準B4(J) 水準B1(A) 水準B2(B) 水準B3(I) 水準B4(J) 水準B1(A) 水準B2(B) 水準B3(I) 水準B4(J)
    サンプル1 実験01-1 実験02-1 実験03-1 実験04-1 実験05-1 実験06-1 実験07-1 実験08-1 実験09-1 実験10-1 実験11-1 実験12-1 実験13-1 実験14-1 実験15-1 実験16-1
    サンプル2 実験01-2 実験02-2 実験03-2 実験04-2 実験05-2 実験06-2 実験07-2 実験08-2 実験09-2 実験10-2 実験11-2 実験12-2 実験13-2 実験14-2 実験15-2 実験16-2
    ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・ ・・
    サンプル5 実験01-5     実験04-5   実験06-5         実験11-5         実験16-5
    平均値                                



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


    表3 ラット親子の遺伝子型と子供の体重(g)を
    検討する実験計画 「スタック形式」


    データは「表シート」をマウスでコピーし入力.

    変数の「母親と子の遺伝子型」のデータ形式はfactor.

       litgen motgen weight
    1       A      A   61.5
    2       A      A   68.2
    3       A      A   64.0
    4       A      A   65.0
    5       A      A   59.7
    6       A      B   55.0
      ・・・・・・・
    56      J      I   61.4
    57      J      J   44.8
    58      J      J   51.5
    59      J      J   53.0
    60      J      J   42.0
    61      J      J   54.0
    
    (1) 「スタック形式」での操作の手順

    > x <- read.table ("clipboard", header=TRUE )
    > x 
       litgen motgen weight
    1       A      A   61.5
    2       A      A   68.2
    3       A      A   64.0
    4       A      A   65.0
    5       A      A   59.7
    6       A      B   55.0
      ・・・・・・・
    57      J      J   44.8
    58      J      J   51.5
    59      J      J   53.0
    60      J      J   42.0
    61      J      J   54.0
    
    > str (x) 
    'data.frame':   61 obs. of  3 variables:
     $ litgen: Factor w/ 4 levels "A","B","I","J": 1 1 1 1 1 1 1 1 1 1 ...
     $ motgen: Factor w/ 4 levels "A","B","I","J": 1 1 1 1 1 2 2 2 3 3 ...
     $ weight: num  61.5 68.2 64 65 59.7 55 42 60.2 52.5 61.8 ...
    


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

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

    表3 ラット親子の遺伝子型と子供の体重(g)を検討する実験計画「対応なし・くり返し数が同じではない」 「スタック形式」の緑の部分をマウスで選択し,これをコピーする.

    次に,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が理解した,サンプル・分析値・水準の3変数名のリストと変数ごとのカテゴリの数とカテゴリ名が示される.

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


    注意: 子供の遺伝子型と母親の遺伝子型の因子のデータ形式が [Factor] 形式でないときは,

    x$litgen <- factor( x$litgen )

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


  3. 1因子ごとの基礎統計


    a. 因子AB別の基礎統計量

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

    by ( x$weight, x$litgen, summary)

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


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

    plot.design(x) # 因子ごとの平均値のプロット図

    plotMeans( x$weight, x$litgen, error.bars="se")

    plotMeans( x$weight, x$motgen, error.bars="se")

    と記述し,因子(親と子の遺伝子型)ごとの体重の平均値の推移図を表示する.

    c. グラフの観察

    (1) litgen,x$weightとx$litgenのグラフでは3群の平均値に差を観察できない.子供の遺伝子型は体重の平均値の差に影響を与えていないようだ.

    (2) motgen,x$weightとx$motgenのグラフでは3群の平均値に差を観察できる.母親の遺伝子型は体重の平均値の差に影響を与えていると想像できる.

    これらの想像が可能かを分散分析で検定する.



  4. 2因子間の基礎統計 2因子間AXBの交互作用の可能性検討


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

    どの因子間の水準の組み合わせの時に,子供の体重が重いかをグラフで検討する.

    plotMeans( x$weight, x$litgen, x$motgen, error.bars="se")

    plotMeans( x$weight, x$motgen, x$litgen, error.bars="se") # 軸を変えてみる

    b. グラフの観察

    これらのグラフより,子供の遺伝子型の違いよりも親の遺伝子型の違いにより,子供の体重の差が現れていることが観察できる.

    2本の線が互いに交差するlitgen X motgenのグラフからは因子間の交互作用があると想像できる.

    これらの想像が可能かを分散分析で検定する.



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


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

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

    2つの因子,母親と子供の遺伝子型のそれぞれの水準の分散が等しいか,あるいは等しくないかを検定する.2水準のデータの分散は等質であるとの帰無仮説を立てる.

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

    bartlett.test ( x$weight ~ x$litgen )

    bartlett.test ( x$weight ~ x$motgen )

    と記述し,「スタック形式」のデータで,分散の等質性の検定を行う.

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

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


    > bartlett.test ( x$weight ~ x$litgen )
    
            Bartlett test of homogeneity of variances
    
    data:  x$weight by x$litgen 
    Bartlett's K-squared = 6.1503, df = 3, p-value = 0.1045
    
    > bartlett.test ( x$weight ~ x$motgen )
    
            Bartlett test of homogeneity of variances
    
    data:  x$weight by x$motgen 
    Bartlett's K-squared = 4.0342, df = 3, p-value = 0.2578
    

    分析結果

    2つの因子,母親と子供の遺伝子型は,ともに,p値>.05であるので帰無仮説が採択され,因子AとBの各水準の分散は等しいとする.



  6. 「因子Aに対応なしX因子Bに対応なし・標本数は異なる」 aov関数による二元配置分散分析


    子供の遺伝子型(因子A)で子供の体重に差があるか,また,母親の遺伝子型(因子B)で子供の体重に差があるか,

    そして,因子間の交互作用AXBがあるかを有意水準5%で検討する.

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

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



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

    summary ( aov ( weight ~ litgen * motgen, data=x ) )

    summary ( aov ( x$weight ~ x$litgen * x$motgen) ) # 同じ内容の分析をする.

    と記述し,「スタック形式」のデータで,「繰り返しがある二元配置の分散分析」を行う.

    > summary ( aov ( weight ~ litgen * motgen, data=x ) )
    
                  Df  Sum Sq Mean Sq F value   Pr(>F)   
    litgen         3   60.16   20.05  0.3697 0.775221   
    motgen         3  775.08  258.36  4.7632 0.005736 **
    litgen:motgen  9  824.07   91.56  1.6881 0.120053   
    Residuals     45 2440.82   54.24                    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    分散分析表が表示される.これを整形すると,以下のようになる.
    -----------------------------------------------------------------            (変動) (不偏分散) (分散比) 変動要因   自由度 平方和 平均平方和 F値 P値 ----------------------------------------------------------------- 子供の遺伝子型 3 60.16 20.05 0.3697 0.775 母親の遺伝子型 3 775.08 258.36 4.7632 0.006 ** 子供X母親 9 824.07 91.56 1.6881 0.120 誤差 45 2440.82 54.24 ----------------------------------------------------------------- ** p<.01 /font>
    分析結果

    有意水準5%で検定を行ったところ,

    (1) 子供の遺伝子型による差はP値>.05であるから,帰無仮説を採択する.子供の遺伝子型による差はないとする.

    (2) 母親の遺伝子型による差はP値<.05であるから帰無仮説を棄却する.母親の遺伝子型による差はあるとする.

    (3) 子供の遺伝子型X母親の遺伝子型の交互作用はP値>.05であるから帰無仮説を採択する.交互作用があるとはいえない.

    以上より,母親の遺伝子型の違いは,子供のラットの体重に影響を与えているといえる.

    母親の遺伝子型Bで体重が重くなり,遺伝子型Jで体重が軽くなること推察できる.次に,これを多重分析で確認する..



  7. 「交互作用なし・因子Aに対応なしX因子Bに対応なし・標本数は異なる」場合の多重比較


    二元配置では,交互作用がなければ,因子A・因子Bについて多重比較を検討すればよい.交互作用があれば,因子A・因子B の多重比較は不要になる.


    分散分析では,「母親の遺伝子型」のどの水準とどの水準に差があるかは,明らかにしてくれない.

    このため多重比較の事後検定(下位検定)を行う.

    テューキー(Tukey HSD)の方法は対応なしで標本数が同じ場合に用いることになっている.

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


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

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

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

      (1) 母親の遺伝子型と子供の体重

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

      (2) 子供の遺伝子型と子供の体重 分散分析で効果が無かったことがわかるので,この分析はしなくてもよい.

      tukey ( x$weight, x$litgen )

      > source("http://aoki2.si.gunma-u.ac.jp/R/src/tukey.R", encoding="euc-jp")
      > tukey ( x$weight, x$motgen )
      $result1
              n    Mean Variance
      Group1 16 55.4000 97.79600
      Group2 14 58.7000 52.44769
      Group3 16 53.3625 41.63850
      Group4 15 48.6800 39.65600
      
      $Tukey
                  t           p
      1:2 1.1800192 0.641809664
      1:3 0.7541441 0.874519865
      1:4 2.4468412 0.079747356
      2:3 1.9085916 0.235968142
      2:4 3.5284999 0.004508864
      3:4 1.7049604 0.330631015
      
      $phi
      [1] 57
      
      $v
      [1] 58.39511
      
      > tukey ( x$weight, x$litgen )
      $result1
              n     Mean  Variance
      Group1 17 55.11176  74.55235
      Group2 15 54.66667  50.88952
      Group3 14 52.90714 127.08841
      Group4 15 52.97333  34.46638
      
      $Tukey
                   t         p
      1:2 0.14924478 0.9988010
      1:3 0.72558804 0.8864639
      1:4 0.71703243 0.8899254
      2:3 0.56241148 0.9427535
      2:4 0.55083493 0.9459505
      3:4 0.02115702 0.9999966
      
      $phi
      [1] 57
      
      $v
      [1] 70.87666
      

      分析結果

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

      母親の遺伝子型の[$Tukey]の欄で有意水準5%以下のpを探すと,

      [2:4 p=0.004508864]のB-J間に5%水準で有意な差があることがわかる.遺伝子型の記号は1(A),2(B),3(I),4(J)であるから,

      母親の遺伝子型の水準ではB型とJ型の間に有意な差がある.B>J

      子供の遺伝子型の水準では有意な差がない.

      子供の体重は,母親の遺伝子型Bで重くなり,遺伝子型Jで軽くなることを示している.



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

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

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



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

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

      (1) 母親の遺伝子型と子供の体重

      TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "motgen", ordered = TRUE ) # テューキー(Tukey HSD)の方法

      plot(TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "motgen" )) # 母親の遺伝子型の多重比較のグラフ作成

      (2) 子供の遺伝子型と子供の体重 分散分析で効果が無かったことがわかるので,この分析はしなくてもよい.

      TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "litgen", ordered = TRUE ) # テューキー(Tukey HSD)の方法

      plot(TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "litgen" )) # 子供の遺伝子型の多重比較のグラフ作成

      と記述し,テューキー(Tukey HSD)の方法で多重比較を行う.

      > TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "motgen", ordered = TRUE )
        Tukey multiple comparisons of means
          95% family-wise confidence level
          factor levels have been ordered
      
      Fit: aov(formula = weight ~ litgen * motgen, data = x)
      
      $motgen母親の遺伝子型
              diff        lwr      upr     p adj
      I-J 4.670593 -2.3905240 11.73171 0.3035490
      A-J 6.566168 -0.4949498 13.62729 0.0767540
      B-J 9.896537  2.5954489 17.19762 0.0040509
      A-I 1.895574 -5.0507207  8.84187 0.8853702
      B-I 5.225943 -1.9641552 12.41604 0.2266493
      B-A 3.330369 -3.8597294 10.52047 0.6078581
      
      > plot(TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "motgen" ))
      
      > TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "litgen", ordered = TRUE )
        Tukey multiple comparisons of means
          95% family-wise confidence level
          factor levels have been ordered
      
      Fit: aov(formula = weight ~ litgen * motgen, data = x)
      
      $litgen子供の遺伝子型
                diff       lwr      upr     p adj
      J-I 0.06619048 -7.234897 7.367278 0.9999948
      B-I 1.75952381 -5.541564 9.060611 0.9174672
      A-I 2.20462185 -4.886102 9.295346 0.8402108
      B-J 1.69333333 -5.480769 8.867436 0.9219576
      A-J 2.13843137 -4.821470 9.098333 0.8448689
      A-B 0.44509804 -6.514804 7.405000 0.9982079
      
      > plot(TukeyHSD (aov (weight ~ litgen* motgen, data=x ), "litgen" ))
      

      分析結果

      2つの因子についての多重比較の結果をグラフで示す. このグラフの中で0を含まない水準間の差が有意であると読める.

      (1) 母親の遺伝子型の水準ではB型とJ型の間に[B-J p=0.0040509] に有意な差がある.B>J

      (2) 子供の遺伝子型の水準では体重に有意な差がない.

      子供の体重は,母親の遺伝子型Bで重くなり,遺伝子型Jで軽くなることを示している.



    3. ボンフェローニ(Bonferroni)法

      (1) 母親の遺伝子型と子供の体重

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

      pairwise.t.test(x$weight, x$motgen, p.adjust.method="bonferroni", paired=TRUE ) # これではなく

      pairwise.t.test(x$weight, x$motgen, p.adjust.method="bonferroni" )


      と記述し,Bonferroni法で多重比較を行う.

      > pairwise.t.test(x$weight, x$motgen, p.adjust.method="bonferroni" )
      
              Pairwise comparisons using t tests with pooled SD 
      
      data:  x$weight and x$motgen 
      
        A     B     I    
      B 1.000 -     -    
      I 1.000 0.368 -    
      J 0.105 0.005 0.562
      
      P value adjustment method: bonferroni

      分析結果

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

      母親の遺伝子型の水準B型とJ型の間に[B-J p=0.005] に有意な差がある.B>J

      子供の体重は,母親の遺伝子型Bで重くなり,遺伝子型Jで軽くなることを示している.



    4. Holm法

      (1) 母親の遺伝子型と子供の体重

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

      pairwise.t.test(x$weight, x$motgen, p.adjust.method="holm", paired=TRUE ) # これではなく

      pairwise.t.test(x$weight, x$motgen, p.adjust.method="holm")


      と記述し,Holm法で多重比較を行う.

      > pairwise.t.test(x$weight, x$motgen, p.adjust.method="holm")
      
              Pairwise comparisons using t tests with pooled SD 
      
      data:  x$weight and x$motgen 
      
        A     B     I    
      B 0.486 -     -    
      I 0.486 0.245 -    
      J 0.088 0.005 0.281
      
      P value adjustment method: holm

      分析結果

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

      母親の遺伝子型の水準B型とJ型の間に[B-J p=0.005] に有意な差がある.B>J

      子供の体重は,母親の遺伝子型Bで重くなり,遺伝子型Jで軽くなることを示している.

くり返しのある二元配置分散分析 (<対応なしX対応なし・標本数は異なる)     [ 目次に戻る ]