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

前口上 | 目次 | 第1章 | 第2章 | 第3章 | 第4章 | 第5章 | 第6章 | 第7章 | 第8章 | 第9章 | 第10章
第11章 | 第12章 | 第13章 | 第14章 | 第15章 | 第16章 | 第17章 | 第18章 | 第19章 | 第20章 | 付録

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9

## 12.6 周期共分散分析

### (1) 周期共分散分析の原理

時系列データに限らず、関連性を持つ一連のデータに特定の関数を当てはめた時、その関数が表すグラフにはレベルを表すパラメーターとパターンを表すバラメーターが存在します。 例えば一連のデータに周期関数を当てはめた周期回帰曲線ではメサーがレベルを表し、振幅と位相がパターンを表します。 また一連のデータに直線y=α+βxを当てはめた回帰直線では定数α(y切片)がレベルを表し、回帰係数β(傾き)がパターンを表します。

回帰直線を当てはめる群が複数ある時、共分散分析を利用して各群の回帰直線が平行かどうか、つまり回帰係数の値が同じかどうかを検討することができます。 そして各群の回帰直線が平行なら各群の目的変数の修正平均値——回帰直線を利用して説明変数の影響を取り除いた平均値——を比較することができます。 その場合、各群の修正平均値の差は各群の平行な回帰直線の定数の差と一致するので、修正平均値の比較は平行な回帰直線の定数の比較に相当します。 したがって共分散分析は複数の回帰直線のレベルとパターンを比較する手法と考えることができます。 (→第8章　共分散分析)

この共分散分析の原理を周期回帰曲線に応用すると、複数の群における周期回帰曲線のレベル(メサー)とパターン(振幅と位相)を比較することができます。 その手法を周期共分散分析(PERCOVA：periodic analysis of covariance)と呼ぶことにします。

### (2) 単純な周期共分散分析

表12.1.1を疾患Aの2名の被験者のデータとし、これに疾患Bの3名の被験者のデータを追加したものが表12.6.1だとします。

表12.6.1 A群とB群の収縮期血圧
群 | 被験者番号 | 0時 | 1時 | 2時 | 3時 | 4時 | 5時 | 6時 | 7時 | 8時 | 9時 | 10時 | 11時
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
A | A1 | 107 | 95 | 93 | 112 | 82 | 114 | 105 | 123 | 135 | 135 | 140 | 137
A | A2 | 100 | 90 | 91 | 122 | 92 | 110 | 106 | 124 | 155 | 165 | 165 | 160
B | B1 | 111 | 104 | 118 | 113 | 136 | 123 | 118 | 128 | 145 | 136 | 150 | 164
B | B2 | 110 | 113 | 111 | 133 | 118 | 142 | 138 | 135 | 143 | 148 | 151 | 136
B | B3 | 123 | 121 | 121 | 118 | 136 | 125 | 131 | 144 | 152 | 149 | 134 | 140
A | 平均値 | 103.5 | 92.5 | 92 | 117 | 87 | 112 | 105.5 | 123.5 | 145 | 150 | 152.5 | 148.5
B | 平均値 | 114.7 | 112.7 | 116.7 | 121.3 | 130 | 130 | 129 | 135.7 | 146.7 | 144.3 | 145 | 146.7
群 | 被験者番号 | 12時 | 13時 | 14時 | 15時 | 16時 | 17時 | 18時 | 19時 | 20時 | 21時 | 22時 | 23時
A | A1 | 147 | 138 | 115 | 161 | 160 | 123 | 142 | 155 | 135 | 131 | 129 | 123
A | A2 | 157 | 148 | 145 | 151 | 160 | 133 | 122 | 145 | 155 | 135 | 120 | 101
B | B1 | 142 | 143 | 116 | 139 | 148 | 134 | 133 | 144 | 141 | 128 | 133 | 103
B | B2 | 133 | 123 | 134 | 125 | 133 | 116 | 134 | 131 | 129 | 112 | 98 | 120
B | B3 | 139 | 132 | 126 | 123 | 111 | 128 | 137 | 143 | 146 | 141 | 119 | 144
A | 平均値 | 152 | 143 | 130 | 156 | 160 | 128 | 132 | 150 | 145 | 133 | 124.5 | 112
B | 平均値 | 138 | 132.7 | 125.3 | 129 | 130.7 | 126 | 134.7 | 139.3 | 138.7 | 127 | 116.7 | 122.3

表12.6.1のA群とB群の平均値について、基本周期を24時間とし、周期成分1と周期成分2だけの周期回帰式を当てはめると次のようになります。

○A群
周期回帰式：y = 128.9 + 25.4・cos(15t - 210.5) + 10.4・cos(30t - 262.6)
メサー：A0 = 128.9　　周期成分1の振幅：A1 = 25.4 　　周期成分2の振幅：A2 = 10.4
周期成分1の位相：θ1 = 210.5　　周期成分1の位相時間 = 14.0時間
周期成分2の位相：θ2 = 262.6　　周期成分2の位相時間 = 8.8時間
重寄与率：R2 = 0.810(81.0％)
○B群
周期回帰式：y = 130.5 + 10.0・cos(15t - 171.0) + 7.8・cos(30t - 240.0)
メサー：A0 = 130.5 　周期成分1の振幅：A1 = 10.0 　　周期成分2の振幅：A2 = 7.8
周期成分1の位相：θ1 = 171.0　　周期成分1の位相時間 = 11.4時間
周期成分2の位相：θ2 = 240.0　　周期成分2の位相時間 = 8.4時間
重寄与率：R2 = 0.820(82.0％)

