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

12.9 時系列共分散分析

(1) 単純な時系列共分散分析

周期回帰分析と同様に、周期共分散分析も周期関数を一般的な関数にして一般化することができます。 その手法を時系列共分散分析(time series analysis of covariance)と呼ぶことにしましょう。 例として表12.6.1の平均値に最も簡単な一次関数つまり直線を用いた時系列共分散分析を適用すると次のようになります。 (注1)

表12.9.1 時系列共分散分析表
要因平方和自由度平均平方和(分散)F値
群差30.8802130.88020.13273
共通回帰2167.1412167.149.31489
修正群差30.8802130.88020.13273
全体回帰2167.1412167.149.31489
非平行性1145.1611145.164.92215
残差10236.844232.654
全体13579.947
○A群
群別時系列回帰式:y=109.66 + 1.6763t
共通時系列回帰式:y=117.775 + 0.970688t
群別時系列回帰式の寄与率:r2=0.288(28.8%)
○B群
群別時系列回帰式:y=127.493 + 0.265072t
共通時系列回帰式:y=119.379 + 0.970688t
群別時系列回帰式の寄与率:r2=0.035(3.5%)
○群差:2群の24時間平均値の差の検定
FA=0.13273(p=0.7174)<F(1,44,0.05)=4.062 … 有意水準5%で有意ではない
○共通回帰:2群の時系列回帰直線が平行と仮定した時の回帰の検定
Fβc=9.31489(p=0.0038)>F(1,44,0.05)=2.062 … 有意水準5%で有意
○修正群差:2群の時系列回帰直線が平行と仮定した時の24時間平均値の差の検定(レベルの検定)
FAA=0.13273(p=0.7174)<F(1,44,0.05)=4.062 … 有意水準5%で有意ではない
○全体回帰:2群を合わせて時系列回帰式を計算した時の回帰の検定
FβT=9.31489(p=0.0038)>F(1,44,0.05)=2.062 … 有意水準5%で有意
○非平行性:2群の時系列回帰直線が平行かどうかの検定(パターンの検定)
FD=4.92215(p=0.0317)>F(1,44,0.05)=2.062 … 有意水準5%で有意
図12.9.1 A群とB群の時系列回帰直線

表12.9.1および図12.9.1と第6節の表12.6.2および図12.6.1を比較すると、時系列回帰直線の適合度が低いため残差分散が大きくなり、検定結果の有意確率が大きくなっています。 表12.6.1のデータは日内変動を検討するためのものなので、適合度が低いのは致し方ありません。 この手法を用いると、時系列で繰り返し測定された血圧の平均的なレベルと時間的なトレンドを分離して群間比較することができます。 そのため繰り返し測定データの評価方法のひとつとして利用することができます。 (→4.3 繰り返しのある多標本・多時期の計量値 (5) 繰り返し測定データの評価方法)

(2) 二元配置型時系列共分散分析

次に表12.6.1を被験者と時期を要因にした二元配置型データと捉えて、直線を用いた二元配置型時系列共分散分析を適用すると次のようになります。 (注2)

表12.9.2 二元配置型時系列共分散分析表
要因平方和自由度平均平方和(分散)F値
群差74.1125174.11250.384155
個体残差578.7713192.924
個体652.8834163.2211.437
全体回帰3957.0313957.0334.8286
非平行性2748.3712748.3724.1904
ズレ合計22733.144516.6624.5475
残差7839.469113.614
全体37930.8119
○A群
群別時系列回帰式:y=109.66 + 1.6763t
共通時系列回帰式:y=119.397 + 0.829565t
平均値の変動に対する群別時系列回帰式の寄与率:r2=0.288(28.8%)
○B群
群別時系列回帰式:y=127.493 + 0.265072t
共通時系列回帰式:y=121.002 + 0.829565t
平均値の変動に対する群別時系列回帰式の寄与率:r2=0.035(3.5%)
○群差=修正群差:2群の24時間平均値の差の検定(レベルの検定)
FA=0.384155(p=0.5793)<F(1,3,0.05)=10.128 … 有意水準5%で有意ではない
○個体:被験者ごとの24時間平均値の差の検定
FSUB=1.43662(p=0.2312)<F(4,69,0.05)=2.505 … 有意水準5%で有意ではない
○全体回帰=共通回帰:2群を合わせて時系列回帰式を計算した時の回帰の検定
FβT=34.8286(p=1.214×10-7)>F(1,69,0.05)=3.980 … 有意水準5%で有意
○非平行性:2群の時系列回帰直線が平行かどうかの検定(パターンの検定)
FD=24.1904(p=5.686×10-6)>F(1,69,0.05)=3.980 … 有意水準5%で有意
○ズレ合計:2群の時間ごとの平均値と時系列回帰直線のズレの検定
FLOF=4.5475(p=1.076×10-8)>F(44,69,0.05)=1.552 … 有意水準5%で有意

