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

12.3 分散分析を応用した周期回帰分析

(1) 二元配置型周期回帰分析

第2節で説明した周期回帰分析は表12.1.1の平均値だけを利用する単純な手法です。 しかし個々の被験者のデータを全て利用すれば、より効率的かつ詳細な分析をすることができます。 このデータは第4章の表4.1.6と同じような構造をしているので、被験者と時期を要因にした二元配置型のデータと捉えることができます。 そのためこのデータに二元配置分散分析を適用すれば、被験者間変動と時期変動を検討することができるはずです。 そこで実際に二元配置分散分析を適用すると次のようになります。 (→4.1 多標本の計量値 (2)データに対応がある場合)

表12.3.1 収縮期血圧の分散分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
A(被験者)275.5211275.5212.54499
B(時期)22417.323974.6669.00301
残差2489.9823108.26 
全体25182.847 
要因A(被験者):FA = 2.545(p = 0.1243) < F(1,23,0.05) = 4.279 … 有意水準5%で有意ではない
要因B(時期):FB = 9.003(p = 7.925×10-7) > F(23,23,0.05) = 2.014 … 有意水準5%で有意

上の結果から被検者間変動つまり個人差は比較的小さく、時期変動はかなり大きいことがわかります。 そしてこの時の時期変動は時期ごとの平均値のバラツキであり、図12.2.1の元データの時期変動に相当します。 第2節で説明した単純な周期回帰分析は、元データつまり時期ごとの平均値と周期回帰曲線とのズレを回帰誤差と考え、それに基づいて周期変動を分析する手法です。 しかし二元配置分散分析では時期ごとの平均値に対する個々の被験者のデータのバラツキつまり残差を誤差と考え、それに基づいて時期ごとの平均値の変動を分析しています。

そこでこの二元配置分散分析と周期回帰分析を組み合わせると、平均値の時期変動を周期回帰曲線による周期変動と周期回帰曲線と平均値とのズレに分離し、それらを残差に基づいて個別に分析することができるはずです。 その手法を二元配置型周期回帰分析(two-way layout periodic regression analysis)と呼ぶことにします。 この二元配置型周期回帰分析を表12.1.1のデータに適用すると次のようになります。 (注1)

表12.3.2 収縮期血圧の二元配置型周期回帰分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
周期回帰15536.327768.1471.7545
ズレ6881.0321327.6683.02668
時期22417.323974.6669.00301
被験者275.5211275.5212.54499
残差2489.9823108.26 
全体25182.847 
○周期回帰
周期回帰式:y = 128.9 + 25.4・cos(15t - 210.5)
重寄与率:R2 = 0.693(69.3%)
Fβ = 71.755(p = 1.2982×10-10) > F(2,23,0.05) = 3.422 … 有意水準5%で有意
○ズレ:周期回帰曲線と平均値のズレ(LOF:Lack Of Fit)
FLOF = 3.027(p = 0.0057) > F(21,23,0.05) = 2.036 … 有意水準5%で有意
○時期:二元配置分散分析における要因B
FP = 9.003(p = 7.925×10-7) > F(23,23,0.05) = 2.014 … 有意水準5%で有意
○被験者:二元配置分散分析における要因A
FSUB = 2.545(p = 0.1243) < F(1,23,0.05) = 4.279 … 有意水準5%で有意ではない

表12.3.2の二元配置型周期回帰分析表では、「周期回帰」の平方和と「ズレ」の平方和を合計したものが「時期」の平方和になり、「周期回帰」の自由度と「ズレ」の自由度を合計したものが「時期」の自由度になります。 これは「時期」が平均値の変動を表し、「周期回帰」が周期回帰曲線による周期変動を表し、「ズレ」がそれらの変動の差になるからです。 また普通の分散分析と同様に、この手法における検定はほとんど有意性検定になります。 そのため検定結果よりも周期回帰式そのものや重寄与率を科学的に検討する方が有意義です。

この場合の周期回帰式は周期成分1だけを用いたものであり、平均値だけを利用した単純な周期回帰分析と同じものになります。 また重寄与率は平均値の変動に対する周期変動の割合であり、やはり平均値だけを利用した単純な周期回帰分析と同じ値になります。 ちなみに個々のデータの変動に対する周期変動の重寄与率はR2 = 15536.3/25182.8 = 0.617(61.7%)になります。 しかし周期回帰曲線の適合度の指標としては、この値よりも平均値の変動に対する重寄与率の方が合理的です。

上の結果より周期回帰曲線は平均値の変動とかなり適合しているものの、多少のズレがあることがわかります。 このことは図12.2.1を見れば納得できると思います。 そこで今度は周期成分2を追加して二元配置型周期回帰分析を行うと次のようになります。

