【R×栄養統計】どの群間に差がある？対応のない多標本における多重比較

みなさん，こんにちは。
シンノユウキ（shinno1993）です。

Rで多重比較を行う方法を紹介していきます．ボリュームがあるので，前編と後編とに分けて行っています．後編の今回は，実際にRを用いて多重比較を行っていきます．

多重比較を理解するステップ

検定の多重性の問題を理解する（前回）
多重比較の使い分けについて理解する（前回）
実際に検定してみる ← 今回紹介

前回の記事はこちらです：

【R×栄養統計】多重比較を行う理由と使い分けを紹介｜検定の多重性問題

R言語で栄養統計を行う連載．検定の多重性問題と多重比較の使い分けについて紹介します．

仮説検定の使い分けについては前回も紹介していますので，参照してください．仮説検定としては，以下の３つ

ボンフェローニ補正
チューキーのHSD検定
ダネット検定

をとりあげます．

ではいきましょう！

データのインポート
ボンフェローニ補正
チューキーのHSD検定
ダネット検定
まとめ
連載目次

データのインポート

まずは以下からデータをダウンロードしてください：

[wpdm_package id='2632']

そして，変数dfとして読み込んでおきましょう：

df <- read.csv(
  "sample.csv",
  header = TRUE  
)
df

以下のようなデータが表示されればOKです：

     id  speed energy
1     1   slow   1956
2     2   slow   1654
3     3   slow   1589
4     4   slow   1847
・・・
297 297   fast   2387
298 298   fast   1882
299 299   fast   2371
300 300   fast   1949

300人のデータが収載されており，食べる速さを表したspeedと摂取エネルギーを表したenergyが入力されています．

食べる速さには，slow, medium, fast の3段階が入力されており，この3つのグループでエネルギー摂取量に差がないかどうかを検定します．

本来であれば各群の平均値などの基本統計量を見たりする必要がありますが，以前の記事で一通り行っておりますので，今回は省略します．

ボンフェローニ補正

まずはボンフェローニ補正から紹介します．以下のように書きます：

pairwise.t.test(
x = 検定の対象とするベクトル,
g = グループのベクトル,
paired = 対応の有無,
p.adjust.method = 'p値を調整する方法'
)

x : 検定の対象とするベクトルを指定します．df$energy のように差があるか知りたい変数のベクトルを指定します．
g : グループ分けの対象となるベクトルを指定します．今回のように食べる速さが違う3つのグループに分けて検定したい場合はそれを指定します．df$speed のようにベクトルで指定してください．
paired : 対応の有無を指定します．対応がある場合はTRUEを，ない場合はFALSEを指定してください．
p.adjust.method : p値を調整する方法を指定します．文字列で指定してください．今回はボンフェローニ補正を紹介していますので'bonferroni'を選択しますが，他の方法を利用することもできます．

では，今回のテストデータで実際に利用してみます．以下のようになります：

pairwise.t.test(
  x = df$energy,
  g = df$speed,
  paired = FALSE,
  p.adjust.method = 'bonferroni'
)

これを実行すると，以下のような結果が出力されます：

	Pairwise comparisons using t tests with pooled SD 

data:  df$energy and df$speed 

       fast    medium 
medium 9.1e-13 -      
slow   < 2e-16 < 2e-16

P value adjustment method: bonferroni

今回の検定では，食べる速さ（df$speed）をグループ分けの変数として，摂取エネルギー（df$energy）を比較する変数として利用しています．

そのP値は表形式で出力され，たとえばfastとmediumを比較した結果は「9.1e-13」で，fastとslow，mediumとslowを比較した結果はどちらも「< 2e-16」です．ですので，有意水準を5％とした場合，どのグループ間でも統計的に有意な差がある，という結論になります．

「9.1e-13」という表記にはあまり見覚えがないかもしれませんが，これは9.1×10の-13乗，すなわち0.00000000000091のことです．また「< 2e-16」というのは2×10の-16乗よりも小さい，ということになります．

チューキーのHSD検定

次にチューキーのHSD検定からです．以下のように書きます：