単純な時系列共分散分析の結果と比較すると、群別時系列回帰式と寄与率は同じですが、回帰と非平行性の分散比が非常に大きくなっています。 そしてこの手法では残差からズレ合計を分離して検定することができます。 また表12.9.2と第6節の表12.6.3を比較すると、全体回帰と非平行性とズレ合計以外は同じであることがわかります。 この手法は周期回帰曲線の代わりに直線を当てはめただけなので、変わるのは回帰に関係した部分だけなのです。

(3) 一元配置型時系列共分散分析

次は表12.6.1を測定時点ごとに独立した被験者で測定された一元配置型データと捉えて、直線を用いた一元配置型時系列共分散分析を適用すると次のようになります。 (注3)

表12.9.3 一元配置型時系列共分散分析表
要因平方和自由度平均平方和(分散)F値
群差74.1125174.11250.634
共通回帰3957.0313957.0333.8442
修正群差74.1125174.11250.634
全体回帰3957.0313957.0333.8442
非平行性2748.3712748.3723.5067
ズレ合計22733.144516.6624.41897
時期29512.647627.9285.37063
残差8418.1772116.919
全体37930.8119
○A群
群別時系列回帰式:y=109.66 + 1.6763t
共通時系列回帰式:y=119.397 + 0.829565t
平均値の変動に対する群別時系列回帰式の寄与率:r2=0.288(28.8%)
○B群
群別時系列回帰式:y=127.493 + 0.265072t
共通時系列回帰式:y=121.002 + 0.829565t
平均値の変動に対する群別時系列回帰式の寄与率:r2=0.035(3.5%)
○群差:2群の24時間平均値の差の検定
FA=0.634(p=0.4286)<F(1,72,0.05)=3.974 … 有意水準5%で有意ではない
○共通回帰:2群の時系列回帰曲線が平行と仮定した時の回帰の検定
Fβc=33.8442(p=1.526×10-7)>F(1,72,0.05)=3.974 … 有意水準5%で有意
○修正群差:2群の時系列回帰曲線が平行と仮定した時の24時間平均値の差の検定(レベルの検定)
FAA=0.634(p=0.4286)<F(1,72,0.05)=3.974 … 有意水準5%で有意ではない
○全体回帰:2群を合わせて時系列回帰式を計算した時の回帰の検定
FβT=33.8442(p=1.526×10-7)>F(1,72,0.05)=3.974 … 有意水準5%で有意
○非平行性:2群の時系列回帰直線が平行かどうかの検定(パターンの検定)
 FD=23.5067(p=6.958×10-6)>F(1,72,0.05)=3.974 … 有意水準5%で有意
○ズレ合計:2群の時間ごとの平均値と時系列回帰直線のズレの検定
FLOF=4.41897(p=1.255×10-8)>F(44,72,0.05)=1.545 … 有意水準5%で有意
○時期:平均値の時期変動の検定
FP=5.37063(p=1.120×10-10)>F(47,72,0.05)=1.535 … 有意水準5%で有意

二元配置型時系列共分散分析の結果と比較すると、群別および共通時系列回帰式と寄与率は同じですが、回帰の分散比も非平行性の分散比もズレの分散比も少し小さくなっています。 そして表12.9.3と第6節の表12.6.4を比較すると、共通および全体回帰と非平行性とズレ合計以外は同じであることがわかります。 この手法は周期回帰曲線の代わりに直線を当てはめただけなので、変わるのは回帰に関係した部分だけなのです。

ちなみにこの手法と同じ原理を用い、時期の代わりに用量を用いたものが用量反応解析で用いられる平行線検定法です。 平行線検定法については第13章で説明します。 (→13.2 平行線検定法)