図12.6.1において、青色の折れ線と青色の曲線がA群の平均値の変動と周期回帰曲線であり、赤色の折れ線と赤色の曲線がB群の平均値の変動と周期回帰曲線です。 このグラフと上記の計算結果から、2群の周期回帰曲線のメサーはあまり変わらないものの、B群の振幅はA群の半分以下であり、位相が3時間ほど早いことがわかります。 つまり2群の周期回帰曲線はレベルはあまり変わらず、パターンが異なっていると解釈できます。

このデータに周期共分散分析を適用すると次のようになります。 (注1)

表12.6.2 周期共分散分析表(PERCOVA table)
要因 | 平方和 | 自由度 | 平均平方和(分散) | F値
--- | --- | --- | --- | ---
群差 | 30.8802 | 1 | 30.8802 | 0.459
共通周期 | 8742.84 | 4 | 2185.71 | 32.481
修正群差 | 30.8802 | 1 | 30.8802 | 0.459
全体周期 | 8742.84 | 4 | 2185.71 | 32.481
非平行性 | 2249.13 | 4 | 562.283 | 8.356
残差 | 2557.09 | 38 | 67.2917 |
全体 | 13579.9 | 47 |

○A群
群別周期回帰式：y = 128.9 + 25.4・cos(15t - 210.5) + 10.4・cos(30t - 262.6)
共通周期回帰式：y = 128.9 + 16.9・cos(15t - 199.7) + 8.9・cos(30t - 253.0)
群別周期回帰式の重寄与率：R2 = 0.810(81.0％)
○B群
群別周期回帰式：y = 130.5 + 10.0・cos(15t - 171.0) + 7.8・cos(30t - 240.0)
共通周期回帰式：y = 130.5 + 16.9・cos(15t - 199.7) + 8.9・cos(30t - 253.0)
群別周期回帰式の重寄与率：R2 = 0.820(82.0％)
○群差：2群の24時間平均値の差の検定
FA = 0.459(p = 0.5022) ＜ F(1,38,0.05) = 4.098 … 有意水準5％で有意ではない
○共通周期：2群の周期回帰曲線が平行と仮定した時の周期回帰の検定
Fβc = 32.481(p = 8.603×10-12) ＞ F(4,38,0.05) = 2.619 … 有意水準5％で有意
○修正群差：2群の周期回帰曲線が平行と仮定した時のメサーの差の検定(レベルの検定)
FAA = 0.459(p = 0.5022) ＜ F(1,38,0.05) = 4.098 … 有意水準5％で有意ではない
○全体周期：2群を合わせて周期回帰式を計算した時の周期回帰の検定
FβT = 32.481(p = 8.603×10-12) ＞ F(4,38,0.05) = 2.619 … 有意水準5％で有意
○非平行性：2群の周期回帰曲線の振幅と位相の差の検定(パターンの検定)
FD = 8.356(p = 6.140×10-5) ＞ F(4,38,0.05) = 2.619 … 有意水準5％で有意

「群差」は2群の24時間平均値の差の検定であり、周期回帰式を当てはめない時の2群のレベルの検定になります。 これは周期共分散分析の主目的ではなく、参考程度です。

「共通周期」は2群の周期回帰曲線が平行と仮定した時の周期回帰の検定、つまり平行な周期回帰曲線の振幅が0かどうかの検定です。 この検定結果が有意ではない時は2群のデータに周期回帰式を当てはめることが妥当ではない、つまり2群の周期回帰曲線のレベルとパターンを比較するのは無意味ということになります。 このデータの場合は有意水準5％で有意ですから、2群の周期回帰曲線のレベルとパターンを比較することが可能です。

「修正群差」は2群の周期回帰曲線が平行と仮定した時のメサーの差の検定つまりレベルの検定です。 原則として、この検定は2群の周期回帰曲線がほぼ平行の時だけ正確な結果になります。 このデータの場合は2群の測定時期が同じため群差の検定結果と同一になります。 これは2群の測定時期が同じ時は2群の24時間平均値の差とメサーの差が同一になるからです。

また測定間隔が等間隔の時は24時間平均値と周期回帰式のメサーが同じ値になります。 このデータの場合は2群の測定時期が同じで、しかも等間隔ですから、群差と修正群差の検定結果が一致し、しかも24時間平均値と周期回帰式のメサーが同じ値になります。 このような時は2群の周期回帰曲線が平行ではなくても検定結果が正確になります。

「全体周期」は2群を合わせて周期回帰式を計算した時の周期回帰の検定、つまり全体の周期回帰曲線の振幅が0かどうかの検定です。 これは周期共分散分析の主目的ではなく、参考程度のものです。 このデータの場合は2群の測定時期が同じため、全体の周期回帰曲線の振幅と位相が平行な周期回帰曲線の振幅と位相に一致し、共通周期の検定結果と一致します。

表12.6.2の周期共分散分析表において、共通周期と修正群差の間に隙間があり、全体周期と非平行性の間にも隙間があります。 これは群差と共通周期の平方和の合計と、修正群差と全体周期の平方和の合計が一致するからです。 このことはデータの変動のうちのある部分については、2群の平均値の違いによる変動と平行な周期回帰曲線による変動を合わせたものと解釈できると同時に、修正平均値つまりメサーの違いによる変動と全体の周期回帰曲線による変動を合わせたものとも解釈できることを意味しています。 これらの事情は普通の共分散分析と同様です。

「非平行性」は2群の周期回帰曲線の振幅と位相の差の検定つまりパターンの検定です。 この検定結果が有意の時は2群の周期回帰曲線が平行ではなく、曲線のパターンが異なっていることになります。 このデータの場合は有意水準5％で有意であり、図12.6.1からわかるように2群の周期回帰曲線のパターンが少し異なっています。

