玄関雑学の部屋雑学コーナー統計学入門

8.4 共分散分析と重回帰分析

(1) ダミー変数を利用した重回帰分析

第3節で回帰式を重回帰式にした共分散分析について説明しましたが、共分散分析は説明変数の中に計量尺度のデータと名義尺度のデータが混在した重回帰分析に相当します。 例えば表8.1.1の薬剤を「0:A 1:B」というダミー変数で表し、薬剤と投与前の収縮期血圧を説明変数にし、投与前後の変化量を目的変数にした重回帰分析を適用してみましょう。 (→第7章 重回帰分析)

表8.4.1 薬剤投与前後の収縮期血圧
患者No.薬剤 (0:A,1:B)投与前投与後変化量
10140126-14
20140132-8
30145127-18
40145132-13
50150130-20
60150135-15
70155132-23
80160140-20
91160142-18
101165152-13
111165155-10
121165150-15
131170155-15
141170150-20
151170148-22
161175155-20
171175150-25
181180157-23
191180160-20
201185158-27
目的変数y(変化量):my=-18  SDy=5
説明変数x1(薬剤):mx1=0.6  SDx1=0.5   説明変数x2(投与前収縮期血圧):mx2=162  SDx2=14
単相関係数:rx1y=-0.264  rx2y=-0.605   rx1x2=0.857
重回帰式:y=59.817 + 9.484x1 - 0.514x2
偏回帰平方和:Sb1=114.81  Sb2=254.94
偏回帰係数b1の検定:Fβ1=10.553(p=0.0047)>F(1,17,0.05)=4.451 … 有意水準5%で有意
偏回帰係数b2の検定:Fβ2=23.435(p=0.0002)>F(1,17,0.05)=4.451 … 有意水準5%で有意
重寄与率:R2=0.609(60.9%)  重相関係数:R=0.780
重回帰式の検定:Fβ=13.238(p=0.0003)>F(2,17,0.05)=3.592 … 有意水準5%で有意
表8.4.2 分散分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
回帰288.0142144.00713.238
残差184.9361710.8786 
全体472.9519 

上記の重回帰式の変数x1は薬剤を表すダミー変数なので、この変数に0または1を代入した時の重回帰式はそれぞれA剤投与群またはB剤投与群のx2(投与前収縮期血圧)とy(変化量)の回帰式になります。

○A剤投与群 … x1=0を代入
回帰式:y=59.817 + 9.484×0 - 0.514x2≒60 - 0.514x2
○B剤投与群 … x1=1を代入
回帰式:y=59.817 + 9.484×1 - 0.514x2≒69 - 0.514x2

これらの回帰式は第2節で説明した共分散分析における共通回帰式に相当することがわかると思います。 また表8.2.1の共分散分析表の非平行性の項を残差にプールすると次のような共分散分析表になります。

表8.4.3 共分散分析表(非平行性を残差にプール)
要因平方和自由度平均平方和(分散)F値
群差33.08133.083.040
共通回帰254.941254.9423.435
修正群差114.811114.8110.553
全体回帰173.211173.2115.922
残差184.941710.879
全体472.9519

この共分散分析表と重回帰分析の結果を比べると、共通回帰の検定が偏回帰係数b2(投与前収縮期血圧)の検定に相当し、修正群差の検定が偏回帰係数b1(薬剤)の検定に相当していることがわかると思います。 つまりこの場合の重回帰分析は2群の回帰式が平行と仮定して共分散分析を行っていることに相当します。 またこの表の残差の平均平方和は表8.2.1の共分散分析表の残差の平均平方和よりも小さくなり、検定と推定の精度が少し高くなります。 このように非平行性のF値が1未満の時は非平行性の項を残差にプールした方が効率的です。 (注1)

次にダミー変数化した薬剤と投与前収縮期血圧を掛け合わせて「薬剤×投与前」という項目を作り、この項目まで含めて重回帰分析を適用してみましょう。