表12.3.3 周期成分2を追加した二元配置型周期回帰分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
周期回帰18147.444536.8641.9071
ズレ4269.8819224.7312.07584
時期22417.323974.6669.00301
被験者275.5211275.5212.54499
残差2489.9823108.26 
全体25182.847 
○周期回帰
周期回帰式:y = 128.9 + 25.4・cos(15t - 210.5) + 10.4・cos(30t - 262.6)
重寄与率:R2 = 0.810(81.0%)
Fβ = 41.907(p = 3.0446×10-10) > F(4,23,0.05) = 2.796 … 有意水準5%で有意
○ズレ
FLOF = 2.076(p = 0.0483) > F(19,23,0.05) = 2.061 … 有意水準5%で有意
○時期
FP = 9.003(p = 7.925×10-7) > F(23,23,0.05) = 2.014 … 有意水準5%で有意
○被験者
FSUB = 2.545(p = 0.1243) < F(1,23,0.05) = 4.279 … 有意水準5%で有意ではない
図12.3.1 平均値の変動と周期回帰曲線

図12.3.1の周期回帰曲線と平均値の変動は図12.2.2と同じですが、平均値の誤差範囲を示すために平均値の上下に標準誤差を描いてあります。 このグラフから周期回帰曲線が平均値の誤差範囲を外れる部分が比較的少ないことがわかると思います。 またこの時の周期回帰式と平均値の変動に対する重寄与率は、第2節で周期成分1と周期成分2を用いたものと同じになります。

表12.3.2と表12.3.3の平方和と自由度の欄を比べると、表12.3.3ではズレの平方和と自由度が小さくなり、その分だけ周期回帰の平方和と自由度が大きくなっていることがわかります。 これは表12.3.2ではズレの変動に含まれていた周期成分2の変動分が周期回帰の方に移行したためです。 その結果、ズレの平均平方和Ms(分散V)と分散比Fが小さくなり、ズレの有意確率が大きくなっています。 これは周期成分2の変動分が周期回帰に移行することによってズレが小さくなり、ズレが0である確率が増えたことを表しています。

一方、周期回帰の平均平方和と分散比も小さくなり、周期回帰の有意確率が少し大きくなっています。 これは周期成分2の変動分が周期回帰の中では相対的に小さい、つまり周期成分1の変動分よりも小さく、周期成分2を追加することによって周期回帰式の信頼性がかえって低くなってしまったことを表しています。 しかし有意確率の増加分はわずかであり、周期回帰式の信頼性はまだまだ非常に高いと言えます。

以上のことから表12.1.1のデータについては、周期成分1だけの周期回帰式よりも周期成分1と周期成分2を組み合わせた周期回帰式の方がより合理的なモデルである可能性が高いことがわかると思います。 このように二元配置型周期回帰分析では周期回帰の分散比とズレの分散比のバランスを見て、できるだけ少ない周期成分で、できるだけ信頼性の高い合理的な周期回帰式を検討することが可能です。 これがこの手法の特徴であり、単純な周期回帰分析よりも優れた点です。

(2) 一元配置型周期回帰分析

二元配置型周期回帰分析では、二元配置分散分析と同様に全ての測定時点のデータが揃っている被験者だけが解析対象になります。 しかし測定時点が多いとどうしても欠測値が生じやすく、解析対象から除外される被験者が出てくる可能性が高くなります。 また被験者によって測定間隔が異なり、測定時点数が異なっている場合も有り得ます。

そのような時はデータを同一の被験者で連続測定されたものと考えず、測定時点ごとに独立した被験者で測定された一元配置型のデータと考え、一元配置分散分析と周期回帰分析を組み合わせた手法を適用することができます。 その手法を一元配置型周期回帰分析(one-way layout periodic regression analysis)と呼ぶことにします。 表12.1.1のデータを測定時点ごとに独立した被験者から得られた一元配置型のデータと考え、周期成分1と周期成分2を用いた一元配置型周期回帰分析を適用すると次のようになります。 (注2)

表12.3.4 収縮期血圧の一元配置型周期回帰分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
周期回帰18147.444536.8639.3725
ズレ4269.8819224.7311.95029
時期22417.323974.6668.4585
残差2765.524115.229 
全体25182.847 
○周期回帰
周期回帰式:y = 128.9 + 25.4・cos(15t - 210.5) + 10.4・cos(30t - 262.6)
重寄与率:R2 = 0.810(81.0%)
Fβ = 39.373(p = 3.2637×10-10) > F(4,24,0.05) = 2.776 … 有意水準5%で有意
○ズレ
FLOF = 1.950(p = 0.0614) < F(19,24,0.05) = 2.040 … 有意水準5%で有意ではない
○時期:一元配置分散分析における要因A
FP = 8.459(p = 9.457×10-7) > F(23,24,0.05) = 1.993 … 有意水準5%で有意
図12.3.2 平均値と周期回帰曲線