また繰り返し測定混合効果モデルを利用した分散分析も、この手法とほぼ同じ原理を用いた手法です。 そのためどちらの手法も繰り返し測定データの評価方法のひとつとして利用することができます。 (注4) (→4.3 繰り返しのある多標本・多時期の計量値 (5) 繰り返し測定データの評価方法)


(注1) 表12.6.1の群ごとの平均値を一般化すると次のようになります。 時期t1〜t(pi)は全ての群で同一でもかまいませんし、群ごとに異なっていてもかまいません。

表12.9.4 時系列共分散分析の
一般的データ
時期平均
t1tjt(pi)
1y11y1jy1(p1)m1.
:::::
iyi1yijyi(pi)mi.
:::::
aya1yajya(pa)ma.
全体m.1m.jm.(pi)mT
※piは群ごとに異なっていても良い

時系列共分散分析ではデータyijに対して次のような時系列回帰モデルを当てはめて考えます。 これは重回帰モデルに相当し、直線を用いた時系列共分散分析では説明変数をひとつにして、x1=tとします。 2番目以降の説明変数に第2節の(注1)で説明した周期成分を入れてx2=cos(ω1t)、x3=sin(ω1t)、…、x2m=cos(ωmt)、x2m+1=sin(ωmt)とすると、時間的なトレンドと周期変動を組み合わせた時系列回帰モデルになります。 そのような時系列回帰モデルは、例えば長期間に渡る血圧の変動をトレンドと日内変動に分離して分析する時などに用いることができます。

○群別時系列回帰モデル:群ごとに計算した時系列回帰式

:群別時系列回帰式による推定値 (i=1,…,a、j=1,…,pi)
b0i:Ai群の定数、群によって異なる (i=1,…,a)
bki:Ai群の偏回帰係数、群によって異なる (k=1,…,q、i=1,…,a)
○共通時系列回帰モデル:各群の時系列回帰直線が平行と仮定して計算した時系列回帰式

:共通時系列回帰式による推定値 (i=1,…,a、j=1,…,pi)
bc0i:Ai群の定数、群によって異なる (i=1,…,a)   bck:全群共通の偏回帰係数 (k=1,…,q)
○全体時系列回帰モデル:全群を合わせて計算した時系列回帰式

:全体時系列回帰式による推定値 (j=1,…,pi)
bT0:全体の定数  bTk:全体の偏回帰係数 (k=1,…,q)

以上の時系列回帰モデルに共分散分析の原理を適用します。 まず3通りの時系列回帰モデルによる推定値を用いてデータyijを3通りに分解し、その基本式に対応する平方和と自由度を求めると次のようになります。