表12.6.1のデータについて、群ごとの平均値ではなく個々の被験者のデータに周期共分散分析を適用することもできます。 しかし個々の被験者は周期的な活動以外に突発的な活動もしていることが多く、データに周期変動だけでなく突発的な変動も含まれている可能性があります。 そのため周期回帰分析と同様に、個々の被験者に周期共分散分析を適用すると往々にして周期回帰曲線がうまくフィットしないことがあります。

個々の被験者のデータを平均すると突発的な変動が小さくなり、群に共通する変動——例えば日内変動のような変動が浮かび上がります。 そのため個々の被験者のデータに周期共分散分析を適用するよりも、平均値に適用する方が周期回帰曲線がうまくフィットする可能性が高く、結果の解釈も妥当なものになる可能性が高くなります。

### (3) 周期回帰曲線の一致係数

普通の共分散分析と同様に、周期共分散分析における検定はほとんど有意性検定になります。 そのため検定結果よりも、周期回帰式そのものや重寄与率を科学的に検討する方が有意義です。 しかし周期回帰曲線の非平行性つまりパターンがどの程度異なっているかについて、周期回帰式から判断するには色々と問題があります。

回帰直線のパターンは回帰係数つまり傾きだけのため、「2群の回帰直線のパターンが異なっている」ということは「2群の回帰直線の傾きが異なっている」という意味になります。 しかし周期回帰曲線のパターンには振幅と位相という2つのパラメーターがあるため、「2群の周期回帰曲線のパターンが異なっている」といっても振幅と位相のどちらが異なっているのかわかりません。 そこで、2群の周期回帰曲線の振幅と位相が一致している程度を表す指標があると便利です。 そしてそれらの指標は関数の相関係数と一致係数の考え方に基づいて導くことができます。

詳しい説明は(注2)を見ていただくとして、関数の相関係数は2つの関数の大小関係が同じかどうかの指標であり、周期回帰曲線の場合は2つの周期回帰曲線の位相が一致しているかどうかを表す位相の一致係数になります。 一方、関数の一致係数は2つの関数の値が同じかどうかの指標であり、周期回帰曲線の場合は2つの周期回帰曲線のメサーが一致するように平行移動した時に、周期回帰曲線が一致するかどうかを表す周期回帰曲線全体の一致係数になります。 そしてそれら2つの一致係数から、2つの周期回帰曲線の振幅が一致しているかどうかを表す振幅の一致係数を導くことができます。

表12.6.1のデータについて、それらの一致係数を計算すると次のようになります。 (注2)

位相の一致係数＝周期回帰曲線の相関係数：rθ = 0.778 　　振幅の一致係数：rA = 0.759
周期回帰曲線全体の一致係数：rP = rθ×rA = 0.591

これらの一致係数から、A群とB群の周期回帰曲線は位相と振幅が同じ程度にずれていて、全体としてはかなりずれていることがわかります。

### (4) 二元配置型周期共分散分析

表12.6.1は繰り返し測定した二元配置型のデータですから、繰り返し測定型二元配置分散分析と周期共分散分析を組み合わせて、より詳細な分析をすることができます。 その手法を二元配置型周期共分散分析(two-way layout periodic analysis of covariance)と呼ぶことにします。 それを表12.6.1のデータに適用すると次のようになります。 (注3)

表12.6.3 二元配置型周期共分散分析表(PERCOVA table)
要因 | 平方和 | 自由度 | 平均平方和(分散) | F値
--- | --- | --- | --- | ---
群差 | 74.1125 | 1 | 74.1125 | 0.384
個体残差 | 578.771 | 3 | 192.924 |
個体 | 652.883 | 4 | 163.221 | 1.437
全体周期 | 18504.3 | 4 | 4626.07 | 40.717
非平行性 | 5397.91 | 4 | 1349.48 | 11.878
ズレ合計 | 5536.32 | 38 | 145.693 | 1.282
残差 | 7839.4 | 69 | 113.614 |
全体 | 37930.8 | 119 |

○群差＝修正群差：2群の24時間平均値の差の検定(レベルの検定)
FA = 0.384(p = 0.5793) ＜ F(1,3,0.05) = 10.128 … 有意水準5％で有意ではない
○個体：被験者ごとの24時間平均値の差の検定
FSUB = 1.437(p = 0.2312) ＜ F(4,69,0.05) = 2.505 … 有意水準5％で有意ではない
○全体周期＝共通周期：2群を合わせて周期回帰式を計算した時の周期回帰の検定
FβT = 40.717(p = 2.220×10-16) ＞ F(4,69,0.05) = 2.505 … 有意水準5％で有意
○非平行性：2群の周期回帰曲線の振幅と位相の差の検定(パターンの検定)
FD = 11.878(p = 2.132×10-7) ＞ F(4,69,0.05) = 2.505 … 有意水準5％で有意
○ズレ合計：2群の時間ごとの平均値と周期回帰曲線とのズレの検定
FLOF = 1.282(p = 0.1832) ＜ F(38,69,0.05) = 1.578 … 有意水準5％で有意ではない
※2群の周期回帰式と重寄与率そして周期回帰曲線の一致係数は、平均値を用いた単純な周期共分散分析の結果と同じ

「群差」は2群の24時間平均値の差の検定であると同時に、2群の周期回帰曲線のメサーの差の検定つまりレベルの検定でもあります。 繰り返し測定した二元配置型のデータは2群の測定時期が必ず同じになるため、2群の24時間平均値の差とメサーの差が同一になります。 そのため群差の検定結果とメサーの差つまり修正群差の検定結果が必ず一致します。

