12.2 周期回帰分析

(1) 周期解析

表12.1の収縮期血圧には夜間が低く、昼間は高いという変動があります。 この変動は人間の活動状態を反映したもので、ほぼ24時間単位で変動する周期変動と考えられます。 このような周期変動のことを「日内変動」または「概日リズム(circadian rhythm)」といいます。 人間をはじめとする多くの生物が日内変動のような周期変動をするため、それを研究する学問として「時間生物学(chronobiology)」という分野ができています。

日内変動のような周期変動を、三角関数――懐かしの正弦関数(sine)または余弦関数(cosine)です――に基づいた周期関数で近似して分析する手法のことを「周期解析(periodic analysis)」といいます。 この手法は元々は光や音といった物理現象を分析する目的で開発されたため、時間生物学分野だけでなく物理学分野や工学分野でも盛んに利用されています。

周期解析には様々な手法がありますが、回帰分析の原理を利用して、測定されたデータに最もフィットする周期関数を求める「周期回帰分析(periodic regression analysis)」と、求めた周期関数にはどのような周期成分(periodic component)が含まれているかを分析する「スペクトル解析(spectrum analysis)」が代表的なものです。

物理学分野や工学分野で研究対象にする光や音は振動であり、その周期変動は原理的に三角関数になります。 そのため周期成分の分析や周期変動原理の検討といったことが主目的になり、スペクトル解析が多用されます。 スペクトル解析は光をプリズムで色々な色に分解することに相当する手法なので、光や音の分析に向いているのです。

それに対して人間の血圧の周期変動が原理的に三角関数になるとは考えにくく、むしろ人間の活動状態を反映して生じた擬似周期変動であると考えられます。 そこで時間生物学分野で周期解析を利用する時は、(今のところ)内容不明な周期変動を三角関数によって近似し、その変動パターンを数量的に検討したり、変動をシミュレートしたりすることが主目的になり、周期回帰分析が多用されます。

これは、内容不明な(ことが多い)薬物の用量−反応関係を直線で近似し、それによって薬理作用の変化を数量的に検討したり、用量−反応をシミュレートしたりすることと似ています。

したがって生物の周期変動を三角関数よりも的確に近似できる別の関数があれば、それを利用する方が良いでしょう。 そしてもしも周期変動関数を科学的な理論から決定できるのであれば、その関数を利用するのが理想的です。

(2) 周期回帰分析

周期回帰分析はフーリエ展開(Fourier expansion)という数学的手法を利用するため、「フーリエ解析(Fourier analysis)」または「調和解析(harmonic analysis)」とも呼ばれます。 この手法は時系列解析の一種であると同時に回帰分析の一種でもあり、最小2乗法を利用してデータの時間的な変動を直線の代わりに三角関数で近似します。 このため普通の回帰分析で行うこと、例えば周期回帰式を求めたり、寄与率を求めたりすることは全て同じように行えます。

また周期回帰分析は説明変数が時間で、目的変数が時間変動する測定値なので、説明変数が原因で目的変数が結果という因果関係がはっきりしています。 (→5.1 相関係数と回帰直線　(2) 回帰分析)

通常、周期回帰分析では次のような周期回帰モデルを仮定します。

A₀は定数であり、周期解析ではメサー(MESOR、Midline Estimating Statistic Of Rhythm)と呼ばれることがあります。 定数は直線回帰式ではy切片になりますが、周期回帰分析では1基本周期あたりの周期関数の平均値になり、データの測定間隔が等間隔の時は1基本周期あたりのデータの平均値と一致します。

振幅は各周期成分がメサーの上下に周期変動する時の最大変動幅であり、メサーから各周期成分の最大値までの距離になります。 このため、各周期成分の最大変動幅は振幅の2倍になります。

角周波数は各周期成分の変動速度を表す値で、基本周期Tを1回転つまり360度で表した時に、各周期成分が単位時間あたりに回転する角度を表します。 例えば基本周期T=24時間とすると、周期成分1の角周波数は(360°/24hr)×1=15°/hrになり、1時間あたり15度だけ回転することに相当します。

位相は各周期成分が開始時から最大値になるまでの時間を角度で表した値であり、この値を角周波数で割ったものが位相時間、つまり開始時から最大値になるまでの時間になります。