表8.4.4 薬剤投与前後の収縮期血圧
患者No.薬剤 (0:A,1:B)投与前投与後薬剤×投与前変化量
101401260-14
201401320-8
301451270-18
401451320-13
501501300-20
601501350-15
701551320-23
801601400-20
91160142160-18
101165152165-13
111165155165-10
121165150165-15
131170155170-15
141170150170-20
151170148170-22
161175155175-20
171175150175-25
181180157180-23
191180160180-20
201185158185-27
目的変数y(変化量):my=-18  SDy=5
説明変数x1(薬剤):mx1=0.6  SDx1=0.5   説明変数x2(投与前収縮期血圧):mx2=162  SDx2=14
説明変数x3(薬剤×投与前収縮期血圧):mx3=103  SDx3=86
単相関係数:rx1y=-0.264  rx2y=-0.605   rx3y=-0.302  rx1x2=0.857   rx1x3=0.998  rx2x3=0.882
重回帰式:y=62.892 + 4.405x1 - 0.535x2 + 0.032x3
偏回帰平方和:Sb1=0.174  Sb2=99.334   Sb3=0.2335
偏回帰係数b1の検定:Fβ1=0.015(p=0.9037)<F(1,16,0.05)=4.494 … 有意水準5%で有意ではない
偏回帰係数b2の検定:Fβ2=8.605(p=0.0097)>F(1,16,0.05)=4.494 … 有意水準5%で有意
偏回帰係数b3の検定:Fβ3=0.020(p=0.8887)<F(1,16,0.05)=4.494 … 有意水準5%で有意ではない
重寄与率:R2=0.609(60.9%)  重相関係数:R=0.781
重回帰式の検定:Fβ=8.323(p=0.0015)>F(3,16,0.05)=3.239 … 有意水準5%で有意
表8.4.5 分散分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
回帰288.247396.0828.323
残差184.7031611.544 
全体472.9519 

上記の重回帰式の変数x1は薬剤を表すダミー変数であり、x3はx1(薬剤)×x2(投与前収縮期血圧)です。 そのためこの2つの変数に0または1を代入した時の重回帰式は、それぞれA剤投与群またはB剤投与群のx2とy(変化量)の回帰式になります。

○A剤投与群 … x1=0を代入
回帰式:y=62.892 + 4.405×0 - 0.535x2 + 0.032×0×x2≒63 - 0.535x2
○B剤投与群 … x1=1を代入
回帰式:y=62.892 + 4.405×1 - 0.535x2 + 0.032×1×x2≒67 - 0.503x2

これらの回帰式は第2節で説明した共分散分析における群別回帰式に相当することがわかると思います。 そして表8.2.1の共分散分析表と重回帰分析の結果を比べると、非平行性の検定が偏回帰係数b3(薬剤×投与前収縮期血圧)の検定に相当していることがわかると思います。 つまりこの場合の重回帰分析は2群の回帰式が非平行と仮定して共分散分析を行っていることに相当します。 (注2)

(2) 交互作用項目の意味

この重回帰式のx3(薬剤×投与前収縮期血圧)は薬剤と投与前収縮期血圧の交互作用(interaction)を表す項目であり、投与前収縮期血圧が変化量に与える影響が薬剤によって異なっている程度を表します。 投与前収縮期血圧が変化量に与える影響が薬剤によって異なっているということは、投与前収縮期血圧と変化量の回帰式の傾きが薬剤投与群ごとに異なっている、つまり非平行であることを表し、結局のところ共分散分析における非平行性に相当します。

そしてこの交互作用は、この例題のように一方が名義尺度のデータで他方が計量尺度のデータという時に限らず、計量尺度のデータ同士、名義尺度のデータ同士でも全く同じようにして計算することができます。

普通の重回帰分析では説明変数同士に相関はあるものの交互作用はないと仮定して計算します。 しかしこの交互作用がないという仮定はある意味で暗黙の仮定であり、はっきり意識することは少ないと思います。 ところが実際には説明変数に性別のような項目が含まれていると、目的変数に与える他の説明変数の影響が男と女で異なる、つまり性と他の説明変数の間に交互作用がある時があります。 そのような時は説明変数間の交互作用を考慮するために、説明変数同士を掛け合わせた交互作用項目も含めて重回帰分析を行う必要があります。

また医学分野などでは2つの項目の積を計算したり、比を計算したりして、総合的な指標を作ることがしばしばあります。 例えば肥満度の指標としてよく用いられるBMI(Body Mass Index)は次のような計算式で求めます。

BMI= 体重(kg) 身長(m) 2
※例えば体重=60kg、身長=160cm(1.6m)の時
BMI= 60 1.6 2 =23.4375