「個体残差」は群差の検定のための誤差であり、群を要因として24時間平均値を一元配置分散分析法で比較する時の残差に相当します。 「個体」は被験者ごとの24時間平均値がばらついているかどうかの検定です。 これは個体差を取り除いて効率的な分析を行うためのものであり、周期共分散分析の主目的ではありません。 「個体残差」と「個体」の間に隙間があるのは、群差と個体残差の平方和を合計したものが個体の平方和になるからです。 そしてこの部分は個体の変動を全変動とし、群を要因Aとした一元配置分散分析に相当します。 これは繰り返し測定型二元配置分散分析と同様です。

「全体周期」は2群を合わせて周期回帰式を計算した時の周期回帰の検定、つまり全体の周期回帰曲線の振幅が0かどうかの検定であると同時に、2群の周期回帰曲線が平行と仮定した時の周期回帰の検定つまり共通周期の検定でもあります。 繰り返し測定した二元配置型のデータは2群の測定時期が必ず同じになるため、全体の周期回帰曲線の振幅と位相が平行な周期回帰曲線の振幅と位相に一致し、共通周期の検定結果と一致します。

「非平行性」は2群の周期回帰曲線の振幅と位相の差の検定つまりパターンの検定です。 この検定結果が有意の時は2群の周期回帰曲線が平行ではなく、パターンが異なっていることになります。

「ズレ合計」は2群の時間ごとの平均値の変動と周期回帰曲線とのズレ(LOF：Lack Of Fit)を合計した値が0かどうかの検定です。 この検定結果が有意の時は時間ごとの平均値の変動と周期回帰曲線のズレが大きいことになります。

表12.6.3より平均値の変動と周期回帰曲線のズレはあまり大きくなく、平均値の変動を周期回帰曲線でうまく近似できることがわかります。 そして平均値だけを用いた単純な周期共分散分析結果と同様に、2群の周期回帰曲線のメサーはあまり変わらず、パターンだけが異なっていることもわかります。

このように二元配置型周期共分散分析ではズレの検討を行うことができるため、できるだけ少ない周期成分で、できるだけ信頼性の高い合理的な周期回帰曲線を用いてレベルとパターンの検討することが可能です。 これがこの手法の特徴であり、平均値を用いた単純な周期共分散分析よりも優れた点です。 そして実際の研究現場で得られる時系列データは繰り返し測定した二元配置型データが多いため、この手法が最も有用だと思います。

### (5) 一元配置型周期共分散分析

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

そのような時は、第3節の分散分析を応用した周期回帰分析と同様に、データを同一の被験者で連続測定されたものと考えず、測定時点ごとに独立した被験者で測定された一元配置型のデータと考え、一元配置分散分析を応用した周期共分散分析を適用することができます。 その手法を一元配置型周期共分散分析(one-way layout periodic analysis of covariance)と呼ぶことにします。

表12.6.1のデータを測定時点ごとに独立した被験者から得られた一元配置型のデータと考え、周期成分1と周期成分2を用いた一元配置型周期共分散分析を適用すると次のようになります。 (注4)

表12.6.4 一元配置型周期共分散分析表(PERCOVA table)
要因 | 平方和 | 自由度 | 平均平方和(分散) | F値
--- | --- | --- | --- | ---
群差 | 74.1125 | 1 | 74.1125 | 0.634
共通周期 | 18504.3 | 4 | 4626.07 | 39.57
修正群差 | 74.1125 | 1 | 74.1125 | 0.634
全体周期 | 18504.3 | 4 | 4626.07 | 39.57
非平行性 | 5397.91 | 4 | 1349.48 | 11.542
ズレ合計 | 5536.32 | 38 | 145.693 | 1.246
時期 | 29512.6 | 47 | 627.928 | 5.371
残差 | 8418.17 | 72 | 116.919 |
全体 | 37930.8 | 119 |

○群差：2群の24時間平均値の差の検定
FA = 0.634(p = 0.4286) ＜ F(1,72,0.05) = 3.974 … 有意水準5％で有意ではない
○共通周期：2群の周期回帰曲線が平行と仮定した時の周期回帰の検定
Fβc = 39.567(p = 2.220×10-16) ＞ F(4,72,0.05) = 2.499 … 有意水準5％で有意
○修正群差：2群の周期回帰曲線が平行と仮定した時のメサーの差の検定(レベルの検定)
FAA = 0.634(p = 0.4286) ＜ F(1,72,0.05) = 3.974 … 有意水準5％で有意ではない
○全体周期：2群を合わせて周期回帰式を計算した時の周期回帰の検定
FβT = 39.567(p = 2.220×10-16) ＞ F(4,72,0.05) = 2.499 … 有意水準5％で有意
○非平行性：2群の周期回帰曲線の振幅と位相の差の検定(パターンの検定)
FD = 11.542(p = 2.704×10-7) ＞ F(4,72,0.05) = 2.499 … 有意水準5％で有意
○ズレ合計：2群の時間ごとの平均値と周期回帰曲線とのズレの検定
FLOF = 1.246(p = 0.2091) ＜ F(38,72,0.05) = 1.571 … 有意水準5％で有意ではない
○時期：平均値の時期変動の検定
FP = 5.371(p = 1.120×10-10) ＞ F(47,72,0.05) = 1.571 … 有意水準5％で有意
※2群の周期回帰式と重寄与率そして周期回帰曲線の一致係数は、平均値を用いた単純な周期共分散分析の結果と同じ

「群差」、「共通周期」、「修正群差」、「全体周期」、「非平行性」の解釈は平均値を用いた単純な周期回帰分析と同様であり、「ズレ合計」の解釈は二元配置型周期回帰分析と同様です。

