Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2018-03-17

心不全の心血管イベント発症率を予測するモデル()

読んだ。

The impact of creating mathematical formula to predict cardiovascular events in patients with heart failure.

Sci Rep. 2018 Mar 5;8(1):3986.

プレスリリース

COI:なし

 

慢性心不全患者の血液検査、診察所見、検査データなどから、心血管イベントによる再入院もしくは死亡の確率を高精度に予測する数理モデルを作った、という話。

 

モデル探索では、402個あるデータのうち252個を選んで、167人の患者で心血管イベントが起こるまでの日¥tau について、パラメータセット¥bf{x} により

¥tau=f(¥bf{x}) …※

という単なる回帰モデルを考える。これにより、患者i のイベントが起こる確率は指数分布でモデル化して

p_i(t)=¥frac{1}{¥tau_i}¥exp(-¥frac{t}{¥tau_i})

になる。このモデルパラメータ¥tau はそのままカプランマイヤーパラメータなので、曲線が引けるようになる。

(152人の患者に16人足して167人が対象、というのが意味不明)

パラメータセットは多すぎるのでL1 正則化によって最終的に50個のパラメータを選び出してきている。このパラメータセットと、回帰係数は表2にまとまっている。

 

validation セットの213人でこの50個モデルを使ってKM曲線(というよりパラメータが決まっているのでもはや関数としての曲線)を実際のデータと比較して、goodness-of-fit 検定をして決定係数を有意かどうかの指標、にしているらしい。

これで比較して「KM と予測曲線はsignificantly に近く、決定係数はP=0.0784」と図2に書いてある。

 

P の意味するものが不明確である。p値であるとするならば、goodness-of-fit で棄却されなかったことを持って「当てはまりがよかった」というのは、帰無仮説検定としてはダメである。

P が決定係数であるとするならば、一般に決定係数が意味するところを考慮すると、推定曲線は実際のデータから推定されたKM を7.84% しか説明しない、ということになる。

いずれにしろ、「心血管イベントの発生確率を予測」には不十分と思う。

 

さて、予測する数理モデル、と言っているので、実際のモデルを眺めてみる。やっていたことは、※で表される式の正則化で、心血管イベントまでの日数を推定しているので、(適当に変換がはいって)係数が正ならば予後が悪い、負ならば予後がいいほうへ寄与する。例えば、虚血性心疾患(IHD)に起因する心不全予後不良(10.842)、家族が多いと予後良好(-0.525)などが本文中で言及されている。

これらの中で、どう考えてもおかしいのがBNP(-0.826)である。退院時にはBNP が高いほうが心血管イベントに対しては予後がいい、という結果になっている。循環器内科としてこれを見逃してるのはやばいと思う。BNP が高いと予後不良と同グループが言っているが、これはいったいどういうことなのだろうか() また、下腿浮腫(-0.692)というのもなんだかなあ、という印象。

筆者らはdiscussion で「各パラメータは複雑に絡み合っているため、ここで捉えるというより50個のパラメータネットワークで考えるべき」みたいなことを言っているが、50個のパラメータセットを提唱するならばそれなりのネットワークを考察するべきではなかろうか。そうでなければ「心房細動があって脳梗塞があってCOPDARCRP/AST/BNP が高くて抗炎症薬、抗甲状腺薬、胃腸薬、PPI、鎮静薬を飲んでいると予後良好」ということになるが、これは一体どういうことなのだろうか。

 

ということで、上記で言っていることは多重共線性に起因すると思われるのだが、L1 正則化では効果があるか、ゼロか、で変数選択が行われがちなので、相関がある変数セットのなかでは、ただひとつの変数選択される、らしい。

(ぐぐってブログをいくつか流し見した程度なので理論的な保証はここではまったくない)

しかし、ある変数が3つ以上の相関している変数セットに含まれているとき、いったいどの変数正則化で選ばれるのかは、不勉強なのでよくわからないし、選ばれた50個の変数を見る限り、浮腫とBNP は相関高いと思われるし、なんだかなあ、という印象。

 

個々の心血管イベントを予測、という割には、途中でKM と予測曲線の分布の差をKL の最小にするようなパラメータを求めていたりして、本当に個人の予測なのか、集団として生存曲線をfit させているのか、ちょっとよく理解しきれなかった。

2017-11-11

面白そうな統計の話・論文

Crowdsourced research: Many hands make tight work

黒人サッカー選手が警告を受けやすいかを解析すると、様々な結果が出てしまうという話。

 

Avoidable waste in the production and reporting of research evidence

適切な疑問を設定すること、適切な研究デザイン・手法を採用すること、結果にアクセス可能にすること、バイアスのない報告をすること、を守らないと実に85%の研究は無駄、というような話。

 

The emerging field of mobile health

生体データを日常的に連続して取れるようになっている話。

 

Timing of onset of cognitive decline: results from Whitehall II prospective cohort study

どんな認知機能が加齢に応じて衰えていくかの話。reasoning (論理的思考)の落ちがやばい。

 