身長を平方した値は人間の体表面積とほぼ比例するため、BMIは単位体表面積あたりの体重を表す指標になります。 つまり端的に言えば体型が球に近いほどBMIが大きくなるということです。 このBMIがある疾患の重症度相当の値に影響を与えていて、その関係が次のような直線回帰式で近似できたとします。

直線回帰式:y(重症度)=10 + 5x(BMI)

この直線回帰式は、実は重症度に関する体重と身長の平方の逆数の交互作用を表す式になります。 例えば身長が150cmの時と190cmの時の体重と重症度の直線回帰式は、それぞれ次のようになります。

○身長が150cmの時 … 身長の値として1.5を代入
y( 重症度 )=10+5 × (BMI)=10+5 × 体重 1.5 2 10+2.222 × 体重
○身長が190cmの時 … 身長の値として1.9を代入
y( 重症度 )=10+5 × (BMI)=10+5 × 体重 1.9 2 10+1.385 × 体重

一方、体重が50kgの時と80kgの時の身長の平方の逆数と重症度の直線回帰式は、それぞれ次のようなります。

○体重が50kgの時 … 体重の値として50を代入
y( 重症度 )=10+5 × (BMI)=10+5 × 50 身長 2 10+250 × 1 身長 2
○体重が80kgの時 … 体重の値として80を代入
y( 重症度 )=10+5 × (BMI)=10+5 × 80 身長 2 10+400 × 1 身長 2
図8.4.1 体重と重症度の関係 図8.4.2 身長と重症度の関係

このようにBMIと重症度の関係が直線で近似できるということは、体重と重症度の関係は直線で近似でき、身長の平方の逆数と重症度の関係も直線で近似でき、体重と身長の平方の逆数の間に交互作用があるということを意味します。 つまり身長が高くなるほど体重が重症度に与える影響は弱くなる、あるいは体重が重くなるほど身長の平方の逆数が重症度に与える影響は強くなるということになります。 逆にい言えば体重と重症度の関係が直線で近似でき、身長の平方の逆数と重症度の関係も直線で近似でき、体重と身長の平方の逆数の間に交互作用がある時に限って、BMIと重症度の関係を直線で近似することができるということになります。

これはほとんど意識されていませんが、実は大きな意味を持つことです。

(注1) 表8.4.1を一般化したものに重回帰モデルを当てはめ、その解を求めてみましょう。 (→7.3 変数の選択 (注1))

重回帰モデル:
x1:群1=0(n1例)、群2=1(n2例)、n1 + n2=r
データ行列

単純積和行列

(1,1)要素で掃き出し

(2,2)〜(3,3)要素で掃き出し

b0=my(1) - b2mx2(1)   b1=(my(2) - b2mx2(2)) - (my(1) - b2mx2(1))
my(1)、mx2(1):群1のyとx2の平均値   my(2)、mx2(2):群2のyとx2の平均値

Sx2y(1)、Sx2y(2)、Sx2x2(1)、Sx2x2(2):群1と群2の積和と平方和
SR=Syy - n1(my(1) - my)2 - n2(my(2) - my)2 - b2(Sx2y(1) + Sx2y(2))

この結果と第2節の(注1)の結果を比べると、b0が群1の共通回帰式の定数項に相当し、b1が群2と群1の共通回帰式の定数項の差に相当し、b2が共通回帰式の回帰係数に相当し、SRが非平行性をプールした残差平方和に相当することがわかると思います。 表8.4.1のデータについて実際に計算すると次のようになります。 (→8.2 共分散分析結果の解釈 (注1))




b0=59.817  b1=9.484  Sb1=114.81   SEb1=2.920
Fβ1=10.553(p=0.0047)>F(1,17,0.05)=4.451
b2=-0.5144  Sb2=254.94  SEb1=0.106
Fβ2=23.435(p=0.0002)>F(1,17,0.05)=4.451
Fβ=13.238(p=0.0003)>F(2,17,0.05)=3.592
※この時の分散分析表は本文中の表8.4.2になる。

(注2) 共分散分析は次のような手順で複数の重回帰分析を行い、その結果をまとめたものに相当します。

1.群だけを説明変数にした重回帰分析

第2節(注1)の表8.2.3のデータについて、群A1〜Aaを(a-1)個のダミー変数で表して説明変数にし、yに関する重回帰分析を行ったもの。