振幅と位相は回帰係数に相当する値であり、これらの値とメサーを最小2乗法を利用して求めます。 したがって周期回帰分析は、単回帰分析ではなく重回帰分析に相当します。 重回帰分析には、未知パラメーターつまり回帰係数と定数の合計数がデータ数以下でなければ計算できないという制限があります。 周期回帰分析でも同様の制限があり、各周期成分ごとに2つの未知パラメーターつまり振幅と位相があるため、周期成分数mは(データ数)/2以下でなければなりません。

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

上記の周期回帰式をグラフ化したものが、図12.2の赤い曲線で表した周期回帰曲線です。 この周期回帰曲線は位相時間つまり14時に最大値になり、その値はメサー＋振幅＝128.9+25.4=154.3になります。

周期回帰の検定は、周期回帰式の振幅が0かどうかの検定です。 この検定もこれまでの検定と同様に、たとえ有意になっても、

「データから求められた振幅の値は信頼できるから、母集団の振幅はほぼ確実に0ではない」

つまりは、

「データから得られた振幅の値が信頼できる」

ということを意味しているだけです。

この例題の場合は周期回帰の検定結果が有意ですから、振幅の値が信頼できます。 そしてその振幅の値は25.4であり、重寄与率が約69％もあるので、医学的に見て収縮期血圧は周期変動があると言えそうです。

(3) 周期成分の合成

図12.2を見ると、周期回帰曲線が最大になる14時頃に実際のデータが少し低下しているため、この部分では周期回帰曲線とデータのズレが比較的大きいことがわかります。 この14時頃の低下は、昼食時に休憩を取り、それによって血圧が低下することを反映していると考えられます。 このことから、人間の活動には昼間と夜間という24時間周期のものと、午前中と午後という12時間周期のものがあり、それが血圧に反映されていると考えられます。

そこで今度は24時間周期の周期成分1と、12時間周期の周期成分2を組み合わせた周期回帰式を当てはめてみましょう。 _(注2)

上記の周期回帰式をグラフ化したものが、図12.3の赤い曲線で表した合成周期回帰曲線です。 この合成周期回帰曲線は、青い破線で表した24時間周期の周期成分1と、同じく青い破線で表した12時間周期の周期成分2を合成したものです。 12時間周期の周期成分2を追加することによって14時頃に血圧が少し低下し、図12.2の周期回帰曲線よりも元データにうまくフィットしていることがわかります。

周期成分が2つ以上含まれた周期関数では、最大値と最小値は各周期成分の振幅と位相の組み合わせによって微妙に変化するため、周期成分が1つの周期関数と違って簡単に計算することができません。 それらの値は周期回帰式に基づいて実際に周期回帰曲線を描き、そのグラフから探索する必要があります。

実際、図12.3を見ると、合成周期回帰曲線には極大点と極小点がそれぞれ2つあり、それらの時間は周期成分1の位相時間(約14時間)とも周期成分2の位相時間(約9時間)とも異なることがわかります。 そして最大値になる時間とその時の値をグラフから読み取ると、最大値になるのは11時頃であり、最大値は約150であることがわかります。

(4) 周期回帰分析の注意点

原理的には、周期成分をもっと追加すれば周期回帰曲線はデータにもっとフィットし、12個の周期成分を用いれば周期回帰曲線とデータが完全に一致する、つまり重寄与率が100％になります。 しかしいくらデータとの一致性が高いからと言って、科学的に妥当な解釈ができない周期成分を含めるのは意味がありません。

例えば表12.1の収縮期血圧については、24時間周期の周期成分1は昼間と夜間の活動を表し、12時間周期の周期成分2は午前中と午後の活動を表すというように、医学的に妥当な解釈ができます。 しかしそれ以上周期が短い周期成分については、通常は医学的に妥当な解釈をすることが難しいと考えられます。

周期回帰式に限らず回帰分析では、モデルにする関数を複雑にすればするほどデータとの一致性は高くなります。 しかし科学的に妥当な解釈ができないモデル関数では意味がありませんし、たとえ解釈できるとしても、あまりに複雑なモデル関数は結果の解釈が難しくなります。

そのようなことから、普通の回帰分析では最も単純な直線モデルつまり一次関数を利用し、周期回帰分析では周期成分1だけの周期関数を利用するのが一般的です。 このため日内変動の場合は24時間周期の周期成分1だけを利用するのが一般的であり、表12.1のようにデータが多数の時点で測定されている場合でも、せいぜい12時間周期の周期成分2を追加する程度です。

