Excel (Microsoft) を用いた最小自乗余弦曲線法の 開発
著者 市丸 雄平, 小林 美佳子, 沈 恵芳, 安原 朋美
雑誌名 東京家政大学研究紀要 2 自然科学
巻 42
ページ 11‑17
発行年 2002
出版者 東京家政大学
URL http://id.nii.ac.jp/1653/00010709/
Exce1(Microsoft)を用いた最小自乗余弦曲線法の開発
市丸 雄平, 小林 美佳子,沈 恵芳,
(平成13年10月4日受理)
安原 朋美
Development of Least−Square Cosine Spectrum Analysis System using Excel(Microsoft)as a Method for the Biologic Time Series Analysis in the Field of Nutrition.
Yuhei lcHiMARu, Mikako KoBAYAsHI, Hui Fang SHEN and Tomomi YAsuHARA
(Received on October 4,2001)
キーワード 最小自乗余弦曲線,時間栄養学,胃電図,エネルギー摂取リズム
Key words Least square analysis, cosine curve, biological rhythm electrogastrogram energy intake
はじめに
生体リズムは,大きくウルトラディアンリズム,サー カディアンリズム,インフラディアンリズムに分類され ている.このリズムを推定する方法としてフーリエ解析 法,最大エントロピー法,最小自乗余弦法などがある.
データ取得において,フーリエ解析あるいは最大エント ロピーはサンプル間隔が等間隔であることが前提である.
また,サンプリング数は2の整数乗であり,一般に連続 関数を前提としている,しかし,生体現象は非等間隔サ ンプリングが困難なことが多いこと,またかならずしも 連続でない点を考慮すると,限られた時間幅のリズムの 振幅および周期を正確に求める方法としては最小自乗余 弦法がすぐれている1)〜4).しかし,最小自乗余弦法は,
周期成分が既定値となっているたあ,われわれは,これ を発展させ,UNIX下でC言語を用いて,最小自乗余弦 スペクトル法を開発し,周波数の推定も行った5).しか し,この使用に際して,習熟を必要とした.今日,より 簡易な表計算のアプリケーションとして,Exce1(Micr osoft社)が広汎に使用されるようになってきた. Excel はデータを対話的に取り扱う方法として優れ,とくにア プリケーション・プログラムあるいは数学・統計関数を 用いると,多くの数値演算処理が可能となる.しかし,
現時点では,最小自乗余弦法あるいは,最小自乗余弦ス
ペクトル法による数値演算法は搭載されていない.今回,
Exce1を用いて最小自乗余弦曲線法を開発し,時間栄養 学の解析手段として応用を試みたので報告する.
臨床栄養第2研究室
方 法
1)Exce1のVisual Basicによる最小自乗余弦曲線法の 数学的背景(図1)。
基本的に,推定する余弦曲線を Y=M十Acos(ωt一θ)十M
とし,最小自乗法をもちいた.ここで,Mはメサー, A は振幅,ωは周波数となる.ωはサーカディアンリズム の場合24時間を単位としたが,ωは2π/24となる.一 方,一年を月で表した場合は2π/12となる.
2)図1であらわした,数式をExce1のVisual Basicで 書き換えると,図2のようになる.
ここで,MyDataはオリジナルのデータ系列で,
Exce1(Microsoft)のsheetをクリックすることより求 めることができる.そのデータはmyData(3, i)に格 納される.Mytimeは,時間のデータであり,指示
(lnputBox関数)にしたがって,その領域(range)を入 力する.これらのデータの正弦成分はmyData(2, i)
に,余弦成分はmyData(3, i)に格納される. myPeriod は周期である.周期成分も指示(InputBox関数)にし たがって,領域を入力すると,最小自乗余弦スペクトル を求めることができる.三元一次方程式は,逆行列を 用いて,その解を求めた.その数値解析の部分は
市丸 雄平・小林 美佳子・沈 恵芳・安原 朋美
最小自乗余弦曲線の数値演算:
最小自乗法の理論的背景: 最小自乗余弦曲線の解法:求める余弦曲線は YニM+Acos(ωt一θ) ・……・・(1)
で表すことができる。ここで、時間t(i)におけるYの値をy(i),誤差をE(i)とすると、Y(i)は Y(i)=M+Acos(ωt(i)一θ)十e(i)・・……・(2)
となる。ここで、Y(i)はt(i)における実測値、 Mはメサー、 Aは振幅、ωは周期、θはアクロフェー ズ(極値)、e(i)は残差である。この式は、
Y(i)=M+Acosωt(i)・cos(θ)十Asinωt(i)。sinθ+e(i)……(3)
と変形することができる。ここで、誤差e(i)は
e(i)ニY(i)−M−Acosωt(i)・cosθ一Ashlωt(i)・sinθ……(4)
となり、誤差の自乗和
Σ(e(i)2)ニΣ((Y(i)−M−Acosω t(i)・cosθ一Asinωt(i)・sinθ)2)・……(5)
が最小になるようにM、Acosθ、Asin(θ)を設定すればよい。
ここで、βニAcosθ、γ =Asin eとすると、(5)式は、
Σ(e(i)2)=Σ((Y(i)−M一βcosωt(i)一γsh1ωt(i))2)…・・…(6)
と変形され、この式をM、β、γで偏微分し、その値が0となるように、M、β、γの値を求める。
(6)式をそれぞれ、M、β、γで偏微分すると、
δΣ(e(i)2)/δM= Σ(2(Y(i)−M一βcosωt(i)一γ sin to t(i))x(−1))・……(7)
δΣ(e(i)2>/δβ= Σ(2(Y(i)−Mこβcosωt(i)一γsinωt(i))x(−cosωt(i)))...(8)
δΣ(e(i)2)/δγ= Σ(2(Y(i)−M−一 B cos ul t(i)一γsinωt(i))x(−sinωt(i)))…(9)
となり、次の3元1次方程式を解くと、M,β、γが得られる。
Σ(Y(i))= ΣM一Σcosωt(i)・β 一Σ sin a) t(i)・γ……(10)
Σ(Y(i)・cosωt(i))ニΣM・cosωt(i)一Σcosωt(i)・cosωt(i)・β一Σsinωt(i)・cosωt(i)・γ…・・(11 Σ(Y(i)・sinωt(i))=ΣM・sin to t(i)一Σcosωt(i)・sinωt(i)・β一Σstnωt(i)・sinωt(i)・・Y・・…(12)
図1最小自乗余弦法の数学的背景
mylnvmatの配列で示した.この式よりM,α,βが 得られる.さらに,αおよびβの値より,次式を用いて 振幅(A),およびω(Acrophase)を求めることができ る.すなはち,
ω=arctan(α/β)A。再
となり,α,βの符号によりωまたはω+πをその解と
して求めることができる.
得られた余弦曲線の有意性にっいては,確率を直接法 でもとめた6).式は
P.・=(Σ(y,yh。t)2/Σ(y−yb。,)2)(n−3)/2
となる.ここでyhatは最適余弦曲線で求めた推定値で あるybarは測定データの平均値である.
Visual BASIC for Applicationによる数値演算
For i=OTo ihyDatanurnber−1 myClockニmytime(i)
myData(0, i>「1
myData(1, DニCos(2*PAI*myClock/myPeriod)
myData(2, i)ニSin(2*PA1*myClock/myPeriod)
myData(3, i)=Cells(myrow 1+i, mycolumn l)
Sheets(mysheetname).Cells(i+2,1)ニmyClock
Sheets(mysheetname).Cells(i+2,2)ニSheets(mysheetname).CeUs(myrow1+i, mycolllmll l)
Next I
逆行列演算
For iニ1To 3
For j=1To 4 mylnvmat(i, j)=O For kニOTo myNum
mylnvmat(i, D;mylnvmat(i, j)+myD ata(i−1, k)*myData(j−1, k)
Next k Next j Next i
図2最小自乗余弦曲線をExce1のVisual Basic Applicationを用いたマクロのプログラム
80 60 40 20
−20
−40
−60
−80
図3−1 胃電図の原形
29
>ll
311
藪{1
8111417202326293235384144
周期(秒)
図3−2 最小自乗余弦スペクトル
市丸 雄平・小林 美佳子・沈 恵芳・安原 朋美
結 果
鑓 , 層…華纏」醗蝋麗難蓑灘藤灘灘蓑羅1妻霧懸菱鑓蒙識諜鑓
.iiミ・
蝶鑛蕪灘螺欝難馨難慧懇i難灘翻灘畿
RgG7 叢・醐§掘
灘
灘鶏蕪,§・灘,醤鑑繋・霧蓑灘華.繋鱗・難蕪譲騨舗鑛彗璽E滋 2@ −22 052i O.89 Q988177
董 羅 、 ゜°@伽 響・鱒 ゜ ,° ・・3° ・ 9° ■ .・ . °° 、・9 9 ,°,°9り,鯛 7,: 響゜ り ,. 曹 曹゜曹..f響
、・ 、
Di騰
, r , ρ
@ 4: −219 144; −02213188597 ,・ ・ , , ▼ く 9,。, ● ,, ▼ 6, ・ , 19 , .., , D.電,,,
@ 5 −22 368 138・0549806 .韓;. 。,,,..__. 響●7■■77
℃、E雛鍵
?E .潤。o・. ・ ., 7●9
@ 6・ −218 324: 062 4012865
@ 7 −1 94 1022 一つ58i OOO879 7冒甲薗,,
@ 9
ミゴ ζ゜°ζ 、 , , u , , . , ,電
@ 0280973 , 、
Ψ、?Eく
i:豊 、
・ , ・ 9 ,
C_ユ9 ,.二217三 _5.B4、 . −a29..⊆㌧〜,5Z2窪.
曾蓬サ
・甲
@ 11 −・258 6.54 −172 257.891 、
、 嚇 、
m・
, 、 , 7 .,
│157 1415 −043 95E−05 . r77 , 、・・ o ,,. 一 幽 ●o ・ ● .
、魏・
13・, 9 7 0 ・ ,
@ 14 −183
@ 15 −22
塁議.・駕︑ 1083 −127 0005203
@ 258 −111 0745032・,7,7,,▼●曹,曾暫o,,騨▼響.響響.▼
@ 16 一書59 122 −031 0001341
難懇. 1
P8 −0◎6 3199・ −143 423E−27 ♂7,0●o響●r・顧 難; ,●7,,曾●o●曾,甲・7.9@ 19 −068 3738 −1 97° 214E−43
@ −25
、
V 噂
…冒ナ・
D
■●9,尋,■・璽 ■・
@ 21
嚢{
窪q1.興
@ 23 r565 3705. −377 154E−39
蕪聾灘・o■⁝
尊
秩D.
く香A灘︐﹂..響99
Q8 −3 1755 129 493E−07 甲 . , . 量 , , , , 8 7 , ● 「 , , . , , 「
,
777
Q9 −254 1437 11 735E−05 、 、 、、 、、、 、、、、
.ii曝 }噤Aい…、甲
@ ,
9ζ゜・… °
^
,、 亀A、、、 一A、−
@ 30 −223° 1155 095・、OQg2283 , 曹 、 , , 璽
@ 3! −207. . 9」2 . . Q」8@ Og22031
@ 32. −206、 719 −229 1021741 Q 、_.3i…_._二aユヱ.___.§.貫塾、_ 二ala. 1〜塾9農53.
墨,.,・,.σ
図3−3 最小自乗余弦スペクトルを求あるためのシート
80 60 40
320
こ}o
藩一2・
:19
−80
1痔間 (秒)
図3−4 胃電図原波形余弦曲線 y=−5.02十42.04cos(2πt/21十2.99)
ExcelのVlsual Baslc for Applicationのマクロを 用い記述した最小自乗余弦曲線法の関数を時間栄養学に 応用した.
1)ウルトラディアンリズム
ゥルトラディアンリズムの例として,図3に胃電図現 象を示した.3−1に胃電図の原波形,3−2には最小自 乗余弦スペクトル,3−3には最小自乗余弦スペクトル 解析のためのディジタルデータ系列,3−4には,最小 自乗余弦スペクトルよりもとめた最適余弦曲線と原波形 を示した.最小自乗余弦スペクトルはシートの1行1列 より32行1列まで,レンジを設定すると,2−5行目に それぞれの周期に対応するMesor値,振幅, Acrophase,
および確率が,計算の結果自動的に提示される.3−1 に対応する最適余弦曲線の周期は,21秒であり,その Mesor値は5.02,振幅は41.04であった.これにより,
次式に示される余弦曲線は現波形のもっとも適合する余 弦曲線(最適余弦曲線)であることが推定された.
y;−5.02十42,04cos(2πt/21)
2)サーカディアンリズム(図4)
血圧および心拍は24時間リズムを示すことがよく知ら れている.図4は携帯型24時間血圧計で計測した収縮期 血圧および脈拍の24時間変動にっいてしめしたもので
ある.
血圧:Y=118+9cos(2πt/24−3.69)
心拍:Y=59+8cos(2πt/24−3.68)
血圧および心拍の極値は14:00と,同位相にあった,
3)インフラディアンリズム
第5図は栄養摂取状況を2週間にわたって計測し,エ ネルギー摂取量を計算したものに,最適余弦曲線をあて はめたものである.余弦曲線は,
y=・1019.2十10.9 cos(2π/7−O.47)
となり,エネルギー摂取は統計学的に7日間のリズム変 動を認めた.
考察
生体は種々リズム成分が認められるこのリズムは,外 的環境に内的環境を積極的あるいは受動的に同調させて いるものと解釈される。したがってリズムを解析するこ とにより生体内の制御システムを推測することも可能で ある.リズムの解析は従来より,フーリエ解析,最大工
1
8:1 8:1
$IQi5
1:l l:1
0
088
8寸
確率スペクトル
Nゆ
寸 卜 cりト・ 寸 σ》
@お . ●
o
寸◎o ゆ
o◎︾ oco
d
▼一 卜oo
卜o
呂o
d d
Nゆ
oり ◎◎
d
◎o
o
お
o
めoo o①oq
}oPo o
3
4
5 67 8 9
時間(日)
10 11 12 13
Y
ミ耗H
エネルギー摂取の週内リズム
8 1011121314
時間(日):(y=1 Ol 92+109pos(2π/7−O・47)P:0・0035)
図4 収縮期血圧と心拍の日内リズム 収縮期血圧の最適余弦曲線:Y−118+9cos(2πt/24−3.69)
心拍の最適余弦曲線:Y三59+8cos(2πt/24−3.68)
市丸 雄平・小林 美佳子・沈 恵芳・安原 朋美
収縮期血圧と心拍の日内変動
180 160
癒140 澤120
.bSloo
80 60 目 40 20 0
0.25 2 4.25 6 8.25
10 11.813.515.8
時計時間(時)18 19.8 21 22.8
図5 エネルギー摂取量(栄養計算で求めたもの)の確立スペクトル(上段)と7日間周期の変動(下段)
ントロピー法,余弦曲線法が使用されてきたが,余弦曲 線法はリズムとくに生体のリズム成分を検討する方法と
して広く使用されている.われわれは,最近,この方法 を用いて,ラットのエネルギー摂取比率において,リズ ムがあることを示し,とくに,糖尿病ラットにおいては,
脂質と炭水化物のエネルギー比に位相のずれがあること をあきらかにし,その病態生理にっいて検討した6).こ の現象は,リズム解析によってはじあて明らかにできた ことである.しかし,解析システムの構築はすでに開発 していたUNIXを用いて行なってきたため,汎用性に 欠けることが問題点であった.今回,表計算ソフト
(Excel:Microsoft)を用いておこなった表計算システ ムの特徴は,対話型であり,データ選択が容易であるこ とが挙げられる.今回のシステムでは,データ領域,時 間領域,求める周期の領域(スペクトル計算),および 結果データの表示領域を対話的に表示した.結果におい て示されたように,時間領域の表示は任意の時間が可能 であり,ウルトラディアン,サーカディアン,およびイ ンフラディアンディアンリズムのいずれの解析も可能と
した.
ウルトラディァンリズムの例とした胃電図は従来より 20秒のリズムがあることが報告されている7).今回の解 析で,胃電図が21秒のリズムであることが推定された.
また,今回の胃電図は連続的に解析が可能であるため,
日内変動を加味した,胃電図周期にっいても検討中であ
る.
まとめ
Excel(Microsoft社)のVisual Basic for Applica−
tionを利用して,生体リズム(ウルトラディァン・リズ ム,サーカディアンリズム,およびインフラディアンリ ズム)の解析システムを構築した.このシステムは,対 話的にリズムを解析することが可能である.このシステ ムを時間栄養学に応用する可能性として,胃電図,血圧・
心拍,およびエネルギー摂取のリズム解析例を提示した.
参考文献
1)Nelson W., Tong Y.L, Lee K−K., Halberg F.:
Methods for cosinor・rhythmometry.
Chronobiologia 6:305−322,1979.
2)Vagnucci A.H., Wong A.K.C., and Liu T,S.:
Time series analysis of hormonal patterns in human plasma.
Computer and Biomedical Research 7:513−532,
1974.
3)Halberg R, Lagoguay M。, and Reinberg A.:
Human circannual rhythms over a broad spec−
trum of physiological processes.
Inter. J. Chronobiology 8:225−268,1982.
4)Monk T.H。, and Fort A.
℃osina Acosine curve fitting Program suit−
able for small computers.
Inter. J. Chronobiology,8:193−224,1981
5)Ichimaru Y.
Multivariate cosine spectrum analysis for am−
bulatory blood pressure and heart rate.
Therapeutic Res.14:194−201,1993
5)佐々木 隆周期成分の探索,時間生物学,
pp.312−332,朝倉書店,東京
6)Ichikawa M., Kanai S., Ichimaru Y., Funakoshi A.,Miyasaka K.:The diurnal rhythm of energy expenditure differs between obese and glucose −intolerant rats and streptozotocin−
induced rats.
J.Nutrition.,130:2562−2567,2000
7)Alvarez W.C.:The electrogatrogram and what it shows. J.A.MA.:78:1116−1119,1922.
Abstract
We have developed a system fbr the analysis of bioIogical rhythln by using Visual Basic fbr Application(Microso丘).
By using dialogue method or by simply selecting the range of data−, time−and periodicity−table written by Exce1, we can easily analyze ultradian−, circadian−or infradian−rhythm of the biological data. Furthermore, mesors, amplitudes, acrophases fbr each period were written automatically on the table which enables us to visualize tho data on chart easily. In this paper, we have shown illustrative examples of time−sequential伽一series of infradian rh抽m of electrogastrogram,
circadian rhythm of heart rate and systolic blood pressure, and in丘adian rhythm of fbod intake.