TukeyHSD(
x = aov(対象のベクトル~グループ分けベクトル),
conf.level = 信頼区間,
)

x：分散分析のモデルを指定します．これはaov関数で指定すると便利です．aov関数は(対象とする値のベクトル~グループ分けのベクトル)のように指定します．今回の場合はaov(df$energy~df$speed)のようにすると良いでしょう．
conf.level：チューキーのHSD検定では，それぞれのグループ間の平均値の差を，信頼区間でもって表示させることができます．デフォルトは0.95です．この場合，グループ間の平均値の差は95％の確率でこの範囲に収まる，という解釈ができます．

今回の目的に合わせる場合，以下のようなコードになります：

TukeyHSD(
    x = aov(df$energy~df$speed),
    conf.level = 0.95,
)

結果は以下のように出力されます：

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = df$energy ~ df$speed)

$`df$speed`
               diff       lwr       upr p adj
medium-fast -186.92 -244.5633 -129.2767     0
slow-fast   -419.68 -477.3233 -362.0367     0
slow-medium -232.76 -290.4033 -175.1167     0

diffはそれぞれのグループの平均値の差を，lwr・uprは信頼区間の下限と上限をそれぞれ示しています．p adjは調整したp値を示しています．

一例として真ん中のslow-fastについて説明すると，食べるのが遅いグループ(slow)は速いグループ(fast)よりもエネルギー摂取量が平均値で420キロカロリー程度低く，95％信頼区間は-477~-362の範囲です．調整されたp値は限りなく0に近く，統計的に有意な差がある，と解釈できます．

ダネット検定

最後はダネット検定です．以下のように書きます：

install.packages("SimComp")
library("SimComp")SimTestDiff(
data = データ,
grp = グループ分けのベクトル名,
resp = 対象とするベクトル名,
type = 'Dunnett',
base = 基準とするグループのインデックス
)

ダネット検定を使用するためにはSimCompと呼ばれるライブラリを使う必要があります．初めて使用する方はインストールしてください．その後，libray("SimComp")で読み込みます．

data：対象とするデータを選択します．dfのように，グループ分けベクトルも，検定の対象とするベクトルも含まれたデータを選択してください．
grp：グループ分けのベクトル名を指定します．'speed'のように文字列で指定してください．
resp：検定の対象とするベクトル名を指定します．'energy'のように文字列で指定してください．
type：比較に用いる方法を指定します．今回はダネット検定を使用しますので，'Dunnett'のように文字列で指定してください．
base：基準とするグループのインデックスを指定します．1や2などの整数で指定します．グループ分けのベクトルはアルファベット順で上から連番が振られます．その連番を指定します．例えばfast, medium, slowの３つのグループに分けられる場合，アルファベット順だと1.fast, 2.medium, 3.slow となります．これで1のfastを基準としたい場合は1を，3のslowを基準としたい場合は3を指定します．

今回のデータの場合，実際のコードは以下のようになります：

install.packages("SimComp")
library("SimComp")

SimTestDiff(
  data = df,
  grp = 'speed',
  resp = 'energy',
  type = 'Dunnett',
  base = 1
)

結果は以下のように出力されます：

Test for differences of means of multiple endpoints 
Assumption: Heterogeneous covariance matrices for the groups 
Alternative hypotheses: True differences not equal to the margins 
 
                 comparison endpoint margin estimate statistic degr.fr p.value.raw
medium - fast medium - fast   energy      0   -186.9    -7.536   197.9           0
slow - fast     slow - fast   energy      0   -419.7   -17.373   197.8           0
              p.value.adj
medium - fast           0
slow - fast             0

食べる速さが速い人とそれ以外の人とを比べています．estimateが平均値の差でp.value.adjが調整されたp値です．限りなく0に近いという結果になりました．

まとめ

今回は実際にRのコードを用いて多重比較を行いました．３つの方法を紹介しましたが，目的に応じて上手に使い分けてください．

データのインポート

ボンフェローニ補正

チューキーのHSD検定

ダネット検定

まとめ

連載目次