共通周期と修正群差の間に隙間があるのは、群差と共通周期の平方和の合計と、修正群差と全体周期の平方和の合計が一致するからです。 これは単純な周期共分散分析と同様です。 また全体周期と非平行性の間と、ズレ合計と時期の間にも隙間があります。 これは群差＋共通周期＋非平行性+ズレ合計、または修正群差＋全体周期＋非平行性+ズレ合計の平方和の合計が時期の平方和に一致するからです。 そして「時期」以下は時期を要因Aにした一元配置分散分析に相当します。

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

ただし各時点ごとに独立した被検者から得られたデータの時、または欠測値が多い時は、平均値の変動が時間によるものか被検者の違いによるものか厳密には区別できません。 そのため結果の解釈は慎重に行う必要があります。 (→4.3 繰り返しのある多標本・多時期の計量値　(6) 欠測値の処理方法)

以上のような3種類の周期共分散分析のどれを適用するにしても、次のような周期回帰分析の一般的条件を満たしている必要があります。

1. 周期回帰式に組み込んだ周期成分が科学的に解釈可能か？
2. 周期回帰式のメサー、振幅、位相が科学的に納得のいくものか？
3. 周期回帰式の重寄与率が大きな値か？

周期回帰分析と同様に、この場合も時系列データの変動を一般的な関数――例えば1次時間数等――で回帰し、共分散分析を行うことができます。 そこで時系列データの変動を一般的な関数で回帰して共分散分析を行う手法を時系列共分散分析と呼ぶことにすると、周期共分散分析はその一種ということになります。 (→第8章　共分散分析)

(注1)平均値だけを用いた単純な周期回帰分析は、群ごとの平均値ではなく被験者ごとのデータに対して適用することもできます。 被験者ごとのデータに対して適用する時は、以下の説明の「群」を「被験者」と読み替えてください。

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

表12.6.5 周期共分散分析の
一般的データ
群 | 時期 | 平均
--- | --- | ---
t1 | … | tj | … | t(pi)
1 | y11 | … | y1j | … | y1(p1) | m1.
: | : | … | : | … | : | :
i | yi1 | … | yij | … | yi(pi) | mi.
: | : | … | : | … | : | :
a | ya1 | … | yaj | … | ya(pa) | ma.
全体 | m.1 | … | m.j | … | m.(pi) | mT

周期共分散分析ではデータyijに対して次のような周期回帰式を当てはめて考えます。

○群別周期回帰モデル：群ごとに計算した周期回帰式
k = 1,…,m ≦ p/2　　j = 1,…,pi (群ごとに異なっていても良い)
：群別周期回帰式による推定値 (i = 1,…,a、j = 1,…,pi) 　　A0i：Ai群のメサー、群によって異なる (i = 1,…,a)
Aki：Ai群の振幅、群によって異なる (k = 1,…,m、i = 1,…,a) 　　θki：Ai群の位相、群によって異なる (k = 1,…,m、i = 1,…,a)
○共通周期回帰モデル：各群の周期回帰曲線が平行と仮定して計算した周期回帰式
：共通周期回帰式による推定値 (i = 1,…,a、j = 1,…,pi) 　　Ac0i：Ai群のメサー、群によって異なる (i = 1,…,a)
Ack：全群共通の振幅 (k = 1,…,m) 　　θck：全群共通の位相 (k = 1,…,m)
○全体周期回帰モデル：全群を合わせて計算した周期回帰式
：全体周期回帰式による推定値 (j = 1,…,pi) 　　AT0：全体のメサー
ATk：全体の振幅 (k = 1,…,m)　　θTk：全体の位相 (k = 1,…,m)

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

群別周期回帰：
全体周期回帰：
共通周期回帰：
○全体変動：(yij - mT)
総例数： 平方和： 自由度：φT = N - 1
○群別周期回帰変動：
平方和：
β1i = a1i β2i = b1i…　β(2k-1)i = aki β(2k)i = bki…　β(2m-1)i = amiβ(2m)i = bmi
xj0 = 1
：
：
tan-1(x)関数は-π/2〜π/2の値を返すので、次の図のようにしてθkを求める。
自由度： (2m = piの時は(2m-1)a)
○共通周期回帰変動：
平方和：
自由度：φβc = 2m (2m = pの時は(2m-1))
○全体周期回帰変動：
平方和：
自由度：φβT = 2m (2m = pの時は(2m-1))
○群差変動：(mi. - mT)
平方和： 自由度：φA = a - 1
○修正群差変動：
平方和： 自由度：φAA = a - 1
○非平行性変動：
平方和：
自由度： (2m = pの時は(2m-1)(a-1))
○残差変動：
平方和： 自由度：

これらの変動の関係をグラフで表すと図12.6.2のようになり、平方和間の関係を模式図にすると図12.6.3のようになります。 そして平方和と自由度を表12.6.6のような周期共分散分析表にまとめます。

表12.6.6 周期共分散分析表(PERCOVA table)
要因 | 平方和 | 自由度 | 平均平方和(分散) | F値
--- | --- | --- | --- | ---
群差 | SA | φA | VA=SA/φA | FA=VA/VR
共通周期 | Sβc | φβc | Vβc=Sβc/φβc | Fβc=Vβc/VR
修正群差 | SAA | φAA | VAA=SAA/φAA | FAA=VAA/VR
全体周期 | SβT | φβT | VβT=SβT/φβT | FβT=VβT/VR
非平行性 | SD | φD | VD=SD/φD | FD=VD/VR
残差 | SR | φR | VR=SR/φR |
全体 | ST | φT |