Are Observational Studies More Informative Than Randomized Controlled Trials in Hypertension?

高血圧研究において、RCT よりも観察研究のほうがよさそうな場合とは?

RCT が真の臨床状況を反映していない場合、報告バイアス、メタアナリシスの収集バイアスなど。

 

HASE: Framework for efficient high-dimensional association analyses

高次元データの相関を解析するときに、重複した組み合わせを除外して…とか言っていた。積んでる。

 

Voxelwise genome-wide association study (vGWAS).

MRI のボクセル単位でデータを取得してGWAS をするという話。よくわかっていない。

 

A polymorphic DNA marker genetically linked to Huntington’s disease

家系解析でハンチントン病の原因遺伝子を4番染色体に見つけた話。

Lin-Gettig syndrome

Joubert syndrome

Rubinstein-Taybi syndrome

 

Mendelian randomization

2017-11-09

 

NGS データの解析

no title

2013-01-11

MikuHatsune2013-01-11

世界で一番わかりやすい心電図ベクトルを目指す

心電図というものは各心筋の電気活動を心臓という臓器で見た全体がもつ電位ベクトルというものを、体の各部位に取り付けられた電極で観測してグラフ化しているという。

電極の付け方で、胸部誘導という、6つの電極をつけるやり方があるのだが、これは胸を水平に切ったとして、そのときの心臓の電気活動を記録する。

ベクトル心電図はカージオイドとなり、こんな感じの記録が取れる。

心臓の形から実際にどんなカージオイドっぽい形(関数)が取れるのかはめんどくさすぎるので、今回はa=1とした簡単なカージオイドからベクトル心電図→胸部誘導心電図の変換をしてみる。

極座標表示ではx軸から反時計回りを角の正方向とするが、胸部誘導の置き方は第3, 4象限になる。

各誘導と心ベクトルのなす角をもううろ覚えな内積の公式からゴリ押しで求めつつ、原点から各誘導に向かうベクトルに正射影ベクトルを作り、原点から誘導に正射影ベクトルが向いているなら正の電位が、反対を向いているなら負の電位が観測される、という仕組みである。

#極座標で考えつつxy座標にもする。
Vs <- c(4/3, 3/2, 19/12, 5/3, 11/6, 0) * pi #胸部誘導の置く角度
ra <- 4 #胸部誘導を置く遠さ。あまり関係ない。
Vs_xy <- cbind(ra*cos(Vs), ra*sin(Vs)) #胸部誘導の座標
theta <- seq(0, 2, length=10000) * pi #心ベクトルの回転
a0 <- 1 #心ベクトルの長さ
#心ベクトルの座標
#カージオイドを採用する。
x0 <- a0*(1+cos(theta))*cos(theta)
y0 <- -a0*(1+cos(theta))*sin(theta)

#カージオイドを回転させたかったら
kaiten <- function(vec, psi){
	rot <- matrix(c(cos(psi), sin(psi), -sin(psi), cos(psi)), 2, 2)
	c(rot %*% vec)
}
rot0 <- 0 # 1/8*pi
#回転させないけれども0度回転させておく。
xy <- apply(cbind(x0, y0), 1, kaiten, psi=rot0)
x0 <- xy[1, ]
y0 <- xy[2, ]

#カージオイドがなぜかうまく描けない。描けるところをとってくる。
p0 <- 10 * (seq(ceiling(length(theta) / 10)) - 1) + 1
card0 <- c((length(p0)/2):1, rev(tail(seq(p0), length(p0)/2)))
x1 <- x0[p0][card0]
y1 <- y0[p0][card0]

#心ベクトルがx軸となす角度を求める。
acos1 <- acos(x1 / sqrt(rowSums(cbind(x1, y1) ^2)))
acos2 <- rep(0, length(y1))
for(l in seq(y1)){
	acos2[l] <- c(acos1[l], 2*pi - acos1[l])[1 + (y1[l] < 0)]
}
#V1-V6までの誘導と、心ベクトルがなす角を一気に求める。
ECGs <- t(cos(mapply(rep, acos2, length(Vs)) - Vs)) * sqrt(rowSums(cbind(x1, y1) ^ 2))

#プロットのレイアウト
L0 <- matrix(c(8, 2, 1, 3, 6, 4, 7, 5), nr=2, nc=4)
#胸部誘導のの付け方は赤黄緑茶黒紫と順番が決まっている。
leadcolor <- c("red", "yellow", "green", "brown", "black", "violet")
for(i in seq(x1)){
	#png()
	layout(L0)
	par(mar=c(4, 2, 2, 1)* 0.3)
	plot(x1, y1, type="l", xlab="", ylab="", axes=FALSE)
	abline(v = 0, h = 0, lty=2)
	arrows(x0=0, y0=0, x1=x1[i], y1[i], length=0.1)
	points(x1[i], y1[i], col=4, pch=16, cex=2)
	for(k in seq(Vs)){segments(0, 0, Vs_xy[k, 1], Vs_xy[k, 2], lty=3, col=leadcolor[k])}
	
	#各誘導の電位
	for(j in seq(Vs)){
		plot(ECGs[, j], type="l", xlab="", ylab="", axes=FALSE, ylim=range(ECGs), col=leadcolor[j])
		legend("topleft", legend=paste("V", j, sep=""), bty="n", cex=1.8)
		points(i, ECGs[i, j], col=4, pch=16, cex=2)
	}
	#空欄埋め用
	plot(0, type="n", axes=FALSE, xlab="", ylab="")
	#dev.off()
}