これは群を要因Aにした一元配置分散分析に相当し、この時の回帰平方和が共分散分析における群差の平方和SAになり、回帰自由度が群差の自由度φAになります。 表8.4.1のデータについて実際に計算すると次のようになります。

表8.4.6 分散分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
回帰33.075133.0751.35345
残差439.8751824.4375 
全体472.9519 
SA=33.075  φA=1
2.共変数だけを説明変数にした重回帰分析

第2節(注1)の表8.2.3のデータについて、共変数xだけを説明変数にし、yに関する重回帰分析を行ったもの。

これは共変数に関する全体回帰式を計算したものに相当し、この時の回帰平方和が共分散分析における全体回帰の平方和Sβ0になり、回帰自由度が全体回帰の自由度φβ0になります。 表8.4.1のデータについて実際に計算すると次のようになります。

表8.4.7 分散分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
回帰173.2071173.20710.4014
残差299.7431816.6524 
全体472.9519 
Sβ0=173.207  φβ0=1
3.群と共変数を説明変数にした重回帰分析

第2節(注1)の表8.2.3のデータについて、群を(a-1)個のダミー変数で表したものと共変数xを説明変数にし、yに関する重回帰分析を行ったもの。

これは共変数に関する全体回帰式と修正群差を計算したものに相当すると同時に、共通回帰式と群差を計算したものにも相当します。 そのためこの時の回帰平方和は共分散分析における全体回帰の平方和Sβ0と修正群差の平方和SAAの合計になると同時に、共通回帰の平方和Sβと群差の平方和SAの合計にもなります。 そして回帰自由度は全体回帰の自由度φβ0と修正群差の自由度φAAの合計になると同時に、共通回帰の自由度φβと群差の自由度φAの合計にもなります。

表8.4.1のデータについて実際に計算すると、本文中の表8.4.2のようになります。 そしてこの結果と表8.4.6と表8.4.7の結果を利用して、共通回帰と修正群差の平方和と自由度を求めることができます。

Sβ0 + SAA=Sβ + SA=288.014   φβ0 + φAAβ + φA=2
SAA=288.014 - Sβ0=288.014 - 173.207=114.807   φAA=2 - φβ0=2 - 1=1
Sβ=288.014 - SA=288.014 - 33.075=254.939   φβ=2-φA=2 - 1=1
4.群と共変数と交互作用を説明変数にした重回帰分析

第2節(注1)の表8.2.3のデータについて、群を(a-1)個のダミー変数で表したものと共変数x、そして群と共変数を掛けた交互作用変数を説明変数にし、yに関する重回帰分析を行ったもの。

これは共変数に関する群別回帰式と群差を計算したものに相当し、この時の回帰平方和が共分散分析における全体回帰の平方和Sβ0と修正群差の平方和SAAと非平行性の平方和SDの合計になり、回帰自由度が全体回帰の自由度φβ0と修正群差の自由度φAAと非平行性の自由度φDの合計になり、残差SRがそのまま残差SRになります。

表8.4.1のデータについて実際に計算すると、本文中の表8.4.5のようになります。 そしてこの結果と表8.4.2の結果を利用して、非平行性の平方和と自由度を求めることができます。

Sβ0 + SAA + SD=288.247   φβ0AA + φD=3
SD=288.247 - (Sβ0 + SAA)=288.247 - 288.014=0.233   φD=3 - (φβ0 + φAA)=3 - 2=1
SR=184.703  φR=16

以上の結果から表8.2.1の共分散分析表を作成することができます。

表8.2.1 共分散分析表(ANCOVA table)
要因平方和自由度平均平方和(分散)F値
群差33.08133.082.865
共通回帰254.941254.9422.084
修正群差114.811114.819.945
全体回帰173.211173.2115.004
非平行性0.2310.230.020
残差184.701611.54
全体472.9519

この手順を応用すれば、色々な多変量解析手法について共分散分析的な解析を行うことができます。 例えばこの手順をロジスティック回帰分析に応用すればロジスティック回帰曲線に関する共分散分析相当の解析を行うことができますし、周期回帰分析に応用すれば周期回帰曲線に関する共分散分析相当の解析を行うことができます。 (→第10章 ロジスティック回帰分析12.6 周期共分散分析)