図12.3.2は図12.3.1とほとんど同じですが、各測定時点のデータが独立していることを表すために平均値を●で表し、平均値と平均値を折れ線で結んでありません。 またこの時の周期回帰式と平均値の変動に対する重寄与率は二元配置型周期回帰分析と全く同じです。 そして表12.3.3と表12.3.4を比べると、個人差を表す「被験者」という要因が残差に含まれ、残差の平均平方和が少し大きくなっていることがわかります。 その結果、周期回帰の有意確率とズレの有意確率が少し増え、周期回帰式の信頼性が少し低くなっています。

表12.1.1のデータは例数が2例しかなく、しかも個人差が比較的小さいので二元配置型周期回帰分析と一元配置型周期回帰分析の結果はあまり変わりません。 しかし普通はもっと例数が多く、しかも個人差が大きいので、多少の欠測値があったとしても二元配置型周期回帰分析の方が効率的です。 一元配置型周期回帰分析を適用した方が良い場合は、本当に各時点ごとに独立した被験者から得られたデータか、あるいは測定間隔が被験者ごとに異なるので二元配置型にすると欠測値が非常に多くなってしまうデータを解析する時だけです。

ただし各時点ごとに独立した被検者から得られたデータの時または欠測値が多い時は、平均値の変動が時間によるものか被検者の違いによるものか厳密には区別できません。 そのため結果の解釈は慎重に行う必要があります。 また時系列データには周期関数以外の関数――例えば1次時間数等――を当てはめることもできます。 そこで時系列データの変動を一般的な関数で回帰する手法を時系列回帰分析と呼ぶことにすると、周期回帰分析はその一種ということになります。 (→4.3 繰り返しのある多標本・多時期の計量値 (6) 欠測値の処理方法)


(注1) 表12.1.1を一般化すると次のようになります。

表12.2.5 二元配置型周期回帰分析
の一般的データ
被験者時期平均
t1tjtp
1y11y1jy1pm1.
:::::
iyi1yijyipmi.
:::::
nyn1ynjynpmn.
平均m.1m.jm.pmT

二元配置型周期回帰分析ではデータyijを次のように分解して考えます。

基本式:
:周期回帰式による推定値
周期回帰式:

周期回帰式を行列とベクトルで表現し、最小2乗解を求めると次のようになります。

  β1 = a1 β2 = b1 … β2k-1 = ak β2k = bk … β2m-1 = am β2m = bm
xj0=1
  
 :
  
 :
  

        

        
tan-1(x)関数は-π/2 〜 π/2の値を返すので、次の図のようにしてθkを求める。
図12.2.3 フーリエ係数と位相の関係

二元配置型周期回帰分析の基本式に対応する平方和、自由度、分散、そして分散分析表と検定は次のようになります。

総例数:N = np
全体:  φT = N - 1   
時期:  φP = p - 1   
周期回帰:      
φβ = 2m(2m < pの時)、2m - 1(2m = pの時)   
ズレ:   φLOF = φP - φβ   
被験者:   φSUB = n - 1  
残差:
φR = φT - φP - φSUB = φP×φSUB   
重寄与率:
表12.2.6 二元配置型周期回帰分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
周期回帰SβφβVβFβ=Vβ/VR
ズレSLOFφLOFVLOFFLOF=VLOF/VR
時期SPφPVPFP=VP/VR
被験者SSUBφSUBVSUBFSUB=VSUB/VR
残差SRφRVR 
全体STφT 
○周期回帰の検定
帰無仮説 H0:β1 = … = β2m = 0
Fβ > F(φβR,α)の時、有意水準100α%で有意
○ズレの検定
帰無仮説 H0:周期回帰曲線が時期ごとの平均値の変動と一致している
FLOF > F(φLOFR,α)の時、有意水準100α%で有意
○時期の検定
帰無仮説 H0:m.1 = … = m.p = mT
FP > F(φPR,α)の時、有意水準100α%で有意
○被験者の検定
帰無仮説 H0:m1. = … = mn. = mn.
FSUB > F(φSUBR,α)の時、有意水準100α%で有意

表12.1.1のデータについて、基本周期を24時間とし、周期成分1だけの周期回帰モデルを用いて実際に計算してみましょう。