以上のことから、周期回帰分析結果を評価する時の注意点をまとめると次のようになります。

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

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

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

もちろん周期回帰分析を適用する前処理として、前節で説明した移動平均法を利用してデータの変動を平滑化し、周期回帰曲線をよりフィットしやすくするという方法もあります。 しかし個々の被験者のデータを平均することによって突発的な変動が小さくなるため、平均値の変動はある程度平滑化されています。 そのため平均値に周期回帰分析を適用する場合は、前処理として移動平均法を利用することは少ないようです。

(注1)　n個の時系列データ{y₁,…,y_j,…,y_n}は連続的に変化する時系列データを離散的に測定したものであり、測定間隔Δt→0、時期数n→∞の時、連続関数y(t)になると考えられます。 一般に、任意の連続関数y(t)は次のような無限級数によって表現できます。

：メサー　　：振幅
T：基本周期　　：周期成分kの周期
：角周波数　　：位相、

これを「フーリエ級数(Fourier series)」といい、y(t)をフーリエ級数で表現することを「フーリエ展開(Fourier expansion)」、係数a_k、b_kを「フーリエ係数(Fourier coefficient)」といいます。 (→「ベクトルと行列」　10.フーリエ展開)

フーリエ級数をk=mまでとし、それ以上を誤差εとした式が周期回帰モデルになります。

k=1,…,m≦n/2　　t=t₁,…,t_j,…,t_n

この周期回帰モデルにおいて、

……


　　：

　　：

と置くと、次のように重回帰モデルに相当します。

重回帰分析と同様に、最小2乗法による最小2乗解と各種のパラメーターは次のようになります。 (→7.1 重回帰モデル　(注1)、7.3 変数の選択　(注1))

tan^-1(x)関数は-π/2〜π/2の値を返すため、次の図のようにしてθ_kを求める。

全体平方和：
全体自由度：φ_y=n-1
周期回帰平方和：
　
　
周期回帰自由度：φ_β=2m(2m＜nの時)、2m-1(2m=nの時)
周期回帰分散：
残差平方和：
残差自由度：φ_R=φ_y-φ_β
残差分散：
重寄与率：

表12.2 周期分散分析表

要因	平方和SS	自由度φ	平均平方和Ms(分散V)	分散比F
回帰	Sβ	φβ	Vβ	Fβ=Vβ/VR
残差	SR	φR	VR
全体	Syy	φy

・周期回帰の検定

帰無仮説H₀：A₁=…=A_k=0
F_β≧F(φ_β,φ_R,α)の時、有意水準100・α％で有意

周期回帰の検定は、周期回帰誤差εが近似的に正規分布することを必要とします。 これは目的変数そのものの正規性ではなく、周期回帰曲線と実際のデータの誤差の正規性であることに注意してください。 また説明変数xは確率変数ではなく研究者が任意の値を指定することができる変数のため、誤差がなく、正規性を必要としないことにも注意する必要があります。

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

周期回帰モデル：
　　t=0,1,2,…,23

このデータのように測定間隔が等間隔の場合、cos(x)とsin(x)が直行するため[Ｘ'Ｘ]が対角行列になり、計算が簡単になります。

表12.3 周期分散分析表

要因	平方和SS	自由度φ	平均平方和Ms(分散V)	分散比F
周期回帰	7768.14	2	3884.07	23.707
残差	3440.52	21	163.834
全体	11208.7	23

…有意水準5％で有意

(注2)　表12.1の平均値について、周期成分2を追加した周期回帰モデルを用いて実際に計算してみましょう。

周期回帰モデル：
　　　　t=0,…,23

測定間隔が等間隔の場合、[Ｘ'Ｘ]が対角行列になるため、周期成分1と周期成分2を独立に計算できるようになります。 このため、メサーと周期成分1の振幅と位相は変わりません。

表12.4 周期分散分析表

要因	平方和SS	自由度φ	平均平方和Ms(分散V)	分散比F
周期回帰	9073.72	4	2268.43	20.188
残差	2134.94	19	112.365
全体	11208.7	23

…有意水準5％で有意

最終更新日：2010年6月1日 　第1節へ 第3節へ

　webmaster@snap-tck.com 　Copyleft (C) 2000-2012 SNAP(Sugimoto Norio Art Production)