12.2 周期回帰分析

(1) 周期解析

表12.1.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)　　A_k：振幅(amplitude) 　　：角周波数(angular frequency)
t：時間　　T：基本周期(fundamental period) 　　：周期成分kの周期(period)
θ_k：位相(acrophase) 　　：位相時間 　　m≦(データ数)/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.1の平均値について基本周期を24時間とし、周期成分1だけの周期回帰式を当てはめると次のようになります。 _(注1)

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

周期回帰の検定は周期回帰式の振幅が0かどうかの検定であると同時に、重寄与率が0かどうかの検定でもあります。 重回帰分析と同様に周期回帰分析も記述統計学的手法のため、推測統計学的手法である検定とは相性が良くありません。 それでも有意症患者のために、回帰誤差が近似的に正規分布するという性質を利用してこのような検定を行います。 しかしこの検定で検出差を指定するのは困難なので、たいていは有意性検定になります。 そのため実質的な意味はほとんどないと考えた方が良いと思います。

そこで周期回帰分析では、検定結果よりも振幅と重寄与率の値そのものを科学的に検討する方が有意義です。 この例題の場合は振幅の値が25.4であり、重寄与率が約69％もあるので、医学的に見て収縮期血圧は周期変動していると言えそうです。

(3) 周期成分の合成

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

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

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

周期成分が2つ以上含まれた周期関数では、最大値と最小値は各周期成分の振幅と位相の組み合わせによって微妙に変化するため、周期成分が1つの周期関数と違って簡単に計算することはできません。 それらの値は周期回帰式に基づいて実際に周期回帰曲線を描き、そのグラフから探索する必要があります。 実際、図12.2.2を見ると、合成周期回帰曲線には極大点と極小点がそれぞれ2つあり、それらの時間は周期成分1の位相時間(約14時間)とも周期成分2の位相時間(約9時間)とも異なることがわかります。 そして最大値になる時間とその時の値をグラフから読み取ると、最大値になるのは11時頃であり、最大値は約150であることがわかります。

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

原理的には周期成分を追加すればするほど周期回帰曲線は元データによりフィットし、12個の周期成分を用いれば周期回帰曲線とデータが完全に一致する、つまり重寄与率が100％になります。 しかしいくらデータとの一致性が高いからといって、科学的に妥当な解釈ができない周期成分を含めるのは意味がありません。 例えば表12.1.1の収縮期血圧については、24時間周期の周期成分1は昼間と夜間の活動を表し、12時間周期の周期成分2は午前中と午後の活動を表すというように医学的に妥当な解釈ができます。 しかしそれ以上周期が短い周期成分については、通常は医学的に妥当な解釈をすることが難しいと考えられます。

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

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

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

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

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

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

また周期回帰分析を適用する前処理として、第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)といいます。 フーリエ級数をk=mまでとし、それ以上を回帰誤差εとした式が周期回帰モデルです。 (→「ベクトルと行列」　10.フーリエ展開)

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

この周期回帰モデルにおいて、次のように置くと重回帰モデルになります。

　　β₁=a₁　β₂=b₁　…　β_2k-1=a_k 　　β_2k=b_k　…　β_2m-1=a_m　β_2m=b_m
x₀=1
　　
　：
　　
　：
　　

　　 　　 　　

最小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.1 周期分散分析表

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

○周期回帰の検定
帰無仮説 H₀：A₁=…=A_k=0
F_β＞F(φ_β,φ_R,α)の時、有意水準100α％で有意

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

ちなみに時系列データに回帰モデルを当てはめた時、回帰誤差に系列相関(自己相関)があるかどうかを検討するための指標としてダービン・ワトソン比(Durbin-Watson ratio)というものがあります。 この値が2に近ければ系列相関はなく、周期変動を回帰モデルによってうまく説明できていると考えることができます。 また定数項だけの回帰モデルを当てはめた時、つまり時期変動がないというモデルのダービン・ワトソン比を求めると、時系列データに系列相関があるかどうかを検討することができます。

ダービン・ワトソン比：　(0≦DW≦4)
ε_j=ρ ε_j-1 + ε_εj の時
h統計量によるρ=0の検定：＞t(∞,α)の時、有意水準100α％で有意
ε_j：回帰誤差　　ε_εj：ε_jの系列誤差 　　ρ：系列相関係数　　V(γ)：y_j=α + β x_j + γ y_j-1 + ε_jにおけるγの分散
DW≒2 → ρ≒0：系列相関なし、ε_jはランダム。
DW≪2 → ρ＞0：正の系列相関あり、ε_jは正が連続で続いた後に負が連続で続くという傾向がある。
DW≫2 → ρ＜0：負の系列相関あり、ε_jは正と負が交互に出現するという傾向がある。

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

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

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

　　 　　

A₀=128.938　　
　　θ₁=π + 0.5322944=3.67389(210.498°)
S_yy=11208.7　　φ_y=23　　S_R=3440.52 　　φ_R=23 - 2=21　　
S_β=11208.7 - 3440.52=7768.14　　φ_β=2 　　
　　

表12.2.2 周期分散分析表

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

上記のように周期回帰の検定結果が有意で、重寄与率が約69％あり、ダービン・ワトソン比が2に近い値なので系列相関はあまりありません。 この結果と図12.2.1から周期回帰曲線によって周期変動をある程度は説明できると考えられます。

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

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

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

　　 　　

A₀=128.938　　 　　
　　θ₁=π+0.5322944=3.67389(210.498°)
　　θ₂=π+1.441365=4.582957(262.584°)
S_yy=11208.7　　φ_y=23　　S_R=2134.94 　　φ_R=23-4=19　　
S_β=11208.7-2134.94=9073.72　　φ_β=4 　　
　　

表12.2.3 周期分散分析表

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

上記のように周期回帰の検定結果が有意で、重寄与率が約81％に増え、ダービン・ワトソン比がさらに2に近づきました。 この結果と図12.2.2から周期回帰曲線の適合度がより高くなったと考えられます。

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

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