統計ブログはじめました!

各専門分野の統計技術、方法、テクニックなどを気ままに分かり易く例題をもとに解説します。

医学と統計(49)

2010-07-15 11:01:56 | 日記・エッセイ・コラム

情報統計研究所へのアクセスはここからお気軽に

自己相関係数について。
前回は移動平均と時差相関を取り上げました。今回は、自己相関について考えてみたいと思います。時差相関では、末梢血中好塩基球数(Ba)と血清IgE値の2群間のタイムラグでしたが、ここでの
自己相関は心電図波形から計測したR-R間隔時間(msec:以下、RRと言う)を例題としています。通常、健常者であれば呼気と吸気時でRRが異なります。
図1は健常者の通常の呼吸時におけるRRの頻度を幹葉図で示しました。
 図1 通常呼吸時RR(SAS-JMP)
Photo

そして、
強制的にゆっくりとした深い腹式呼吸(約10回/分)をおこなった時のRRの頻度を図2に示しました。
 図2 強制呼吸時RR(SAS-JMP)
Photo_2

通常呼吸時と強制呼吸時の統計量は次の通りでした。
                          通常呼吸時        強制呼吸時
平均値  =  915.65 msec       914.67 msec
標準偏差=    41.52 msec          73.75 msec
最頻値  =   920     msec        830      msec

ここでの、
通常呼吸時と強制呼吸時におけるRRの自己相関係数は、図3と図4のようになりました。
これをコレログラムと言います。
 図3 通常呼吸時のコレログラム(SAS-JMP)
Photo_3

 図4 強制呼吸時のコレログラム(SAS-JMP)
Photo_4

強制呼吸時のコレログラムは規則的に循環変動していることが分かります。

自己相関係数の計算は前回の時差相関係数と同じで、RRを一つずつずらして計算したものですので、MSエクセル関数の「SUMSQ」と「SUMPRODUCT」を使って、前回の時差相関係数と同様の計算をすれば良いでしょう。

「R」なら、acf ( time series data ) で図5の結果を得ることができます。
しかし、
acf ( time series data ) の”time series data” は時系列データオブジェクトでなければなりません。もし、
RRのデータが例えば、次の様に並んでいるだけなら、
 dat <-  c ( 860,930,890,911,1000,891,・・・・ )

これを、
例えば、t_dat <-  ts ( dat, frequency=1 ) としてデータオブジェクトを作成しておき、
  acf ( t_dat )

とすれば、図5のような自己相関プロットが得られます。
 図5 「R」による自己相関プロット図
Acfplor

SAS-JMP なら、「分析→モデル化→時系列分析」で求めることが出来ます。

RRは副交感神経機能検査などで用いられており、ここでの自己相関分析も一つの検査方法として考えてみて下さい。

 


医学と統計(48)

2010-07-04 11:14:29 | 日記・エッセイ・コラム

情報統計研究所(統計分析のご相談ご依頼)のアクセスはここから。

前回の末梢血中好塩基球数( Ba ) と血清IgE値( IgE ) の実測データを用いて、実際に、時差相関係数を求めてみましょう。
Ba と IgE の測定単位が違いますので、次により、標準化(正規化変換)を行います。すなわち、
 Ba ( Xi ) の mean±SD=125.35±32.62、
 IgE ( Yi ) の mean±SD=583.48±183.42、
 
から、
 ( Xi-125.35 ) / 32.62 、( Yi-583.48 ) / 183.42

により標準化(正規化変換)したデータを図1に示します。
図1 標準化した時系列データ
Normalydata

図1のデータについて、図2のようにデータを一つずつずらせて
 r[k]=Sxy[k]/SQRT( Sxx[0] ・Syy[0] )

を計算しますと、時差[k](k=Lag Time)の相関係数が求まります。。
図2 時差相関係数の求め方
Lagcorrenormaly

BaとIgEの時差相関係数は次の通りです。
 r[0]=0.537、r[1]=0.549、r[2]=0.552、・・・、r[10]=-0.149

Ba とI gE の時差はLag=2 で最大であり、Ba の増加と IgE の上昇には2日の時差があるようです。
さて、MSエクセルの関数を使っての計算が面倒なら・・・、「R」を使用すれば簡単に求める事が出来ます。
「R」では ccf() が用意されており、
 ccf ( X , Y )

で図3 の出力結果が得られます。
図3 「R」による時差相関プロット
Ccfbaige 

時系列データでは、単相関で相関が見られなくても、時差を与える事で相関関係や循環性などの情報が得られる事もあります。

次回は自己相関について、心拍(R-R間隔)を取り上げたいと思います。