t検定

(水準数2のカテゴリ変数が連続変数に与える影響)

t検定は、2群の平均値の差を検定する方法である。

ある実験処理を施した処理区と、処理をしなかった対照区との間で、ある連続変数(体重など)を比較したい、といった状況でよく使われる。

ところで、これは「処理」をカテゴリ変数と解釈すれば、処理区と対照区という2つの水準をもつカテゴリ変数「処理」が、連続変数に与える影響を検証する手法であるとも言い換えられる。

要点

例題1

実験

ある植食性昆虫の幼虫の成長に栄養の添加が与える影響を調べるため、幼虫を何も処理を与えない餌を与えるグループ(Control)と栄養を添加した餌を与えるグループ(Treated)に分けて飼育し、その後の体重を測定したところ、以下のようなデータが得られた(単位はmg)。

データファイル(csv)

> larval.data

      Weight Treatment

1   737.3546   Control

2   818.3643   Control

3   716.4371   Control

4   959.5281   Control

5   832.9508   Control

6   717.9532   Control

7   848.7429   Control

8   873.8325   Control

9   857.5781   Control

10  769.4612   Control

11 1502.3562   Treated

12 1277.9686   Treated

13 1075.7519   Treated

14  757.0600   Treated

15 1424.9862   Treated

16 1191.0133   Treated

17 1196.7619   Treated

18 1388.7672   Treated

19 1364.2442   Treated

20 1318.7803   Treated

作図によってデータの傾向をつかむ

まずは作図によってデータの傾向をつかんでみよう。

par(cex=1.5, las=1, lwd=2, family="sans") #見栄えをよくするための細工

plot(Weight ~ Treatment, larval.data) #モデル式の右辺がfactor変数のときは、plot関数でも箱ひげ図が描かれる

図1

処理を施したグループの方が体重が1.5倍ほどになっているようにみえる。また、箱の高さがControlとTreatedでは異なっており、Treatedの方が分散が大きくなっていることがうかがえる(異分散性)。t検定は母集団が正規分布であることを仮定しているのだが、ここではとりあえずその仮定は満たされているものとする。

t検定の実践(Welchのt検定)

Rでは、t.test関数によってt検定を行う。

plot関数と同じように、モデル式をカッコ内に書くことによって行うことができる。

t.test(Weight ~ Treatment, larval.data)

> t.test(Weight ~ Treatment, larval.data)


  Welch Two Sample t-test


data:  Weight by Treatment

t = -6.0627, df = 11.355, p-value = 7.136e-05

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 -594.4284 -278.6690

sample estimates:

mean in group Control mean in group Treated 

             813.2203             1249.7690 


まず見るべきは、t = -6.0627, df = 11.355, p-value = 7.136e-05の部分である。t値=-6.06, 自由度=11.36であり、P値=7.1x10^(-5)なので、ControlとTreatedの体重が有意に異なっていることがわかる。したがって、餌への栄養添加によって幼虫の体重は有意に増加すると結論づけられる。結果のレポートにはt値、自由度、P値を添えればよい。


なお、t.test関数はデフォルトではWelchのt検定を行う。これは、2群の等分散を仮定しないt検定である。通常はWelchの方法を使うことが推奨されているので、設定をいじる必要は普通はない。ちなみに、等分散を仮定したt検定を行うには、引数var.equalTRUEにする。

t.test(Weight ~ Treatment, larval.data, var.equal=TRUE)

> t.test(Weight ~ Treatment, larval.data, var.equal=TRUE)


  Two Sample t-test


data:  Weight by Treatment

t = -6.0627, df = 18, p-value = 9.914e-06

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 -587.8262 -285.2712

sample estimates:

mean in group Control mean in group Treated 

             813.2203             1249.7690 


Welchのt検定の時と比べて、自由度とP値が異なっていることに注意する必要がある。

改めて作図

t検定は平均値を比較するものである。平均値を比較する場合、作図の基本は「平均値±標準誤差」である。

Rで平均値±標準誤差の図を描く工程の基本は以下である。

体重のデータの棒グラフを、上の工程を使って描いてみよう。

#1

se <- function(x) sd(x)/sqrt(length(x)) #標準誤差(SE)を計算する関数定義


#2

m <- tapply(larval.data$Weight, larval.data$Treatment, mean) #TreatmentごとにWeightの平均を計算

e <- tapply(larval.data$Weight, larval.data$Treatment, se) #TreatmentごとにWeightのSEを計算


