12.5 スペクトル解析

(1) パワースペクトル

第2節の(3)で24時間周期成分と12時間周期成分を合成した時、24時間周期成分の振幅と位相の値は、この周期成分だけで周期回帰式を求めた時の値と変わりませんでした。 周期回帰分析は重回帰分析に相当し、振幅と位相は偏回帰係数に相当するので、複数の周期成分を合成すればお互いの影響で振幅と位相の値は少し変化します。 しかしこの例のように時系列データが等間隔で測定されていると各周期成分は独立になり、振幅と位相の値がお互いに影響を受けなくなります。

各周期成分が独立ということは、それらを独立に検討することができるということであり、各周期成分ごとの重寄与率を求めることができます。 例えば第2節の(3)で求めた周期回帰式について、24時間周期成分と12時間周期成分の重寄与率を求めると次のようになります。 _(注1)

上の結果から24時間周期成分の重寄与率は約69％とかなり高いことと、12時間周期成分の重寄与率は約12％と比較的低いこと、そしてそれらを合計した値が全体の重寄与率約81％になることがわかります。 また各周期成分ごとの検定は各周期成分の振幅が0かどうかの検定であり、重回帰分析における偏回帰係数の検定と同じようなものです。

第2節で説明したように合成できる周期成分の最大数は(データ数)/2であり、このデータの場合は24/2 = 12個になります。 そこで上記の周期成分以外の10個の周期成分について振幅と位相、そして重寄与率を計算すると次のようになります。 この結果から24時間周期と12時間周期以外の周期成分はどれも重寄与率が小さく、周期変動にほとんど寄与していないことがわかります。 また12個の重寄与率を合計すると100％になり、元のデータと周期回帰曲線がぴったりと一致します。

これらの重寄与率は各周期成分の平方和から求めるので振幅を平方した値と比例します。 そこで周期解析では振幅を平方した値のことをパワー(power)と呼び、各周期成分が周期変動に対してどの程度寄与しているか、つまり各周期成分がどの程度含まれているかを表す目安にします。 そして周期成分の番号を横軸にし、パワーを縦軸にプロットしたグラフをパワースペクトル(power spectrum)と呼びます。 上の結果をパワースペクトルで表すと図12.5.1のようになります。 ただし1つの周期成分はsin関数とcos関数を合成したものなので、このグラフでは振幅を平方して2分の1にした値をパワーとしてプロットしています。

名前からわかるように、周期解析におけるパワースペクトルは光をプリズムで色々な波長の色に分けた時のスペクトルに相当します。 そして周期成分の番号は1基本周期あたり三角関数の波がいくつあるかを表す値になるので周波数(frequency)と呼ばれます。 ただし一般には周波数を単位時間当たりの波数で表すので、この場合は周期成分番号を24で割って1時間当たりの値にするのが普通です。

一方、電気や電波を扱う工学分野では1秒当たりの周波数を用いるのが普通であり、単位としてヘルツ(Hz)を用います。 そして周期が長い周期成分のことを低周波(LF：Low Frequency)、周期が短い周期成分のことを高周波(HF：High Frequency)といいます。 それに対して力学的な運動を扱う物理学分野では、周波数の代わりに振動数(英語ではどちらもfrequency)という用語を用います。 そしてご存知のように音波の場合は周期が長い周期成分のことを低音、周期が短い周期成分のことを高音といいます。 ちなみに振幅を平方した値をパワーと呼ぶ理由は、振幅を平方した値と振動エネルギーが比例するからです。

(2) スペクトル解析

図12.5.1のパワースペクトルは、スペクトルが線で表されるので線スペクトルと呼ばれます。 しかし元のデータの測定間隔を無限小にしてデータ数が無限大になると、理論的には周期成分の数も無限大になり、無限個の線スペクトルがつながって滑らかに変化する連続スペクトルになると考えられます。 その時、元のデータを時間tと収縮期血圧yの関係を表す周期関数y = f(t)で表すと、連続スペクトルは周期成分の番号つまり周波数fと振幅の平方つまりパワーpとの関係を表すスペクトル関数p = f(f)になります。

このように時間領域における周期関数を、周波数領域におけるスペクトル関数に変換することをフーリエ変換(Fourier transform)といいます。 そしてフーリエ変換を利用して、周期関数にどのような周期成分がどの程度含まれているか分析する手法をスペクトル解析(spectrum analysis)といいます。 _(注2)

スペクトル解析は主として工学分野で用いられ、時間生物学分野ではあまり用いられません。 時間生物学分野で周期解析を利用する時は、内容不明な周期変動を三角関数によって近似し、その変動パターンを数量的に検討したり、変動を予測したりすることが主目的になるからです。 例えば前述の例でも12個の周期成分の重寄与率または図12.5.1のパワースペクトルを見て、どのような周期成分がどの程度含まれているかを検討することよりも、重寄与率まはたパワーが大きい周期成分の中で医学的に解釈可能なものだけを選択し、収縮期血圧の日内変動の様子を検討する方が多いと思います。

また医学分野では「何でも検定しなければ気がすまない！」という有意症患者が多いので、スペクトルよりも各周期成分の振幅が0かどうかという検定の方を重視しがちです。 しかし全ての周期成分を合成すると、元のデータと周期回帰曲線が一致してしまいズレがなくなります。 そのため平均値だけを利用する単純な周期回帰分析では、各周期成分の検定も全体の周期回帰式の検定もできません。