○群差の検定
帰無仮説 H0：m1. = … = ma. = mT
FA ＞ F(φA,φR,α)の時、有意水準100α％で有意
○共通周期の検定
帰無仮説 H0：βc1 = … = βc(2m) = 0
Fβc ＞ F(φβc,φR,α)の時、有意水準100α％で有意
○修正群差の検定
帰無仮説 H0：βc01 = … = βc0a
FAA ＞ F(φAA,φR,α)の時、有意水準100α％で有意
○全体周期の検定
帰無仮説 H0：βT1 = … βT(2m) = 0
FβT ＞ F(φβT,φR,α)の時、有意水準100α％で有意
○非平行性の検定
帰無仮説 H0：β11 = … = β1a = βc1、… 、β(2m)1 = … = β(2m)a = βc(2m)
FD ＞ F(φD,φR,α)の時、有意水準100α％で有意

表12.6.1のA群とB群の平均値について、基本周期を24時間とし、周期成分1と周期成分2を含んだ周期回帰モデルを用いて実際に計算してみましょう。

○A群の群別周期回帰式 → 第2節の(注2)参照
○B群の群別周期回帰式
A0 = 130.542　　A1 = 9.97588　　A2 = 7.76766 　　θ1 = 2.98392(170.966°) 　　θ2 = 4.18886(240.004°)
ST = 821533 - 48×129.742 ≒ 13579.9 　　φT = 48 - 1 = 47
Sβc = (-15.8875)×(-381.299) + (-5.67313)×(-136.155) + (-2.61483)×(-62.7559) + (-8.5352)×(-204.848) ≒ 8742.84
φβc = 4
SβT = (-15.8875)×(-381.299) + (-5.67313)×(-136.155) + (-2.61483)×(-62.7559) + (-8.5352)×(-204.848) ≒ 8742.84
φβT = 4
SA = 24×128.9382 + 24×130.5422 - 48×129.742 ≒ 30.8802 　　φA = 1
SAA = 30.8802 + 8742.84 - 8742.84 = 30.8802 　　φAA = 1
SD = 10991.9 - 8742.84 = 2249.13 　　φD = 8 - 4 = 4
SR = 13579.87 - 30.8802 - 10991.9 = 2557.09 　　φR = 48 - 2×(1+2×2) = 38

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

(注2) 第6章第2節の(注2)で説明したように、2つの変数y1とy2の相関係数をベクトル表示すると次のようになります。

y1の偏差ベクトル： y2の偏差ベクトル：
y1とy2の相関係数：

ここでy1とy2を周期関数とすると、その相関係数は次のようになります。

ｄ1 = y1 - m1 = Ｄxβ1* ｄ2 = y2 - m2 = Ｄxβ2*