#3-4

par(cex=1.5, las=1, lwd=2, family="sans") #見栄えをよくする調整

x <- barplot(m, ylim=c(0, 1400), xlab="Treatment", ylab="Larval weight (mg)") #棒グラフ

arrows(x, m-e, x, m+e, code=3, angle=90, length=0.1) #エラーバー

・工程1

Rには標準誤差を計算する関数が用意されていないので、これを定義する必要がある。標準誤差とは「平均値の標準偏差」のことで、データの標準偏差(不偏分散の平方根)をデータの個数の平方根で割ることで計算される。

function関数は「関数を定義する関数」で、function()のあとにスペースを挟んで定義したい関数を入力する。

・工程2

tapply関数は、カテゴリー別に関数を適用する関数で、第1引数に関数を適用したいベクトル、第2にカテゴリー、第3に適用したい関数を入力する。上の場合、1つ目でControlとTreatmentの平均値をmに、2つ目でControlとTreatmentの標準誤差をeに格納している。

・工程3-4

par関数はグラフの見栄えを調整する関数である。引数cexは文字とシンボルの大きさ、lasは軸ラベルの向き、lwdは線の太さ、familyはフォントを指定する。

barplot関数は棒グラフを描く関数で、第1引数にベクトルまたは行列を入れる。ここではControlとTreatmentそれぞれの平均値を格納したベクトルmを入れている。

またbarplot関数は棒グラフを描くと同時に、それぞれの棒の”x座標”を返す。これがエラーバーを描くのに必要となるので、ここでxに格納している。

arrows関数は矢印を上描きする関数なのだが、矢じりの角度を90度にすることでエラーバーのようなものにできる。第1、第2引数でそれぞれ矢印始点のx座標とy座標を指定し、第3、第4引数でそれぞれ矢印終点のx座標とy座標を指定する。ここでは、x座標は始点・終点ともに先ほど格納した棒のx座標xでよい。y座標は始点が平均-標準誤差(m-e)、終点が平均+標準誤差(m+e)となる。さらに、code=3とすることで始点と終点双方に矢じりが描かれ、angle=90とすることで矢じりの角度が90度になり、lengthで矢じり(エラーバーの横線)の長さを指定できる。

図2.

このような図を作成することができる。

図を保存する方法にはいくつかあるが、Graphics deviceをアクティブにして、メニューバーからファイル→別名で保存と進み、好きな形式で保存すればよい。

この工程は最初はかなりとっつきにくいかもしれないが、慣れれば大したことはないだろう。スクリプトをコピペすればとりあえず描けるし、自分なりに色々いじってみるとより理解が進むだろう。

パッケージ"gplots"を使ったエラーバー付き棒グラフ

Rでは”パッケージ”をインストールすることで、スマートフォンのアプリのように拡張機能を使用することができる。

メニューバーの「パッケージ」→「パッケージのインストール」と進むと、Secure CRAN mirrorsのリストが現れるが、これはどれでもよいので任意のmirrorをクリックする。

すると、パッケージのリストが現れるので、ここから"gplots"を選ぶ。

しばらくするとダウンロード、解凍、と自動的に行われて、gplotsによる拡張機能を使用することができるようになる。

パッケージを使用するには、library関数を用いる。

library(gplots) #パッケージの呼び出し

この命令以降、gplotsの機能を使うことができる。

gplotsには、エラーバーを最初から描いてくれる棒グラフの関数であるbarplot2が用意されている。

plot.ci=TRUEと指定し、ci.lにエラーバーの下限値、ci.uにエラーバーの上限値を指定することでエラーバーつきの棒グラフが描かれる(ciはconfidence intervalの略)。

par(cex=1.5, las=1, lwd=2, family="sans")

barplot2(m, plot.ci=TRUE, ci.l=m-e, ci.u=m+e, ylim=c(0, 1400), xlab="Treatment", ylab="Larval weight (mg)")

図3.

デフォルトではエラーバーの線の太さが1、横線の長さが0.1になる。これを直したければ、それぞれci.lwd、ci.widthを変更すればよい。

例題2

実験

ある植食性昆虫の幼虫の成長に栄養の添加が与える影響を調べたいが、植物の栄養は植物の個体によって大きく変化しうる。

そのため、まず植物を10個体用意し、同じ植物個体から1枚ずつ葉を採取して、1枚には何も処理を与えず、もう1枚には栄養を添加した。それぞれをControlグループとTreatedグループの幼虫に与えて飼育し、その後の体重を測定したところ、以下のようなデータが得られた(単位はmg)。