そこでこのような時は二元配置型周期回帰分析を用いて各周期成分の検定と全体の周期回帰式の検定を行うことが考えられます。 前述の例について二元配置型周期回帰分析を適用すると次のようになります。 _(注3)

この周期回帰分析表から、有意水準5％で有意な周期成分——F値が2.895以上の周期成分——は周期成分1(24時間周期)、周期成分2(12時間周期)、周期成分4(6時間周期)、周期成分6(4時間周期)、周期成分7(3.429時間周期)であり、12個の周期成分を合成した周期回帰式の検定結果は時期変動の検定結果と一致することがわかります。

ただしこの手法における検定も有意性検定なので実質的な意味はほとんどなく、周期と振幅と重寄与率を医学的に検討する方が有意義です。 上記の結果では周期と振幅と重寄与率の値から、収縮期血圧の日内変動は24時間周期の変動に12時間周期の変動が少し加わったものと解釈できることがわかります。 6時間周期や4時間周期や3.429時間周期の成分は、検定結果は有意ですが、振幅と重寄与率が小さくて医学的な解釈も難しいのでとりあえず無視してかまわないでしょう。

(注1)　フーリエ級数をk = mとした周期回帰モデルにおいて、測定間隔が等間隔の時は[Ｘ'Ｘ]が対角行列になるので、各周期成分を独立に計算できるようになると同時に計算が簡単になります。 さらにデータ数が偶数の時は次のように計算がより簡単になります。 (→12.2 周期回帰分析 (注1))

周期回帰モデル：
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
　　
　：
　　
　：
　　

　　　　

　　

　　 　　 　　

全体平方和： 　　全体自由度：φ_y = n - 1
残差平方和： 　　残差自由度：φ_R = φ_y - φ_β 　　残差分散：
○周期成分k
平方和： (ただし2m = nかつk = mの時はS_β(2k) = 0)
自由度：φ_k = 2(2m ＜ nまたはk ＜ mの時)、1(2m = nかつk = mの時)
分散： 　　重寄与率：
検定の帰無仮説 H₀：A_k = 0
F_k = V_k/V_R ＞ F(φ_k,φ_R,α)の時、有意水準100α％で有意
○全周期
平方和： 　　自由度： 　　分散： 　　重寄与率：
検定の帰無仮説 H₀：A₁ = … = A_m = 0
F_β = V_β/V_R ＞ F(φ_β,φ_R,α)の時、有意水準100α％で有意

これ以外にも、sin関数とcos関数の特性を利用して計算をさらに簡略化することができます。 そのような計算法を高速フーリエ変換(FFT：Fast Fourier Transform)といいます。

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

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

A₀ = m_y = 128.938 　　 　　
S_yy = 11208.66　　φ_y = 23 　　S_R = 11208.66 - 9073.72 = 2134.94　　φ_R = 23 - 4 = 19 　　
S₁ = 12×25.4430² = 7768.15　　φ₁ = 2 　　 　　

S₂ = 12×10.4306² = 1305.57　　φ₂ = 2 　　 　　

S_β = 7768.15 + 1305.57 = 9073.72　　φ_β = 2 + 2 = 4 　　 　　

(注2)　フーリエ変換は数学的にはベクトルの直交分解に相当し、スペクトル解析だけでなく色々な分野で利用されています。 例えば量子力学の波動方程式や、画像のJPEG圧縮法などもフーリエ変換を利用しています。 (→「ベクトルと行列」　11.フーリエ変換)

無限個のデータをフーリエ変換すると連続スペクトルが得られますが、実際に得られるデータは有限個です。 そこでスペクトル解析では、有限個のデータから近似的な連続スペクトルを得る方法が色々と開発されています。 しかしスペクトル解析は医学・薬学分野ではあまり利用されないので詳しい説明は省略します。

ちなみに周期解析では時間に対するデータの変動をプロットしたグラフをクロノグラム(chronogram)といい、周期に対する振幅やパワーの変動をプロットしたグラフをペリオドグラム(periodgram)といいます。 したがって周期回帰曲線はクロノグラムの一種であり、パワースペクトルはペリオドグラムの一種です。

(注3)　分散分析を応用した周期回帰分析も、測定間隔が等間隔の時は各周期成分を独立に計算できるようになると同時に計算が簡単になります。 しかし各周期成分を独立に計算するためには全ての被験者の測定時点が同一で、しかも全ての測定時点のデータが揃っている、つまり二元配置型でなければなりません。 そのため各周期成分を独立に計算することができるのは二元配置型周期回帰分析だけになります。 (→12.3 分散分析を応用した周期回帰分析 (注1))

表12.1.1のデータについて、基本周期を24時間とし、周期成分1〜12を含んだ周期回帰モデルを用いた二元配置型周期回帰分析を実際に適用してみましょう。

周期回帰モデル：

A₀ = m_y = 128.938 　　…　
S₁ = 2×12×25.4430² = 15536.3　　φ₁ = 2
　：
S₁₂ = 2×24×1.52083² = 111.021　　φ₁₂ = 1
S_β = 15536.3 + … + 111.021 = 22417.3　　φ_β = 23
S_T = 25182.8　　φ_T = 47　　S_P = 22417.33 　　φ_P = 23　　S_LOF = 0　　φ_LOF = 0
S_SUB = 275.521　　φ_SUB = 1　　S_R = 2489.979 　　φ_R = 23

これらの値を用いて表12.5.1の二元配置周期回帰分析表を作成することができます。

最終更新日：2016年3月6日 　第4節へ 第6節へ

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