群別時系列回帰:
全体時系列回帰:
共通時系列回帰:
○全体変動:(yij - mT)
総例数:   平方和:
自由度:φT=N - 1
○群別回帰変動:
平方和:

        
i=[i'i]-1i'i
     

dxjki=xjk - mx.ki   
自由度:
○共通回帰変動:
平方和:


  
自由度:φbc=q
○全体回帰変動:
平方和:
     
  
自由度:φbT=q
○群差変動:(mi. - mT)
平方和:   自由度:φA=a - 1
○修正群差変動:
平方和:   自由度:φAA=a - 1
○非平行性変動:
平方和:   自由度:
○残差変動:
平方和:   自由度:

表12.6.1のA群とB群の平均値について、x1=tだけを入れた時系列回帰式を用いて実際に計算してみましょう。

○A群の群別時系列回帰式
     
○B群の群別時系列回帰式
     
  
ST=821533 - 48×129.742≒13579.9   φT=48 - 1=47
bc1=2300-1×2232.58≒0.970688
Sbc=0.970688×2232.58≒2167.14  φbc=1

S1y=2232.58  SbT=0.970688×2232.58≒2167.14  φbT=1
SA=24×128.9382 + 24×130.5422 - 48×129.742≒30.8802   φA=1
SAA=30.8802 + 2167.14 - 2167.14=30.8802   φAA=1
SD=3312.3029 - 2167.14≒1145.16  φD=2 - 1=1
SR=13579.9 - 30.8802 - 2167.14 - 1145.16≒10236.8   φR=48 - 2×(1+1)=44

これらの統計量を用いて表12.9.1の時系列共分散分析表を作成することができます。

(注2) 二元配置型時系列共分散分析を適用する一般的データは第6節(注3)の表12.6.7と同じであり、基本式も同様に記述することができます。 この手法では説明変数に時期だけでなく他の項目も入れることにより、それらを共変数にした共分散分析を行うことができます。 例えば年齢などの背景因子を入れると、背景因子で補正して時期変動を群間比較することができます。

表12.9.5 二元配置型時系列共分散分析の一般的データ
被験者時期平均
t1tjtp
11y111y1j1y1p1m1.1
:::::
n1y11(n1)y1j(n1)y1p(n1)m1.(n1)
平均m11.m1j.m1p.m1..
::::::
i1yi11yij1yip1mi.1
:::::
niyi1(n1)yij(n1)yip(n1)mi.(n1)
平均mi1.mij.mip.mi..
::::::
a1ya11yaj1yap1ma.1
:::::
naya1(na)yaj(na)yap(na)ma.(na)
平均ma1.maj.map.ma..
全体Nm.1.m.j.m.p.mT
被験者と時期の二元配置:
群と時期の二元配置:
群別時系列回帰:
全体時系列回帰:
○全体変動:(yijl - mT)
総例数:   平方和:   自由度:φT=Np - 1
○被験者変動:(mi.l - mT)
平方和:   自由度:φSUB=N - 1
○時期変動:(m.j. - mT)
平方和:   自由度:φP=p - 1
○被験者×時期変動:{yijl - (mi.l + m.j. - mT)}
平方和:
自由度:φS×PT - φSUB - φPSUB×φP=(N-1)(p-1)
○群差変動:(mi.. - mT)
平方和:   自由度:φA=a - 1
○被験者残差変動:(mi.l - mi..)
平方和:   自由度:φSRSUB - φA=N - a
○群−時期変動:(mij. - mT)
平方和:   自由度:φAP=ap - 1
○群×時期変動:{mij. - (mi.. + m.j. - mT)}
平方和:
自由度:φA×PAP - φA - φPA×φP=(a-1)(p-1)
○残差=被験者残差×時期変動:{yijl - (mi.l + mij. - mi..)}
平方和:
自由度:φRSR×PT - φSUB - φPA×PSR×φP=(N-a)(p-1)
○群別回帰変動:
平方和:
自由度:
○群別時期変動:(mij. - mi..)
平方和:   自由度:
○群別ズレの変動:
平方和:
自由度:
○全体時回帰変動=共通回帰変動:
平方和:
     
  
自由度:φbTbc=q
○非平行性変動:
平方和:
自由度:

表12.6.1のデータについて、x1=tだけを入れた時系列回帰式を用いて実際に計算してみましょう。

ST=2062810 - 120×129.92≒37930.8   φT=120 - 1=119
SSUB=2025532 - 120×129.92≒652.883   φSUB=5 - 1=4
SP=5×409340 - 120×129.92≒21821.2   φP=24 - 1=23
SS×P=37930.8 - 652.883 - 21821.2=15456.72   φS×P=4×(24-1)=92
SA=2024953 - 120×129.92≒74.1125   φA=2 - 1=1
SSR=652.883 - 74.1125=578.771  φSR=5 - 2=3
SAP=2054392 - 120×129.92≒29512.6   φAP=2×24 - 1=47
SA×P=29512.6 - 74.1125 - 21821.2=7617.287   φA×P=(2-1)×(24-1)=23
SR=SSR×P=37930.8-652.883 - 21821.2 - 7617.287=7839.4
φRSR×P=(5-2)×(24-1)=69
○A群の群別時系列回帰式
     
○B群の群別時系列回帰式
     
  
  
  

S1yT=4770  SbT=0.829565×4770≒3957.03  φbT=1
SD=6705.4 - 3957.03=2748.37   φD=1×(2-1)=1

これらの統計量を用いて表12.9.2の周期共分散分析表を作成することができます。

(注3) 一元配置型時系列共分散分析の一般的データは、表12.9.5において各群の被験者が測定時期ごとに独立したものになります。 そのため群ごとに時期t1〜tpが必ずしも同一ではなく、被験者数niも測定時期ごとに必ずしも同一ではなくなります。 したがって単純な時系列共分散分析と同様に、全体回帰変動と共通回帰変動が同一になるとは限らず、群差と修正群差も同一になるとは限りません。 この手法ではデータを次のように分解し、その基本式に対する平方和と自由度を求めます。

時期の一元配置:(yijl - mT)=(mij. - mT) + (yijl - mij.)
群別時系列回帰:
全体時系列回帰:
共通時系列回帰:
○全体変動:(yijl - mT)
総例数:   平方和:
自由度:φT=N - 1
○時期変動:(mij. - mT)
平方和:
自由度:φP=P - 1   
○残差変動:(yijl - mij.)
平方和:   自由度:φRT - φP=N - P
○群差変動:(mi.. - mT)
平方和:   自由度:φA=a - 1
○群別時系列回帰変動:
平方和:
        
i=[i'i]-1i'i
     
dxjki=nij(xjk-mx.ki)   
  自由度:
○群別ズレの変動:
平方和:
自由度:
○全体回帰変動:
平方和:
     
  自由度:φbT=q
○共通回帰変動:
平方和:


  自由度:φbc=q
○修正群差変動:
平方和:   自由度:φAA=a-1
○非平行性変動:
平方和:
自由度:

表12.6.1のデータについてx1=tだけを入れた時系列回帰式を用いて実際に計算すると、群別時系列化回帰式と共通時系列回帰式は(注2)と同じ値になり、群差と修正群差の平方和と自由度は(注2)の群差の平方和と自由度と同じ値になり、共通回帰と全体回帰の平方和と自由度は(注2)の全体回帰の平方和と自由度と同じ値になり、非平行性とズレ合計と全体の平方和と自由度は(注2)と同じ値になります。 これは表12.6.1のデータが同じ被験者について繰り返し測定したものであり、時期ごとの例数が揃っているからです。

時期と残差の平方和と自由度は次のような値になり、これは一元配置分散分析独自の値です。 これらの統計量と(注2)の統計量を用いて、表12.9.3の周期共分散分析表を作成することができます。

SP=2054393.8 - 120 × 129.92=29512.6   φP=48 - 1=47
SR=37930.8 - 29512.6=8418.2  φR=120 - 48=72

(注4) 線形混合効果モデル(linear mixed effects model)は固定効果(fixed effect)だけから構成された通常の線形モデルに変量効果(random effect)を追加したものであり、次のように記述されます。

線形混合効果モデル:

データベクトル:   誤差ベクトル:
固定効果のデザイン行列:
変量効果のデザイン行列:
固定効果の偏回帰係数ベクトル:   変量効果の偏回帰係数ベクトル:
  
  
  
  
※一般にの各成分は不等分散かつ互いに独立ではない。 しかし通常は次のように等分散かつ互いに独立という仮定を置くことが多い。
  
=の時、線形混合効果モデルは通常の線形モデル――分散分析や共分散分析や重回帰分析で用いられるモデル――になる。

このモデルに最尤法を適用してβγ最良線形普遍推定量(BLUE解)を求めると、次のようになります。 (→7.1 重回帰モデル (注1)9.3 1変量の場合 (注1))

γの密度関数:
εの密度関数:
βγの尤度関数:
最尤法:ℒ(β,γ)の最大化=exp()の中身を(-2)倍したものの最小化 ← 最小2乗法の原理

Qをβγ偏微分して0と置くと
… (1)
… (2)
(1)と(2)の連立方程式を行列で表したものを線形混合効果モデル方程式と呼び、最小2乗法における正規方程式に対応するものになる。
線形混合効果モデル方程式:
(2)より
ここで


にWoodburdyの恒等式を適用 (→「ベクトルと行列」7. 逆行列)
を(1)に代入すると



-1を重みとした重み付け最小2乗法のBLUE解に相当します。 そのためつまりが決まればデータから求められることがわかると思います。 そしては次の尤度関数に最尤法を適用することによって求めます。

尤度関数:
対数尤度関数:

この対数尤度関数の値を最大にする時のが最尤解になります。 しかしこのままでは最尤解が求められないので、普通は色々な制約条件を付けて初期値を推測し、繰り返し計算によって近似解を求めます。 しかしその制約条件はどうしても作為的になり、現実には有り得ないような条件になりがちです。

そのためどうせ現実には有り得ない作為的な条件が必要なら、=かつは等分散で互いに独立という通常の線形モデルと同じ簡単な条件の方が便利です。 したがって線形混合効果モデルは理論的には厳密であるものの、あまり実用的ではなく、従来の線形モデルつまり分散分析や共分散分析や重回帰分析を用いた方が実際的ということになります。

繰り返し測定混合効果モデル(MMRM:Mixed effect Models for Repeated Measures)は表12.9.4の群を被験者にした繰り返し測定データに線形混合モデルを適用したものであり、次のように記述されます。

線形混合効果モデル:
被験者:i=1,…,n  時期:j=1,…,pi
  
        
        
        
        
  

この場合のγiは被験者iの個体差ベクトルであり、全体からの偏差ベクトルになります。 線形混合効果モデルを利用した分散分析ではデータの時期変動を直線で近似し、全体と被験者ごとに回帰直線を当てはめるのが普通です。 そして全体に当てはめた回帰直線を固定効果にし、その回帰直線の回帰係数(定数も含む)と被験者ごとの回帰直線の回帰係数の偏差を変量効果にすると、この時のモデルは次のようになります。

線形混合効果モデル:
β0:全体の回帰直線の定数(固定効果)   β1:全体の回帰直線の回帰係数(固定効果)
βi00i0:被験者iの回帰直線の定数   βi11i1:被験者iの回帰直線の回帰係数
γi0、γi1:被験者iの回帰直線の変量効果
     

これは「被験者ごとの回帰直線が全て異なる」という仮定で組み立てたモデルであり、単純な時系列共分散分析と同様の分析をすることができます。 ただし時系列共分散分析は全ての効果を固定効果として分析するため、検定誤差としてを用います。 それに対してこのモデルでは固定効果である全体の回帰係数の検定誤差にはを用い、変量効果である回帰係数の個人差つまり非平行性の検定誤差には'を用います。

また回帰直線の定数項だけを用いる、つまり全体の定数項と被験者ごとの定数項の偏差を個体差として、それを変量効果にする時もあります。 そのモデルは次のようになり、時系列共分散分析において被験者ごとの回帰直線が平行な時に相当します。

線形混合効果モデル:
β0:全体の回帰直線の定数(固定効果)   βi00i0:被験者iの回帰直線の定数   γi0:被験者iの回帰直線の変量効果

これらのモデルに固定効果として群と時期を追加すると、一元配置型時系列共分散分析に相当する手法になります。 さらに被験者ごとの時期1〜pが全て等しく、固定効果に被験者を追加すると、二元配置型時系列共分散分析に相当する手法になります。 また固定効果も変量効果も定数項だけを用いるモデルで、被験者ごとの時期1〜pが全て等しく、固定効果に被験者と時期と交互作用を追加すると、繰り返し測定型二元配置分散分析に相当する手法になります。 (→4.3 繰り返しのある多標本・多時期の計量値 (5) 繰り返し測定データの評価方法)

このように線形混合効果モデルは、色々な統計手法で用いられる数学モデルを一般化したものです。 そのため線形混合モデル用のソフトウェアを作成しておくと、そのソフトウェアを利用して色々な統計手法を行うことができます。 そこで既存の統計ソフトは線形混合効果モデルをサポートし、まるで万能ツールのように宣伝しています。 そのせいか医学論文などでも、往々にして「統計手法として線形混合効果モデルを用いた」と記述してしまうことがあるようです。

しかし上記の説明からわかるように線形混合効果モデルは数学モデルの名称であり、統計手法の名称ではありません。 線形混合効果モデル用のソフトウェアを用いた時は、具体的なモデルを記述する必要がありますし、それが従来の統計手法に相当する時はモデル名ではなく統計手法の名称を記述するべきです。

また線形混合効果モデルは、普通は最尤法を用いて繰り返し計算による近似計算によって結果を求めます。 そのため従来の統計手法に相当する時は従来の統計手法よりも結果の信頼性が低くなり、しかもその手法独自の有用な情報を出力することができません。 したがってそのような時はその手法専用のソフトウェアを用いるのが賢明です。 「何でもできる万能ツール」というものは、えてして「中途半端で何にも使えない器用貧乏ツール」になってしまうものです。