Rのパッケージ、RHRVを用いて心拍変動解析
今日は心拍変動解析のお話。
心拍変動解析は心電図のRR間隔(RRI)のゆらぎをスペクトル解析をもちいて、
周波数帯域を高周波HF、低周波LFとし、自律神経の副交感神経の活動をHFが、
交感神経の活動をLF/HFが反映することを利用した解析のことです。
今回は実践的にコンピューター言語であるRで解析できるRHRVパッケージが
2019.10.28にupdateされていたので、健忘録としてまとめました。
公式HP:http://rhrv.r-forge.r-project.org/
まずRにインストール
install.packages("RHRV")
library(RHRV)
読み込めるのは欧米フォントのみの様々なファイル形式に対応しています。
でもexample.beatsをNotopadにコピーして、CSVファイルかTXTで保存します。
RHRVでは特殊なデータ構造を構築して解析していくようですが、詳細は公式に。
まず、箱を作るためのおまじないです。
data("HRVData")
data("HRVProcessedData")
hrv.data = CreateHRVData()
hrv.data = SetVerbose(hrv.data, TRUE )
解析している内容としては
①RRIデータを読み込む
②心拍数へ変換する。
③等間隔のデータで解析するために線形補間で再構築する。
④フーリエ変換かウーブレット変換でスペクトル解析する。
⑤HF、LF、VHF、VLFを計算して図示する。
まずRR間隔(RRI)データを読み込みます。
hrv.data = LoadBeatAscii(hrv.data, " example.beats")
(ここでの落とし穴はこのexampleデータの内容でした。
どうもRRIの積算データになっているようです。
もし用意したデータが心拍数データであれば、変換が必要です。
心拍数=60 ÷ RRI(msec) なので、 RRImsec = 60 ÷ 心拍数として、
1拍目から累積データを作成して、読み込んでください。
だれかうまい方法があればお知らせください。)
②心拍数へ変換する。
hrv.data = BuildNIHR(hrv.data)
素の心拍数データを作成してます。
外れ値を削除するためにフィルターを自動でかけます。
hrv.data = FilterNIHR(hrv.data)
③等間隔のデータで解析するために線形補間で再構築します。
hrv.data = InterlolateNIHR(hrv.data, freqhr =4)
デフォルトは4Hzです。
素の心拍数を図示します。
PlotNIHR(hrv.data, main = "niHR")
④ タイムドメイン解析
hrv.data = CreateTimeAnalysis(hrv.data, size =300, interval = 7.8125)
④フーリエ変換かウーブレット変換でスペクトル解析する。
今回はフーリエ変換
hrv.data = CreateFreqAnalysis(hrv.data)
⑤HF、LF、VHF、VLFを計算して図示する。
デフォルトで行うと
hrv.data = CalculatePowerBand(hrv.data, indexFreqAnalysis = 1, size = 300, shift = 30)
PlotPowerBand (hrv.data, indexFreqAnalysis = 1, ymax= 200, ymaxratio = 2.0)
なにが行われているが不明な点が多いけど、高額なsoftが不要で、心拍変動解析ができるので、連続する心拍数さえあれば、いろいろ解析できて、役立つ可能性があります。
■
昨日のつづき
線形回帰で解析します。
> lm=lm(PRdelta~Medicine,BB)> lm=lm(PRdelta~Medicine,BB)> summary(lm)
Call:lm(formula = PRdelta ~ Medicine, data = BB)
Residuals:
Min 1Q Median 3Q Max
-24.0000 -6.0326 0.8696 5.8696 18.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.130 2.078 5.356 2.38e-06 ***
MedicineT 1.870 2.828 0.661 0.512
---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 9.967 on 48 degrees of freedom
Multiple R-squared: 0.009021, Adjusted R-squared: -0.01162 F-statistic: 0.437 on 1 and 48 DF, p-value: 0.5117
古典統計学では95%信頼区間から統計学的有意差なし。p値0.512
次はmcmcでベイズ推定 モデルは Xi1=a+b*Xi1+ei xi2=a+b*Xi2+ei ダミー変数を設定して > model<-" + data{ + int<lower=0> M; + int<lower=0> N; + real x1[M]; + real x2[N]; + } + parameters{ + real<lower=0> s1; + real<lower=0> s2; + real m1; + real m2; + } + model{ + for(i in 1:M) + x1[i] ~ normal(m1, s1); + for(i in 1:N) + x2[i] ~ normal(m2, s2); + }
つぎはbrmsパッケージを用いてマルコフ連鎖モンテカルロ法によるベイズ推定
> brm=brm(
+ formula=PRdelta~Medicine,
+ family=gaussian(),
+ data=BB,
+ )
> print(brm)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: PRdelta ~ Medicine
Data: BB (Number of observations: 50)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept 11.10 2.11 7.06 15.20 1.00
MedicineT 1.88 2.87 -3.76 7.59 1.00
Bulk_ESS Tail_ESS
Intercept 3715 2945
MedicineT 3809 2881
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma 10.13 1.05 8.35 12.31 1.00 3516
Tail_ESS
sigma 2552
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
2群間の平均の差は1.88程度と推定されている。
■
> X2=subset(BB,Medicine=="s")$PRdelta
> M=nrow(subset(BB,Medicine=="T"))
> N=nrow(subset(BB,Medicine=="s"))
t = 0.68164, df = 45.735, p-value = 0.4989
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.652168 7.391298
sample estimates:
mean of x mean of y
13.00000 11.13043
+ data{
+ int<lower=0> M;
+ int<lower=0> N;
+ real x1[M];
+ real x2[N];
+ }
+ parameters{
+ real<lower=0> s1;
+ real<lower=0> s2;
+ real m1;
+ real m2;
+ }
+ model{
+ for(i in 1:M)
+ x1[i] ~ normal(m1, s1);
+ for(i in 1:N)
+ x2[i] ~ normal(m2, s2);
+ }
+ generated quantities{
+ real diff;
+ diff=m1-m2;
+ }
+ "
> print(fit1)
post-warmup draws per chain=500, total post-warmup draws=2000.
s1 12.10 0.04 1.76 9.25 10.86 11.92 13.11 16.15 1884 1
s2 8.23 0.03 1.35 6.09 7.26 8.04 9.02 11.29 1936 1
m1 12.96 0.06 2.36 8.27 11.41 13.00 14.49 17.59 1773 1
m2 11.11 0.04 1.71 7.71 9.99 11.13 12.19 14.47 1660 1
diff 1.85 0.07 2.88 -3.82 -0.06 1.76 3.76 7.69 1800 1
lp__ -134.65 0.04 1.43 -138.30 -135.35 -134.33 -133.59 -132.83 1033 1
+ data{
+ int<lower=0> M;
+ int<lower=0> N;
+ real x1[M];
+ real x2[N];
+ }
+ parameters{
+ real<lower=0> s1;
+ real<lower=0> s2;
+ real m1;
+ real m2;
+ }
+ model{
+ m1~normal(0,10000);
+ m2~normal(0,10000);
+ s1~cauchy(0,5);
+ s2~cauchy(0,5);
+
+ for(i in 1:M)
+ x1[i] ~ normal(m1, s1);
+ for(i in 1:N)
+ x2[i] ~ normal(m2, s2);
+ }
+ generated quantities{
+ real diff;
+ diff=m1-m2;
+ }
+ "
> print(fit0)
post-warmup draws per chain=500, total post-warmup draws=2000.
s1 11.76 0.04 1.73 9.05 10.48 11.56 12.82 15.64 1605 1
s2 7.91 0.03 1.23 5.97 7.06 7.75 8.64 10.82 1797 1
m1 13.05 0.05 2.21 8.64 11.58 13.02 14.54 17.46 1827 1
m2 11.13 0.04 1.67 7.83 10.00 11.14 12.23 14.41 1714 1
diff 1.92 0.06 2.76 -3.51 0.23 2.01 3.68 7.36 2189 1
lp__ -137.69 0.05 1.45 -141.29 -138.42 -137.37 -136.63 -135.87 850 1
■
今日はニューラルネットワークです。
といっても3層パーセプトロンをnnetの構築です。
ユニット層は3としています。
昨日と同様に教師データと評価データを分けて開始です。
#乱数のもとを固定します。
set.seed(1)
#ニューラルネットワークのモデルを構築します。
nn.model<nnet(eGFR2低下群~.,data=eGFR.train,size=3,maxit=300)
# weights: 160
initial value 38.368978
iter 10 value 36.379955
iter 20 value 27.414648
iter 30 value 22.653967
iter 40 value 17.436598
iter 50 value 16.306668
iter 60 value 15.248421
iter 70 value 8.661530
iter 80 value 4.821300
iter 90 value 4.418075
iter 100 value 4.417687
final value 4.417682
converged
#80回を超えたあたりから収束しています。
#ニューラルネットワークで作成した分類モデルを評価データに適応
nn.pred<-predict(nn.model,eGFR.test,type="class")
#ニューラルネットワークによる予測結果を確認
head(nn.pred)
[1] "a" "a" "a" "a" "b" "b"
#予測の正誤をまとめた表を作成
(nn.tab<-table(eGFR.test$eGFR2低下群,nn.pred) )
nn.pred
a b
a 35 1
b 9 10
#予測精度の確認
sum(diag(nn.tab))/sum(nn.tab)
[1] 0.8181818
このモデルでは81.8%でした。
意外と単純ベイズの精度がよかった。
RとRstudioで機械学習
今日はRとRStudioで機械学習の備忘録
最近始めたRでの機械学習で医療系データを解析してみました。
データは健診データで目的変数はeGFRの悪化予測としました。
(R3.5.0を使用して、どんな程度できるかのお試しです。)
まずパッケージの読み込みで(e1071)と(kernlab)を指定します。
各説明変数を教師データと解析データに2分割していきます。
その後、単純ベイズモデルに適応しています。
感度特異度から予測精度を確認します。
実際のRの解析の流れはこんな感じです。
#パッケージの読み込み
library(e1071)
library(kernlab)
#健診データの読み込み
#(省略)
#教師データと評価データに分
n<-seq(1,nrow(Dataset),by=2)
eGFR.train<-Dataset[n,]
eGFR.test<-Dataset[-n,]
#単純ベイズで分類モデルを作成
nb.model<-naiveBayes(eGFR2低下群~.,data=eGFR.train)
#単純ベイズで作成した分類モデルを確認
nb.model
Naive Bayes Classifier for Discrete Predictors
Call:
naiveBayes.default(x = X, y = Y, laplace = laplace)
A-priori probabilities:
Y
a b
0.5357143 0.4642857
Conditional probabilities
#単純ベイズで作成した分類モデルを評価データに適応
nb.pred<-predict(nb.model,eGFR.test)
#単純ベイズによる予測結果を確認
head(nb.pred)
[1] b a a a b a
Levels: a b
#予測の正誤をまとめた表を作成
(nb.tab<-table(eGFR.test$eGFR2低下群,nb.pred))
nb.pred
a b
a 36 0
b 2 17
#予測精度の確認
sum(diag(nb.tab))/sum(nb.tab)
[1] 0.9636364
年齢・性別・血圧・ABI値・CAVI値など
すべて説明変数に入れ込んだらいい感じに
単純ベイズで96.3%の精度が出てます。
意外と臨床的意義があるかもしれないです。
機械学習ってすごいって思いますね。
次回はニューラルネットワークモデルや
サポートベクターマシンも試してみます。