データファイル(csv)

> larval.data2

     Weight2 Treatment2 Plant.ID

1  1253.5344    Control        1

2   916.9530    Control        2

3   613.6278    Control        3

4   135.5900    Control        4

5  1137.4793    Control        5

6   786.5199    Control        6

7   795.1429    Control        7

8  1083.1509    Control        8

9  1046.3664    Control        9

10  978.1704    Control       10

11 1315.5982    Treated        1

12 1222.0460    Treated        2

13  612.9392    Treated        3

14  864.1743    Treated        4

15 1486.3316    Treated        5

16  790.3794    Treated        6

17 1191.3716    Treated        7

18 1554.6483    Treated        8

19 1469.1008    Treated        9

20 1136.5539    Treated       10

Plant.IDの列は、用いた植物個体を示す。IDが同じであれば、同じ植物個体を摂食した幼虫である。

作図によってデータの傾向をつかむ

par(cex=1.5, las=1, lwd=2, family="sans")

plot(Weight2 ~ Treatment2, larval.data2)

図4.

中央値は栄養添加の方が高いものの、かなりばらつきが大きいことがわかる。

t検定の実践(対応のあるt検定)

さきほどと同じように、まずはデフォルトの設定でt検定を行ってみる。

t.test(Weight2 ~ Treatment2, larval.data2)

> t.test(Weight2 ~ Treatment2, larval.data2)


  Welch Two Sample t-test


data:  Weight2 by Treatment2

t = -2.0274, df = 17.999, p-value = 0.05769

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 -589.82623   10.50456

sample estimates:

mean in group Control mean in group Treated 

             874.6535             1164.3143 


P=0.058であり、微妙に有意ではない。したがって、栄養添加は幼虫の成長を促進すると結論を下すことはできない… というのは早計である。


同じ植物個体の葉は栄養価が同じであり、別の個体では栄養価が違うなら、栄養添加の効果を検証するためには、同じ植物個体からの葉同士で幼虫の成長を比較する必要があると考えられる。このような実験デザインを「対応のあるデザイン」と呼び、対応のあるt検定を用いるのが有効である(「局所管理」を参照)。

対応のあるt検定を行うには、t.test内でpaired=TRUEと指定すればよい。また、比較したいデータ同士が、1つ目のデータは植物1、2つ目のデータは植物2…というように同じ順序で並んでいる必要がある。

larval.data2 #データの順番の確認

#ControlとTreatedでは、Plant.IDが同じ順序で並んでいる


t.test(Weight2 ~ Treatment2, larval.data2, paired=TRUE)

> larval.data2 #データの順番の確認

     Weight2 Treatment2 Plant.ID

1  1253.5344    Control        1

2   916.9530    Control        2

3   613.6278    Control        3

4   135.5900    Control        4

5  1137.4793    Control        5

6   786.5199    Control        6

7   795.1429    Control        7

8  1083.1509    Control        8

9  1046.3664    Control        9

10  978.1704    Control       10

11 1315.5982    Treated        1

12 1222.0460    Treated        2

13  612.9392    Treated        3

14  864.1743    Treated        4

15 1486.3316    Treated        5

16  790.3794    Treated        6

17 1191.3716    Treated        7

18 1554.6483    Treated        8

19 1469.1008    Treated        9

20 1136.5539    Treated       10

> #ControlとTreatedでは、Plant.IDが同じ順序で並んでいる

> t.test(Weight2 ~ Treatment2, larval.data2, paired=TRUE)


  Paired t-test


data:  Weight2 by Treatment2

t = -3.9115, df = 9, p-value = 0.003556

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 -457.1801 -122.1416

sample estimates:

mean of the differences 

              -289.6608 


t値=-3.91, 自由度=9, P値=0.004であり、ControlとTreatedの間での体重の差は有意であることがわかった。

このように、データに対応関係がある場合には、対応のあるt検定が有効である。

ランダム効果入りANOVAと対応のあるt検定

おもしろいことに、実験してみると、ランダム効果入りANOVAと対応のあるt検定の(P値の)結果は一致する(ようだ)。

> library(lme4)

> library(car)

> Anova(lmer(Weight2~Treatment2+(1|Plant.ID), larval.data2), test="F")

Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)


Response: Weight2

              F Df Df.res   Pr(>F)   

Treatment2 15.3  1      9 0.003556 **

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

'18 2/2