大阪で働くメディカルデータサイエンティストの日記

医療系dataをRを用いて解析します。古典統計学、ベイズ統計学など使えるものは何でも。

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).

 

95%ベイズ信用区間からは有意差なし

2群間の平均の差は1.88程度と推定されている。

平均の差の推定シリーズ
古典統計学ベイズ統計学の比較してみました。
今回は手持ちのデータで薬剤sと薬剤Tを使用したときに
どの程度parameterであるPRdeltaが変化したかのデータです。
2群間の平均の差をいろんな方法で比較していきます。
まずデータをBBとして読み込みます。
その後、X1に薬剤Tを使用した群、X2に薬剤sを使用した群を割り付けます。
> X1=subset(BB,Medicine=="T")$PRdelta  
> X2=subset(BB,Medicine=="s")$PRdelta
> M=nrow(subset(BB,Medicine=="T"))     
> N=nrow(subset(BB,Medicine=="s"))
まずは通常通りWelch のt検定です。
> t.test(X1,X2)
Welch Two Sample t-test
data: X1 and X2
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
もちろん古典統計学では統計学的有意差なし、
非劣性マージンを”2”としたときに非劣性は証明できません。
 
次に、ベイズ統計学です。
マルコフ連鎖モンテカルロ法で平均の差Welch のt検定です。
データをlist形式にしてdatに格納します。
> dat <- list(M=M, N=N, x1=X1, x2=X2)
モデルを指定します。
今回は平均、分散は一様分布に従うと仮定しています。
> 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);
+ }
+ generated quantities{
+ real diff;
+ diff=m1-m2;
+ }
+ "
> fit1 <- stan(model_code=model, data=dat, iter=1000, chain=4)
> print(fit1)
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
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
 
このモデルでは95%ベイズ信用区間からは有意差は認めません。
2群間の平均の差は1.85程度と推定されています。

さらに別のモデルを指定します。
平均は幅の広い正規分布、分散は半コーシー分布を指定しています。
> model0<-"
+ 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;
+ }
+ "
> fit0 <- stan(model_code=model0, data=dat, iter=1000, chain=4)
> print(fit0)
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
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

このモデルでは95%ベイズ信用区間からは有意差はないと判断されます。
また2群間の平均の差は1.92値程度と推定できます。
 
次回は線形回帰で確認していきたいと思います。

今日はニューラルネットワークです。

といっても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%の精度が出てます。

意外と臨床的意義があるかもしれないです。

機械学習ってすごいって思いますね。

次回はニューラルネットワークモデルや

サポートベクターマシンも試してみます。