周期回帰モデル:


     

A0 = 128.938  
  θ1 = π + 0.5322944 = 3.67389(210.498°)
ST = 1072 + … + 1012 - 48×128.93752 = 823177 - 797994.2 = 25182.8   φT = 48 - 1 = 47
SP = 2×103.52 + … + 2×1122 - 48×128.93752 = 820411.53 - 797994.2 = 22417.33   φP = 24 - 1 = 23   
Sβ = 2×(-21.9228)×(-263.074) + 2×(-12.91267)×(-154.952) = 15536.28   φβ = 2  
SLOF = 22417.33 - 15536.3 = 6881.03  φLOF = 23 - 2 = 21   
SSUB = 24×126.54172 + 24×131.33332 - 48×128.93752 = 798269.721 - 797994.2 = 275.521
φSUB = 2 - 1 = 1  VSUB = 275.521
SR = 25182.8-22417.33 - 275.521 = 2489.979   φR = 47 - 23 - 1 = 23   
表12.3.2 収縮期血圧の二元配置型周期回帰分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
周期回帰15536.327768.1471.7545
ズレ6881.0321327.6683.02668
時期22417.323974.6669.00301
被験者275.5211275.5212.54499
残差2489.9823108.26 
全体25182.847 



(注2) 一元配置型周期回帰分析を適用する一般的なデータは次のようなものです。

表12.2.7 一元配置型周期回帰分析
の一般的データ
被験者時期平均
t1tjtp
1y11y1jy1p
::::
iyi1yijyip
::::
njyn1・1ynj・jynp・p
平均m.1m.jm.pmT
※njは時期ごとに異なっていても良い

一元配置型周期回帰分析の基本式、周期回帰式の最小2乗解、平方和、自由度、分散、そして分散分析表と検定は次のようになります。

基本式:
:周期回帰式による推定値
周期回帰式:
  β1 = a1 β2 = b1 … β2k-1 = ak   β2k = bk … β2m-1 = am   β2m = bm
※これら偏回帰係数の最小2乗解を求める方法は(注1)と全く同じため省略
総例数:
全体:  φT = N - 1   
時期:  φP = p - 1   
周期回帰:      
φβ = 2m(2m < pの時)、2m - 1(2m = pの時)   
ズレ:   φLOF = φP - φβ   
残差:   φR = φT - φP   
重寄与率:
表12.2.8 一元配置型周期回帰分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
周期回帰SβφβVβFβ=Vβ/VR
ズレSLOFφLOFVLOFFLOF=VLOF/VR
時期SPφPVPFP=VP/VR
残差SRφRVR 
全体STφT 
○周期回帰の検定
帰無仮説 H0:β1 = … = β2m = 0
Fβ > F(φβR,α)の時、有意水準100α%で有意
○ズレの検定
帰無仮説 H0:周期回帰曲線が時期ごとの平均値の変動と一致している
FLOF > F(φLOFR,α)の時、有意水準100α%で有意
○時期の検定
帰無仮説 H0:m.1 = … = m.p = mT
FP > F(φPR,α)の時、有意水準100α%で有意

表12.1.1のデータを一元配置型のデータとして扱い、基本周期を24時間とし、周期成分1と周期成分2を含んだ周期回帰モデルを用いて実際に計算してみましょう。

周期回帰モデル:
  
  
  
A0 = 128.938     
  θ1 = π + 0.5322944 = 3.67389(210.498°)
  θ2 = π + 1.441365 = 4.582957(262.584°)
ST = 1072 + … + 1012 - 48×128.93752 = 823177 - 797994.2 = 25182.8   φT = 48 - 1 = 47
SP = 2×103.52 + … + 2×1122 - 48×128.93752 = 820411.53 - 797994.2 = 22417.33   φP = 24 - 1 = 23   
Sβ = 2×(-21.9228)×(-263.074) + 2×(-12.91267)×(-154.952) + 2×(-1.34629)×(-16.1554) + 2×(-10.3434)×(-124.12) = 18147.45
φβ = 4  
SLOF = 22417.33 - 18147.45 = 4269.88  φLOF = 23 - 4 = 19   
SR = 25182.8 - 22417.33 = 2765.47  φR = 47 - 23 = 24   
表12.3.4 収縮期血圧の一元配置型周期回帰分析表
要因平方和SS自由度φ平均平方和Ms(分散V)分散比F
周期回帰18147.444536.8639.3725
ズレ4269.8819224.7311.95029
時期22417.323974.6668.4585
残差2765.524115.229 
全体25182.847