f:id:MikuHatsune:20130110223023g:image

V1では上がって下がる

V3は上がって下がるのが両方とも同じくらいの高さ

V6は下がって上がって下がり、下がるのはほんのちょっと

というのがまあいい感じでできた感じか…

カージオイドの形や傾き具合が調整できたらもっといい感じになりそう。

2011-02-22

MikuHatsune2011-02-22

ん?

アニメを観ていたら、病院シーンが出てきた。

循環器終わったからわかる。

なんだこの心電図は…?

いや待て。これが心電図でないとしよう。

しかしこの数秒後、よくある「ピー」の音とともに患者さんは亡くなってしまった。

これは心電図なのか…?

torsade de pointes…?

左上の180(130?)/60も気になるぞ…?

f:id:MikuHatsune:20110221233953j:image

2011-02-02

MikuHatsune2011-02-02

そんなことより豆まきしようぜ!!

循環器内科・心臓血管外科学を学んでいる。

8割方動脈硬化の話だ。MIだ。

そういえば生理学で、血管内の血流はずり応力による層流になるとかならんとか言ってた気がする…。

NKの拡散の話を聞いて、まあなんか簡単に書いてみた。

 

x,y,z空間で、最初はz=0(時間t=0)から始まる。

xy平面において、血管腔は半径R。

血管は剛体。形は変わらない。

血液を限りなく微小にし、質量1の質点的なものとする。

あるzでのxy平面上の血液は、自分の存在するマスの周囲8マスから、隣接する血液の速度に応じた力を受ける。

つまり、自分の隣の血液が1の速度で動いていたら、fric倍の力を受ける。これを8マス分合算。

速度と書いたけど、xy平面に垂直に移動しているものとする。

高校物理を思い出すと

F=ma

とかあったので、加速度から時間tにおける速度、変位を出す。

ニョキッとしたやつが描きたい。

# パラメータ
n<- 31                 # 考える世界の広さ。
R<- 5                  # 血管半径。
fric<- 0.0001          # 摩擦係数的な。
A<- B<-  diag(0,n)
Astart<- diag(1,n)
Layout<- matrix(c(2,1,1,1,3,
                  1,2,1,3,1,
                  1,1,1,1,1,
                  1,1,1,1,1,
                  1,1,1,1,1),5,5,byrow=T)
for(i in 1:n){
  for(j in 1:n){
    if(((i-ceiling(n/2))^2+(j-ceiling(n/2))^2)<=R^2){
      B[i,j]<- 1       # 0でないところが血管腔。
    }
  }
}
V<- diag(0,n+2)           # 端の影響のため、付け足す。
accel<- diag(0,n)         # 加速度。
pos<- diag(0,n)           # z軸方向への変位。
V[2:(n+1),2:(n+1)]<-  B
C<- B
Z<- 100
for(t in 1:Z){
  V[2:(n+1),2:(n+1)]<- C
  for(a in 2:(n+1)){
    for(b in 2:(n+1)){
      accel[a-1,b-1]<- fric*(V[a-1,b+1]+  V[a,b+1]+V[a+1,b+1]+
                             V[a-1,b  ]+0*V[a,b  ]+V[a+1,b  ]+
                             V[a-1,b-1]+  V[a,b-1]+V[a+1,b-1])
    }
  }
  C<- C+accel*t           # Cは速度だな。
  pos<- pos+C*t+(accel*t^2)/2
  pos[B==0]<- 0
  #file.name=paste("test",t,".png",sep="")
  #png(file=file.name)
    layout(Layout)
    persp(pos,col=2,phi=20,theta=60,zlim=c(0,80000),
          xlab="x",main="position")
    image(accel,main="accel",col=topo.colors(1000))  
    image(C,main="velocity",col=topo.colors(1000))
  #dev.off()
}

左上がz軸上方から眺めた加速度、右上が同じ視点の速度、perspがシミュレーション上の血流。

発展させる予定は今のところなし。

NKがセミナーで言っていたように、xyを一辺に計算できたらもっと計算速度があがるのになぁ。

グラフィックもぎこちないぞ。層っていうくらいだから同じ速度のところを色分けするのもいいと思う。

スクリプト書いてからwikiったが、もっと真面目にモデル化しておいたほうがきれいな図が描けそう。

 

そんなことより豆まきしようぜ!!