測定間隔が等間隔の時は[Ｄx'Ｄx]が対角化し、その対角成分が(n/2)になります。 そのため相関係数は次のように2つのフーリエ係数ベクトルβ1*とβ2*を標準化してお互いに正射影した時の影の射影、つまりβ1*とβ2*のなす角θの余弦になります。

第4節の(注1)で説明したようにフーリエ係数ベクトルを極座標表示した時の偏角が位相になり、動径つまりベクトルの大きさが振幅になります。 そのため2つの周期関数の相関係数は位相が一致している程度を表す位相一致係数と解釈することができます。

また第5節で説明したように振幅を平方した値をパワーと呼び、周期成分ごとにパワーをプロットしたものをパワースペクトルといいます。 また2つの周期関数のフーリエ係数同士を掛け合わせ、それを周期成分ごとにプロットしたものをクロススペクトル(cross spectrum)といいます。

そして2つの周期関数の周期成分kについて、そのクロススペクトルの平方をそれぞれのパワースペクトルの積で割った値のことをコレーレンス(Coherence、関連度)と呼び、それを周期成分の関数としたものをコヒーレンス関数(Coherence function)と呼びます。 これらは2つの周期関数の位相の類似性を表す指標であり、2つの波の干渉しやすさを表す値になります。

周期成分kのコヒーレンス：

フーリエ係数ベクトルの大きさの平方つまりフーリエ係数ベクトル自身の内積はパワースペクトルの合計を表すことになり、2つのフーリエ係数ベクトルの内積はクロススペクトルの合計を表すことになります。 このことから2つの周期関数の相関係数(位相一致係数)を平方した値つまり寄与率は、2つの周期関数のクロススペクトル合計の平方を、それぞれの周期関数のパワースペクトル合計の積で割った値と解釈することができ、周期関数全体のコヒーレンスを表す指標になることがわかります。

一方、第6章第2節の(注2)で説明したように、2つの周期関数y1とy2の一致度を表すエーベルの級内相関係数は次のようになります。

最後の式のrの後ろの項は、2つの周期関数のフーリエ係数ベクトルの大きさの平方の幾何平均と算術平均の比になっています。 この値は2つの周期関数のフーリエ係数ベクトルの大きさが一致している時は1になり、一致していない時は1よりも小さくなり、どちらかのフーリエ係数ベクトルの大きさが0の時は0になります。 したがってこれは2つの周期関数の振幅が一致している程度を表す振幅一致係数と解釈することができます。

そこで周期関数の相関係数つまり位相一致係数をrθで表し、振幅一致係数をrAで表し、エーベルの級内相関係数つまり周期関数全体の一致係数をrpで表すと、これらの指標の間には次のような関係があります。

位相一致係数：
※測定間隔が等間隔ではない時：
振幅一致係数：
※測定間隔が等間隔ではない時：
周期関数全体の一致係数：rp = rθ・rA

ちなみに周期共分散分析において、群数が2つの時の非平行性の平方和は次のようになります。

これは2つのフーリエ係数ベクトルの差ベクトルに関するマハラノビスの平方距離的な値になります。 そして測定間隔が等間隔の時は[Ｓxy1-1+Ｓxy2-1]-1が対角化し、その対角成分が(n/4)になります。 すると非平行性の平方和は2つのフーリエ係数ベクトルの差ベクトルの大きさ||β1*-β2*||を平方して(n/4)をかけた値になります。 このことから群が2つの時の非平行性の検定は本質的に2つのフーリエ係数ベクトルの差ベクトルの大きさが0かどうか、つまり2つのフーリエ係数ベクトルが一致しているかどうかの検定になり、結局のところ周期関数全体の一致係数が0かどうかの検定になることがわかると思います。

(注1)で求めたA群とB群の周期回帰曲線について、3種類の一致係数を実際に計算してみましょう。

β1*'β1* = (-21.9228)2 + (-12.9126)2 + (-1.34629)2 + (-10.3434)2 = 756.1428
β2*'β2* = (-9.85214)2 + 1.566372 + (-3.88337)2 + (-6.72726)2 = 159.8548
β1*'β2* = (-21.9228)×(-9.85214) + (-12.9126)×1.56637 + (-1.34629)×(-3.88337) + (-10.3434)×(-6.72726) = 270.5715
rp = 0.778×0.759 = 0.591

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

表12.6.7 二元配置型周期共分散分析の一般的データ
群 | 被験者 | 時期 | 平均
--- | --- | --- | ---
t1 | … | tj | … | tp
1 | 1 | y111 | … | y1j1 | … | y1p1 | m1.1
: | : | … | : | … | : | :
n1 | y11(n1) | … | y1j(n1) | … | y1p(n1) | m1.(n1)
平均 | m11. | … | m1j. | … | m1p. | m1..
: | : | : | … | : | … | : | :
i | 1 | yi11 | … | yij1 | … | yip1 | mi.1
: | : | … | : | … | : | :
ni | yi1(n1) | … | yij(n1) | … | yip(n1) | mi.(n1)
平均 | mi1. | … | mij. | … | mip. | mi..
: | : | : | … | : | … | : | :
a | 1 | ya11 | … | yaj1 | … | yap1 | ma.1
: | : | … | : | … | : | :
na | ya1(na) | … | yaj(na) | … | yap(na) | ma.(na)
平均 | ma1. | … | maj. | … | map. | ma..
全体 | N | m.1. | … | m.j. | … | m.p. | mT

二元配置型周期共分散分析では、群ごとの時期平均値mij.に3種類の周期回帰モデルを当てはめます。 その式は(注1)の群別周期回帰モデル、共通周期回帰モデル、そして全体周期回帰モデルのyをmに変えたものになるので(注1)を参照してください。

ただし二元配置型周期共分散分析では全ての群の測定時期が同じため、全体周期回帰モデルと共通周期回帰モデルのメサー以外のフーリエ係数が同じになり、全体周期回帰変動と共通周期回帰変動が同一になります。 また各群の平均値の差は共通周期回帰モデルのメサーの差と同じになり、群差と修正群差が同一になります。 そのため二元配置型周期共分散分析ではデータを次のように分解し、その基本式に対する平方和と自由度を求めます。

被験者と時期の二元配置：
群と時期の二元配置：
群別周期回帰：
全体周期回帰：
○全体変動：(yijl - mT)
総例数： 平方和： 自由度：φT = Np - 1
○被験者変動：(mi.l - mT)
平方和： 自由度：φSUB = N - 1
○時期変動：(m.j. - mT)
平方和： 自由度：φP = p - 1
○被験者×時期変動：{yijl - (mi.l + m.j. - mT)}
平方和：
自由度：φS×P = φT - φSUB - φP=φSUB×φP = (N-1)(p-1)
○群差変動：(mi.. - mT)
平方和： 自由度：φA = a - 1
○被験者残差変動：(mi.l - mi..)
平方和： 自由度：φSR = φSUB - φA = N - a
○群−時期変動：(mij. - mT)
平方和： 自由度：φAP = ap - 1
○群×時期変動：{mij. - (mi.. + m.j. - mT)}
平方和：
自由度：φA×P = φAP - φA - φP=φA×φP = (a-1)(p-1)
○残差=被験者残差×時期変動：{yijl - (mi.l + mij. - mi..)}
平方和：
自由度：φR = φSR×P = φT - φSUB - φP-φA×P = φSR×φP = (N-a)(p-1)
○群別周期回帰変動：
平方和：
※フーリエ係数ベクトルと積和ベクトルの内容は(注1)に準ずる
自由度： (2m=pの時は(2m-1)a)
○群別時期変動：(mij. - mi..)
平方和： 自由度：
○群別ズレの変動：
平方和：
自由度： (2m = pの時は0)
○全体周期回帰変動＝共通周期回帰変動：
平方和：
自由度：φβT = φβc = 2m (2m = pの時は(2m-1))
○非平行性変動：
平方和：
自由度： (2m = pの時は(2m-1)(a-1))

表12.6.8 二元配置型周期共分散分析表(PERCOVA table)
要因 | 平方和 | 自由度 | 平均平方和(分散) | F値
--- | --- | --- | --- | ---
群差 | SA | φA | VA=SA/φA | FA=VA/VSR
個体残差 | SSR | φSR | VSR=SSR/φSR |
個体 | SSUB | φSUB | VSUB=SSUB/φSUB | FSUB=VSUB/VSR×P
全体周期 | SβT | φβT | VβT=SβT/φβT | FβT=VβT/VSR×P
非平行性 | SD | φD | VD=SD/φD | FD=VD/VSR×P
ズレ合計 | ∑SLOFi | ∑φLOFi | VLOF=∑SLOFi/∑φLOFi | FLOF=VLOF/VSR×P
残差 | SSR×P | φSR×P | VSR×P=SSR×P/φSR×P |
全体 | ST | φT |

○群差の検定
帰無仮説 H0：m1.. = … = ma.. = mT
FA ＞ F(φA,φSR,α)の時、有意水準100α％で有意
○個体の検定
帰無仮説 H0：m1.1 = … = ma.na = mT
FSUB ＞ F(φSUB,φSR×P,α)の時、有意水準100α％で有意
○全体周期の検定
帰無仮説 H0：βT1 = … βT(2m) = 0
FβT ＞ F(φβT,φSR×P,α)の時、有意水準100α％で有意
○非平行性の検定
帰無仮説 H0：β11 = … = β1a = βc1、… 、β(2m)1 = … = β(2m)a = βc(2m)
FD ＞ F(φD,φSR×P,α)の時、有意水準100α％で有意
○ズレ合計の検定
帰無仮説 H0：
FLOF ＞ F(φLOF,φSR×P,α)の時、有意水準100α％で有意

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

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 　　φR = φSR×P = (5-2)×(24-1) = 69
※群別周期回帰式は(注1)と同じ
SβT = 5×3700.856 = 18504.28　　φβT = 2×2 = 4
SD = 23902.19 - 18504.28 = 5397.91 　　φD = 2×2×(2-1) = 4

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

(注4)一元配置型周期共分散分析の一般的データは、表12.6.7において各群の被験者が測定時期ごとに独立したものになります。 そのため群ごとに時期t1〜tpが必ずしも同一ではなく、被験者数niも測定時期ごとに必ずしも同一ではなくなります。 したがって単純な周期共分散分析と同様に、全体周期回帰変動と共通周期回帰変動が同一になるとは限らず、群差と修正群差も同一になるとは限りません。

しかし一元配置型周期共分散分析でも群ごとの時期平均値mij.に群別周期回帰モデル、共通周期回帰モデル、そして全体周期回帰モデルを当てはめます。 そしてその式は(注1)の群別周期回帰モデル、共通周期回帰モデル、そして全体周期回帰モデルのyをmに変えたものになります。

この手法ではデータを次のように分解し、その基本式に対する平方和と自由度を求めます。

時期の一元配置：(yijl - mT) = (mij. - mT) + (yijl - mij.)
群別周期回帰：
全体周期回帰：
共通周期回帰：
○全体変動：(yijl - mT)
総例数： 平方和：
自由度：φT = N - 1
○時期変動：(mij. - mT)
平方和：
自由度：φP = P - 1
○残差変動：(yijl - mij.)
平方和： 自由度：φR = φT - φP = N - P
○群差変動：(mi.. - mT)
平方和： 自由度：φA = a - 1
○群別周期回帰変動：
平方和：
dxjki = nij(xjk-mx.ki)
自由度：
○群別ズレの変動：
平方和：
自由度：
○全体周期回帰変動：
平方和：
自由度：φβT = 2m (2m = pの時は(2m-1))
○共通周期回帰変動：
平方和：
自由度：φβc = 2m (2m = pの時は(2m-1))
○修正群差変動：
平方和： 自由度：φAA = a - 1
○非平行性変動：
平方和： 自由度：

表12.6.9 一元配置型周期共分散分析表(PERCOVA table)
要因 | 平方和 | 自由度 | 平均平方和(分散) | F値
--- | --- | --- | --- | ---
群差 | SA | φA | VA=SA/φA | FA=VA/VR
共通周期 | Sβc | φβc | Vβc=Sβc/φβc | Fβc=Vβc/VR
修正群差 | SAA | φAA | VAA=SAA/φAA | FAA=VAA/VR
全体周期 | SβT | φβT | VβT=SβT/φβT | FβT=VβT/VR
非平行性 | SD | φD | VD=SD/φD | FD=VD/VR
ズレ合計 | ∑SLOFi | ∑φLOFi | VLOF=∑SLOFi/∑φLOFi | FLOF=VLOF/VR
時期 | SP | φP | VP=SP/φP | FP=VP/VR
残差 | SR | φR | VR=SR/φR |
全体 | ST | φT |

○群差の検定
帰無仮説 H0：m1.. = … = ma.. = mT
FA ＞ F(φA,φR,α)の時、有意水準100α％で有意
○共通周期の検定
帰無仮説 H0：βc1 = … = βc(2m) = 0
Fβc ＞ F(φβc,φR,α)の時、有意水準100α％で有意
○修正群差の検定
帰無仮説 H0：βc01 = … = βc0a
FAA ＞ F(φAA,φR,α)の時、有意水準100α％で有意
○全体周期の検定
帰無仮説 H0：βT1 = … βT(2m) = 0
FβT ＞ F(φβT,φR,α)の時、有意水準100α％で有意
○非平行性の検定
帰無仮説 H0：β11 = … = β1a = βc1、… 、β(2m)1 = … = β(2m)a = βc(2m)
FD ＞ F(φD,φR,α)の時、有意水準100α％で有意
○ズレ合計の検定
帰無仮説 H0：
FLOF ＞ F(φLOF,φR,α)の時、有意水準100α％で有意
○時期の検定
帰無仮説 H0：m11. = … = map.
FP ＞ F(φP,φR,α)の時、有意水準100α％で有意

表12.6.1のデータについて、基本周期を24時間とし、周期成分1と周期成分2を含んだ周期回帰モデルを用いて実際に計算すると残差以外は(注3)の計算結果と同じ値になります。 そして残差の平方和SRと自由度φRは(注3)の個体残差の平方和SSRと自由度φSR、残差の平方和SSR×Pと自由度φSR×Pをそれぞれ合わせた値になります。 それらの統計量を用いて表12.6.4の周期共分散分析表を作成することができます。