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

datasets"Egyptian: HSAUR"による,古代エジプト人男性の頭蓋骨の年代分析

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


  1. 頭蓋骨「くり返しがあり標本数は同じ」場合の課題

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

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

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



    古代エジプト人男性の年代別の頭蓋骨の形状を分析して,紀元前4000年以降,同じ人種であったのかそれとも交流があり変化したのか.

    5種類の年代分け(因子A)と各年代30例づつ頭蓋骨の4種類の測定部位を測定した.

    表1の測定値から頭蓋骨の形状の年代変化を有意水準5%で検討する.


    表1 古代のエジプト人男性150名の頭蓋骨の形状

    データの出典:Egyptian Skulls  Story Names: Egyptian Skull Development
    Description: Four measurements of male Egyptian skulls from 5 different time periods. Thirty skulls are measured from each time period.
    MB: 頭蓋骨の幅 Maximal Breadth of Skull
    BH: 頭蓋骨ブレグマの高さ Basibregmatic Height of Skull
    BL: 頭蓋骨のバシオンプロスチオン長さ? Basialveolar Length of Skull
    NH: 頭蓋骨の鼻高 Nasal Height of Skull
    Year: Approximate Year of Skull Formation (negative = B.C., positive = A.D.)

    上記を入手できない場合,データセットのパッケージ「HSAUR」を新規にインストール.
    ただし,ここで引用した年代のデータ形式が異なっている.

    library("HSAUR") # でRに読み込む.
    data() # で確認.
    ? Egyptian # でデータの解説が読める
    Egyptian # でデータを読み込む.

     
    頭蓋骨の測定部位
    因子A
    年代
    MB BH BL NH
    水準A1 -4000年 30例
    水準A2 -3300年 30例
    水準A3 -1850年 30例
    水準A4 -200年 30例
    水準A5   150年 30例


    表2 古代のエジプト人男性150名の頭蓋骨の形状

    因子A 年代
    水準 A1
    -4000
    水準 A2
    -3300
    水準 A3
    -1850
    水準 A4
    -200
    水準 A5
    150
    測定部位
    MB BH BL NH MB BH BL NH MB BH BL NH MB BH BL NH MB BH BL NH
    サンプル01                                        
    サンプル02                                        
    ・・・                                        
    サンプル30                                        
    平均値 131.4 133.6 99.17 50.53 132.4 132.7 99.07 50.23 134.5 133.8 96.03 50.57 135.5 132.3 94.53 51.97 136.2 130.3 93.5 51.37



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


    表3 古代のエジプト人男性150名の頭蓋骨の形状
    「スタック形式」

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

    MB	BH	BL	NH	Year
    131	138	89	49	-4000
    125	131	92	48	-4000
    131	132	99	50	-4000
    119	132	96	44	-4000
    136	143	100	54	-4000
    138	137	89	56	-4000
    139	130	108	48	-4000
    125	136	93	48	-4000
        ・・・・・・・・
        ・・・・・・・・
    143	126	88	54	150
    134	124	91	55	150
    132	127	97	52	150
    137	125	85	57	150
    129	128	81	52	150
    140	135	103	48	150
    147	129	87	48	150
    136	133	97	51	150
    
       
    (1) 「スタック形式」での操作の手順



    > x <- read.table ("clipboard", header=TRUE )
    > x 
         MB  BH  BL NH  Year
    1   131 138  89 49 -4000
    2   125 131  92 48 -4000
    3   131 132  99 50 -4000
    4   119 132  96 44 -4000
    5   136 143 100 54 -4000
    6   138 137  89 56 -4000
    7   139 130 108 48 -4000
        ・・・・・・
    144 134 124  91 55   150
    145 132 127  97 52   150
    146 137 125  85 57   150
    147 129 128  81 52   150
    148 140 135 103 48   150
    149 147 129  87 48   150
    150 136 133  97 51   150
    
    > str (x) 
    'data.frame':   150 obs. of  5 variables:
     $ MB  : int  131 125 131 119 136 138 139 125 131 134 ...
     $ BH  : int  138 131 132 132 143 137 130 136 134 134 ...
     $ BL  : int  89 92 99 96 100 89 108 93 102 99 ...
     $ NH  : int  49 48 50 44 54 56 48 48 51 51 ...
     $ Year: int  -4000 -4000 -4000 -4000 -4000 -4000 - 4000 ...
    
     # 測定部位ごとの平均値の推移図
    >  plotMeans( x$MB, x$Year, error.bars="se")
    >  plotMeans( x$BH, x$Year, error.bars="se") 
    >  plotMeans( x$BL, x$Year, error.bars="se") 
    >  plotMeans( x$NH, x$Year, error.bars="se")
    


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

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

    表シートの 表3 古代のエジプト人男性150名の頭蓋骨の形状「スタック形式」の緑の部分をマウスで選択し,これをコピーする.

    次に,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が理解した名前を使用する.


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

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



  3. 1因子Aの基礎統計


    a. 因子Aの基礎統計量 年代別の頭蓋骨の4部位の測定値の基礎統計量

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

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


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

    年代ごとに頭蓋骨の部位の測定値の「平均値の推移図」を作成する.

    plotMeans( x$MB, x$Year, error.bars="se")

    plotMeans( x$BH, x$Year, error.bars="se")

    plotMeans( x$BL, x$Year, error.bars="se")

    plotMeans( x$NH, x$Year, error.bars="se")
    #   と記述し,年代因子と部位寸法の関係を平均値の推移図で表示する.

    > x$Year <- factor( x$Year )
    
    > by ( x$MB, x$Year, summary)
    x$Year: -4000
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      119.0   128.0   131.0   131.4   134.8   141.0 
    ----------------------------------------------------
    x$Year: -3300
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      123.0   130.0   132.0   132.4   134.8   148.0 
    ----------------------------------------------------
    x$Year: -1850
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      126.0   132.2   136.0   134.5   137.0   140.0 
    ----------------------------------------------------
    x$Year: -200
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      129.0   132.2   135.0   135.5   138.8   144.0 
    ----------------------------------------------------
    x$Year: 150
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      126.0   132.2   137.0   136.2   139.0   147.0 
    
    > by ( x$BH, x$Year, summary)
    x$Year: -4000
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      121.0   131.2   134.0   133.6   136.0   143.0 
    ----------------------------------------------------
    x$Year: -3300
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      124.0   129.2   133.0   132.7   136.0   145.0 
    ----------------------------------------------------
    x$Year: -1850
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      123.0   131.0   133.5   133.8   137.0   145.0 
    ----------------------------------------------------
    x$Year: -200
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      120.0   130.0   132.0   132.3   135.8   142.0 
    ----------------------------------------------------
    x$Year: 150
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      120.0   126.0   130.0   130.3   135.0   138.0 
    
    > by ( x$BL, x$Year, summary)
    x$Year: -4000
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      89.00   95.00  100.00   99.17  102.80  114.00 
    ----------------------------------------------------
    x$Year: -3300
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      90.00   97.00   98.50   99.07  101.80  107.00 
    ----------------------------------------------------
    x$Year: -1850
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      87.00   92.25   96.00   96.03   99.75  106.00 
    ----------------------------------------------------
    x$Year: -200
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      86.00   91.25   94.50   94.53   97.75  107.00 
    ----------------------------------------------------
    x$Year: 150
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
       81.0    91.0    94.0    93.5    97.0   103.0 
    
    > by ( x$NH, x$Year, summary) 
    x$Year: -4000
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      44.00   49.00   50.00   50.53   53.00   56.00 
    ----------------------------------------------------
    x$Year: -3300
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      45.00   48.00   50.50   50.23   52.75   56.00 
    ----------------------------------------------------
    x$Year: -1850
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      45.00   48.25   50.00   50.57   52.75   60.00 
    ----------------------------------------------------
    x$Year: -200
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      46.00   50.25   52.00   51.97   53.75   60.00 
    ----------------------------------------------------
    x$Year: 150
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      44.00   48.25   52.00   51.37   54.00   58.00
    
    > plotMeans( x$MB, x$Year, error.bars="se") 
    > plotMeans( x$BH, x$Year, error.bars="se") 
    > plotMeans( x$BL, x$Year, error.bars="se") 
    > plotMeans( x$NH, x$Year, error.bars="se") 
    
    
    
    
    
    
    
    


    c. グラフの観察

    年代の変化とともに4つの部位の寸法の平均値が変化していることが観察できる.

    特に−4000年代の指標が何時大きく変わるのかに注意を払う.



  4. 頭蓋骨 manova 関数による多変量分散分析


    頭蓋骨の4つの部位の測定値をまとめて検討する.

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

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


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

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


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

    result <- manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year) #   と記述し,@多変量分散分析を行う.

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

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

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

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

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


    > summary(result, test="Pillai") 
               Df Pillai approx F num Df den Df    Pr(>F)    
    x$Year      4 0.3533   3.5120     16    580 4.675e-06 ***
    Residuals 145                                            
    ---
    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$Year      4.00 0.6636   3.9009  16.00 434.45 7.01e-07 ***
    Residuals 145.00                                           
    ---
    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$Year      4           0.4818   4.2310     16    562 8.278e-08 ***
    Residuals 145                                                      
    ---
    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$Year      4  0.4251  15.4097      4    145 1.588e-10 ***
    Residuals 145                                             
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    

    分析結果

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

    年代によって頭蓋骨の4つの部位の形状の平均値は異なるといえる.



  5. 頭蓋骨 部位ごとの(1測定値ごと)一元配置分散分析



    年代間の平均値に差があるかを,頭蓋骨の部位ごとに検討する.これは測定値ごとの一元配置の分散分析である.

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

    summary.aov(result)

    > summary.aov(result) 
    
     Response 1 :測定値 MB
                 Df  Sum Sq Mean Sq F value    Pr(>F)    
    x$Year        4  502.83  125.71  5.9546 0.0001826 ***
    Residuals   145 3061.07   21.11                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    
     Response 2 :測定値 BH
                 Df Sum Sq Mean Sq F value  Pr(>F)  
    x$Year        4  229.9    57.5  2.4474 0.04897 *
    Residuals   145 3405.3    23.5                  
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    
     Response 3 :測定値 BL
                 Df Sum Sq Mean Sq F value    Pr(>F)    
    x$Year        4  803.3   200.8  8.3057 4.636e-06 ***
    Residuals   145 3506.0    24.2                      
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    
     Response 4 :測定値 NH
                 Df  Sum Sq Mean Sq F value Pr(>F)
    x$Year        4   61.20   15.30   1.507 0.2032
    Residuals   145 1472.13   10.15  
    
    一元配置の分散分析表が表示される.これを整形すると,以下のようになる.
        
     Response 1 :測定値 MB
    -------------------------------------------------------------
                (変動)  (不偏分散) (分散比)
    変動要因   自由度  平方和   平均平方和    F値      P値    
    -------------------------------------------------------------
    年代         4      502.83     125.71   5.9546  0.0001826***
    誤差       145     3061.07     21.11                      
    -------------------------------------------------------------
    *** p<.001  ** p<.01   * p<.05
    
     Response 2 :測定値 BH
    -------------------------------------------------------------
                (変動)  (不偏分散) (分散比)
    変動要因   自由度  平方和   平均平方和    F値      P値    
    -------------------------------------------------------------
    年代         4       229.9     57.5     2.4474  0.04897*
    誤差       145      3405.3     23.5                  
    -------------------------------------------------------------
    *** p<.001  ** p<.01   * p<.05
    
     Response 3 :測定値 BL
    -------------------------------------------------------------
                (変動)  (不偏分散) (分散比)
    変動要因   自由度  平方和   平均平方和    F値      P値    
    -------------------------------------------------------------
    年代         4       803.3     200.8    8.3057  4.636e-06***
    誤差       145      3506.0      24.2                       
    -------------------------------------------------------------
    *** p<.001  ** p<.01   * p<.05
    
     Response 4 :測定値 NH
    -------------------------------------------------------------
                (変動)  (不偏分散) (分散比)
    変動要因   自由度  平方和   平均平方和    F値      P値    
    -------------------------------------------------------------
    年代         4       61.20     15.30    1.507  0.2032
    誤差       145     1472.13     10.15  
    -------------------------------------------------------------
    *** p<.001  ** p<.01   * p<.05
    
    分析結果

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

    (1) 年代と測定値MBは,P値<.05であるから,帰無仮説は棄却する.頭蓋骨のMBの平均値は年代間で違いがある.

    (2) 年代と測定値BHは,P値<.05であるから,帰無仮説は棄却する.頭蓋骨のBHの平均値は年代間で違いがある.

    (3) 年代と測定値BLは,P値<.05であるから,帰無仮説は棄却する.頭蓋骨のBLの平均値は年代間で違いがある.

    (4) 年代と測定値NHは,P値>.05であるから,帰無仮説は棄却できない.頭蓋骨の鼻の高さNHの平均値は年代差が見られない.



  6. 頭蓋骨の多重比較 どの年代間に差があるか,[-4000]を軸に2つの年代間ごとを比較検討する


    # 年代"-4000","-3300"間の平均値の差の検定
    summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","-3300")))

    # 年代"-4000","-1850"間の平均値の差の検定
    summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","-1850")))

    # 年代"-3300","-1850"間の平均値の差の検定
    summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-3300","-1850")))

    # 年代"-4000","-200"間の平均値の差の検定
    summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","-200")))

    # 年代"-4000","150"間の平均値の差の検定
    summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","150")))


    # 年代"-4000","-3300"間の平均値の差の検定
     > summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","-3300")))
              Df  Pillai approx F num Df den Df Pr(>F)
    x$Year     1 0.02767  0.39135      4     55  0.814
    Residuals 58                                      
    
    
    # 年代"-4000","-1850"間の平均値の差の検定
    > summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","-1850")))
              Df Pillai approx F num Df den Df  Pr(>F)  
    x$Year     1 0.1876   3.1744      4     55 0.02035 *
    Residuals 58                                        
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    
    # 年代"-3300","-1850"間の平均値の差の検定
    > summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-3300","-1850")))
              Df Pillai approx F num Df den Df  Pr(>F)  
    x$Year     1 0.1898   3.2203      4     55 0.01908 *
    Residuals 58                                        
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    
    # 年代"-4000","-200"間の平均値の差の検定
    > summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","-200")))
              Df Pillai approx F num Df den Df    Pr(>F)    
    x$Year     1 0.3030   5.9766      4     55 0.0004564 ***
    Residuals 58                                            
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    
    # 年代"-4000","150"間の平均値の差の検定
    > summary(manova(cbind(x$MB, x$BH, x$BL, x$NH) ~ x$Year, subset=x$Year %in% c("-4000","150")))
              Df Pillai approx F num Df den Df    Pr(>F)    
    x$Year     1 0.3618   7.7956      4     55 4.736e-05 ***
    Residuals 58                                            
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    

    分析結果

    有意水準5%で検定を行う.

    [−4000年と−3300年]には違いが無い.P値>.05であるから,帰無仮説は棄却できない.

    [−4000年と−1850年]には違いがある.P値<.05であるから,帰無仮説は棄却する.

    [−3300年と−1850年]には違いがある.P値<.05であるから,帰無仮説は棄却する.

    [−4000年と−200年]には違いがある.P値<.05であるから,帰無仮説は棄却する.

    [−4000年と150年]には違いがある.P値<.05であるから,帰無仮説は棄却する.

    まとめると,−3300年と−1850年代を境にして頭蓋骨の測定値は年代差が見られる.

    −3300年〜−1850年代に頭蓋骨の形状の変容があり,この時代間に頭蓋骨の形状を変える他人種との交流があったことが推理される.

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