2024-01-01

日本版敗血症診療ガイドライン2020 (J-SSCG2020) の血糖管理でネットワークメタアナリシスをやってみる

R Rpackage Rを使いこなす医学統計麻酔

The Japanese Clinical Practice Guidelines for Management of Sepsis and Septic Shock 2020 (J-SSCG 2020)
「日本版敗血症診療ガイドライン2020 (J-SSCG2020)」正式版公開のお知らせ｜日本集中治療医学会
 日本版敗血症診療ガイドライン2020 (J-SSCG2020)」正式版公開のお知らせ
集中治療学会のほうは「許可なくリンクを貼ることを禁じます」と書いてあるが、ネットに無償公開しているものについて引用の範囲でリンクを貼ることについては法律的になんら問題はないので、APRINで倫理講習を受けましょう、とかアホなこと言ってないで法律的なことを真面目に勉強したほうがいいと思う。

ネットワークメタアナリシスというものがあり、ガイドライン中のいくつかのCQで検討が行われている。
CQ13-2 の「敗血症患者の目標血糖値はいくつにするか？」について、35件の研究が引用されている。
このうち、23番目のKalfon 2014 の論文では、付録の III >180 vs <110 アウトカム①短期死亡について、リスク比が Not estimated となっている。これは、<110 群で1336人にすべきところを136人と誤って解析しているからであるが、これを正しく1336人にして解析すると、論文と同じようなそれっぽい解析結果は出るが、対照 <110 群に対して144-180 群の結果が0.90 [0.83 - 0.97] となっているが、正しい入力データだと0.8946であり有効数字を考慮すると0.90 にならない。

Treatment estimate (sm = 'RR', comparison: other treatments vs '<110'):
            RR           95%-CI     z p-value
<110         .                .     .       .
>180    1.0102 [0.9480; 1.0766]  0.31  0.7538
110-144 0.8841 [0.7144; 1.0942] -1.13  0.2575
144-180 0.8946 [0.8272; 0.9676] -2.78  0.0054

Kalfon の誤った入力データでやってみるともちろん出力は異なるので、本解析と付録で入力データが異なるし、出力も完全再現に至らないので本当どういうことなんだろう（）

出力結果のグラフの点の大きさやエッジの太さを適当にいじったが、本文に書いていないので本当にあっているかは謎。
Control size of points in netgraph in R - Stack Overflow

library(NMA)
library(netmeta)

dat <- read.csv(text="
study,trial,trt,n,d,pubyear
1,Oksanen,110-144,51,18,2007
1,Oksanen,<110,39,13,2007
2,Finfer,144-180,3012,751,2009
2,Finfer,<110,3010,829,2009
3,Preiser,144-180,542,105,2009
3,Preiser,<110,536,125,2009
4,Green,144-180,36,9,2010
4,Green,<110,45,16,2010
5,Coester,144-180,40,4,2010
5,Coester,<110,39,2,2010
6,Cappi,144-180,35,10,2012
6,Cappi,<110,28,5,2012
7,van den Berghe,>180,765,55,2001
7,van den Berghe,<110,783,85,2001
8,Mitchell,>180,35,3,2006
8,Mitchell,<110,35,9,2006
9,van den Berghe,>180,605,242,2006
9,van den Berghe,<110,595,222,2006
10,Iapichino,>180,36,11,2008
10,Iapichino,<110,36,13,2008
11,De La Rosa Gdel,>180,250,96,2008
11,De La Rosa Gdel,<110,254,102,2008
12,Brunkhorst,>180,288,102,2008
12,Brunkhorst,<110,247,98,2008
13,Arabi,>180,257,83,2008
13,Arabi,<110,266,72,2008
14,Savioli,>180,45,13,2009
14,Savioli,<110,45,14,2009
15,Annane,>180,254,109,2010
15,Annane,<110,255,117,2010
16,Arabi,>180,120,45,2011
16,Arabi,<110,120,42,2011
17,Kalfon,>180,1312,447,2014
17,Kalfon,<110,1336,431,2014
18,Wang,>180,44,14,2017
18,Wang,<110,44,12,2017
19,McMullin,144-180,9,4,2007
19,McMullin,110-144,11,6,2007
20,Davies,>180,34,6,1991
20,Davies,110-144,35,6,1991
21,Walters,>180,12,0,2006
21,Walters,110-144,13,1,2006
22,Farah,>180,48,26,2007
22,Farah,110-144,41,19,2007
23,Bruno,>180,15,0,2008
23,Bruno,110-144,31,2,2008
24,Bilotta,>180,49,4,2008
24,Bilotta,110-144,48,4,2008
25,Chan,>180,55,3,2009
25,Chan,110-144,54,2,2009
26,de Azevedo,>180,169,42,2010
26,de Azevedo,110-144,168,38,2010
27,Hsu,>180,57,28,2012
27,Hsu,110-144,55,18,2012
28,Giakoumidakis,>180,107,7,2013
28,Giakoumidakis,144-180,105,1,2013", stringsAsFactors=FALSE)

hf <- setup(study=study, trt=trt, d=d, n=n, measure="RD", ref="<110", data=dat)
n <- nma(hf, eform=TRUE)

nmadat <- as.data.frame(t(sapply(split(dat[,-1], dat$study), unlist))[,c(-1, -10)])
for(i in 4:7){
  nmadat[,i] <- as.numeric(nmadat[,i])
}

p1 <- pairwise(treat=list(trt1, trt2), event=list(d1, d2), n=list(n1, n2), dat=nmadat, sm="RR")
net1 <- netmeta(p1, common = FALSE)
net1
cols <- c("skyblue3", "orange", "red", "green")
netgraph(net1, points=TRUE, cex=1.5, number.of.studies=TRUE, thickness="se.random", col.points=cols, cex.points=table(c(net1$treat1, net1$treat2)), plastic=FALSE, rotate=-90)

2023-12-05

ICUで遭遇する免疫異常

医学麻酔

読んだ。

INTENSIVIST Vol.15 No.4 2023 (特集:ICUで遭遇する免疫異常））

メディカル・サイエンス・インターナショナル

Amazon

INTENSIVIST (15巻4号) | 医書.jp

免疫学的には重要なのだろうがIL-X とかサイトカインとか羅列されて免疫学を概説されてもぶっちゃけ臨床的にはほとんど役に立たない。分子機構の説明など同じようなことが別項で重複している。
○○病のこういうときには薬剤XXをこう使え、みたいな話がほしいのだが結局大量ステロイドかパルス。抗体製剤の使い方は専門的すぎて集中治療医ではなく膠原病内科医の領域になる。
AIDS、膠原病はほとんど診る機会がないので、遭遇してしまったときに参照できればいいと思う。
免疫チェックポイント阻害薬によるirAE は増えててICUで遭遇する機会は今後どんどん増えるだろうが、ステロイドを入れるか、心筋炎だったらVAECMOすることになる（したことはある）が、「ICIによるirAEは可逆的でありICIでそもそもの悪性腫瘍に対する予後延長が見込めるので諦めてはいけない」って書いてあったけど正気か?と思った。

勉強したければ読んでみてもいいが、必読かというと微妙。

2023-11-25

心電図トレーニング

循環器内科学生理学医学麻酔

読んだ。
わからなすぎて不安は自信に全くならなかった。

medicina(メディチーナ) 2022年 8月号特集不安を自信に変える心電図トレーニング~専門医のtipsを詰め込んだ50問

医学書院

Amazon

medicina (59巻9号) | 医書.jp

症例ベースでいえばこちらのほうがわかりやすくて記憶に残りやすかった。

medicina 2016年 4月号特集心電図を詠む心に残る24症例から

医学書院

Amazon

medicina (53巻5号) | 医書.jp

2023-11-11

胸部外科　周術期管理の重要ポイント

医学外科学麻酔

胸部外科 (76巻10号) | 医書.jp
胸部外科 (70巻8号) | 医書.jp
特に真新しいことは書いてなかった。

2023-11-08

Atlas of Critical Care Procedures

医学麻酔

読んだ。

Atlas of Critical Care Procedures (English Edition)

Springer

Amazon

胸腔ドレーンであんなことがあったので手技の確認がてら読んだ。
トロッカーよりペアンで鈍的に挿入するほうが合併症は少ないらしい。

献体を使った実際の手技の写真が多いので実際のページ数より遥かにボリュームは少ない。
知っていたと思っていたが知識が増えたこともあって、例えば

乳頭が第4-5肋間のことが多く、これより下に挿入すると横隔膜損傷のリスクが高くなる
Aラインモニターをフラッシュしてみて、2, 3個程度オシレーションするのが適正なダンピングである。1個しかないとオーバーダンピングでなまった波形になりdicrotic notch が消失し、それ以上だとアンダーダンピングで尖った波形になる
蘇生的開胸術を行ったときの心臓マッサージは、両手で包み込み拍出ができるよう揉むこと
REBOA を挿入するときの体表面ランドマークは、胸骨中部がzone Ⅰ。zone Ⅲ は腎動脈遠位部（L1-2）から総腸骨動脈分岐（L4）で、分岐は臍が目印

神経系で脳室ドレナージ、腰椎穿刺、硬膜外麻酔が書いてあったが、硬膜外麻酔は記述が薄すぎで、しかも某脊椎ブロックの項もあるのにこれは手技のことは一切書いてなかった。
四肢コンパートメント症候群、腹部コンパートメント症候群、VAC療法が異様にたくさん書いてあった。

一通り復習するにはよかった。

2023-11-02

検査と技術

医学循環器内科学麻酔

検査と技術 2018年 3月号増刊号特集感染症クイックリファレンス

医学書院

Amazon

検査と技術 (46巻3号) | 医書.jp
感染症の復習をしようと思ったが、感染症そのものの記載は薄いし抗生剤も薄いし、でも菌体自体の記載は検査法や微生物学的特性などがなんかマニアック。
コロナ前なのでコロナについては記載はない。

検査と技術 2019年 3月号増刊号特集エキスパートが教える心・血管エコー計測のノウハウ

Amazon

検査と技術 (47巻3号) | 医書.jp
これはむちゃくちゃよかった。
ICUエコーが流行っているが、集中治療医が行うエコーは結局のところレベルが低い（）
検査技師が行う心臓エコーは、ちゃんと体位がとれて良い条件で検査するのでよいのは当たり前だが、どのように描出するのが良いのか、体位とプローブを当てる箇所の図、そのときに見えるエコー画像がしつこいほどに毎回並んで掲載されているので、わかりやすい。
パルスや連続波ドップラーのサンプリングをどこに持っていくべきか、どう設定すべきか、正常値と異常値、など網羅されているので、よいと思った。

2023-11-01

中心極限定理を確かめる

R Rを使いこなす数学統計

どんな分布（注）からでもサンプリングされた標本の平均をさらに繰り返して標本の分布を考えたとき、その分布は正規分布になる。
というのでいろいろ分布を変えて試してみる。
Rでは接頭語rdpqに続いて分布名があるので、標準で実装されている関数をapropos関数で取得すると、19個あるようだった。

 [1] "beta"     "binom"    "cauchy"   "chisq"    "exp"      "f"        "gamma"    "geom"     "hyper"   
[10] "lnorm"    "logis"    "nbinom"   "norm"     "pois"     "signrank" "t"        "unif"     "weibull" 
[19] "wilcox"

19個の関数で必要なパラメータを適当に決め、30個サンプリングして平均値を求め、これを3000回繰り返すとする。
関数からサンプリングした標本の標本平均と標本分散はデータから決まるし、関数とパラメータから理論的な平均と分散も計算できるのでこれも計算する。
乱数発生関数に引数をいれるときに無理やりテキスト形式で

rnorm(n, mean, sd)

を実現するために、テキスト処理をしたあとにparse関数とeval関数で無理やり評価して関数を実行している。

ヒストグラムの設定が適当なせいなのか理論的な分布との重ね合わせが合ってないが、だいたいが平均値の分布は正規分布に従っている。
ここで、コーシー分布がおかしいことになっているが、コーシー分布は平均値が定義できない、となっているので、おかしいのが正しい。

library(stringr)
rdpq <- c("r", "d", "p", "q")
fs <- mapply(function(z) str_remove(apropos(sprintf("^%s", z)), sprintf("^%s", z)), rdpq)

f <- fs[[1]]
for(i in 1:(length(fs)-1)){
  for(j in (i+1):length(fs)){
    f <- intersect(f, intersect(fs[[i]], fs[[j]]))
  }
}

n <- 30 # ある分布からのサンプリング数
p <- ev <- replicate(length(f), vector("list"))
names(p) <- names(ev) <- f

p[[f[1]]] <- c(10, 20)    # beta
p[[f[2]]] <- c(20, 0.2)   # binom
p[[f[3]]] <- ""           # cauchy
p[[f[4]]] <- c(5)         # chisq
p[[f[5]]] <- c(2)         # exp
p[[f[6]]] <- c(4, 6)      # f
p[[f[7]]] <- c(0.5, 2)    # gamma
p[[f[8]]] <- c(0.3)       # geom
p[[f[9]]] <- c(5, 10, 12) # hyper
p[[f[10]]] <- c(0.3, 1.5) # lnorm
p[[f[11]]] <- c(10, 3.3)  # logis
p[[f[12]]] <- c(10, 0.25) # nbinom
p[[f[13]]] <- c(6, 3.5)   # norm
p[[f[14]]] <- c(5)        # pois
p[[f[15]]] <- c(10)       # signrank
p[[f[16]]] <- c(4)        # t
p[[f[17]]] <- c(2.8, 9.4) # unif
p[[f[18]]] <- c(4, 3)     # weibull
p[[f[19]]] <- c(4, 7)     # wilcox

niter <- 3000 # n回サンプリングするシミュレーションの試行回数
res <- replicate(niter,
colMeans(
mapply(function(z0, z1){
         eval(parse(text=sprintf("r%s(%d%s)", z0, n, paste0(ifelse(all(z1==""), "", ","), paste0(z1, collapse=",")))))
        }, f, p)
)
)

xq <- seq(0, 1, length=10000)
dens <- mapply(function(z0, z1){
         txt <- sprintf("q%s(%f%s)", z0, xq, paste0(ifelse(all(z1==""), "", ","), paste0(z1, collapse=",")))
         y <- mapply(function(w) eval(parse(text=w)), txt)
         d <- mapply(function(w) eval(parse(text=sprintf("d%s(%f%s)", z0, w, paste0(ifelse(all(z1==""), "", ","), paste0(z1, collapse=","))))), y)
         return(list(q=y, d=d))
        }, f, p, SIMPLIFY=FALSE)

ev[["beta"]] <- c(
  E=p[["beta"]][1]/(p[["beta"]][1]+p[["beta"]][2]),
  V=p[["beta"]][1]*p[["beta"]][2]/((p[["beta"]][1]+p[["beta"]][2])^2*(p[["beta"]][1]+p[["beta"]][2]+1))
  )

ev[["binom"]] <- c(
  E=prod(p[["binom"]]),
  V=p[["binom"]][1]*p[["binom"]][2]*(1-p[["binom"]][2])
  )

ev[["cauchy"]] <- c(
  E=NA,
  V=NA
  )

ev[["chisq"]] <- c(
  E=p[["chisq"]][1],
  V=2*p[["chisq"]][1]
  )

ev[["exp"]] <- c(
  E=1/p[["exp"]][1],
  V=(1/p[["exp"]][1])^2
  )

ev[["f"]] <- c(
  E=ifelse(p[["f"]][2] > 2, p[["f"]][2]/(p[["f"]][2]-2), NA),
  V=ifelse(p[["f"]][2] > 4, 2*p[["f"]][2]^2*(sum(p[["f"]])-2)/p[["f"]][1]/(p[["f"]][2]-2)^2/(p[["f"]][2]-4), NA)
  )

ev[["gamma"]] <- c(
  E=prod(p[["gamma"]]),
  V=p[["gamma"]][1]*p[["gamma"]][2]^2
  )

ev[["geom"]] <- c(
  E=1/p[["gamma"]][1],
  V=(1-p[["gamma"]][1])/(p[["gamma"]][1]^2)
  )

ev[["hyper"]] <- c(
  E=p[["hyper"]][3]*p[["hyper"]][1]/(p[["hyper"]][1]+p[["hyper"]][2]),
  V=p[["hyper"]][3]*p[["hyper"]][1]/(p[["hyper"]][1]+p[["hyper"]][2])*(1-p[["hyper"]][1]/(p[["hyper"]][1]+p[["hyper"]][2]))*(p[["hyper"]][1]+p[["hyper"]][2]-p[["hyper"]][3])/(p[["hyper"]][1]+p[["hyper"]][2]-1)
  )

ev[["lnorm"]] <- c(
  E=exp(p[["lnorm"]][1]+p[["lnorm"]][2]/2),
  V=exp(2*p[["lnorm"]][1]+p[["lnorm"]][2])*(exp(p[["lnorm"]][2])-1)
  )

ev[["logis"]] <- c(
  E=p[["logis"]][1],
  V=pi^2*p[["logis"]][2]^2/3
  )

ev[["nbinom"]] <- c(
  E=p[["nbinom"]][1]*(1-p[["nbinom"]][2])/p[["nbinom"]][2],
  V=p[["nbinom"]][1]*(1-p[["nbinom"]][2])/p[["nbinom"]][2]^2
  )

ev[["norm"]] <- c(
  E=p[["norm"]][1],
  V=p[["norm"]][2]
  )

ev[["pois"]] <- c(
  E=p[["pois"]][1],
  V=p[["pois"]][1]
  )

ev[["signrank"]] <- c(
  E=p[["signrank"]][1]*(p[["signrank"]][1]+1)/4,
  V=p[["signrank"]][1]*(p[["signrank"]][1]+1)*(2*p[["signrank"]][1]+1)/24
  )

ev[["t"]] <- c(
  E=ifelse(p[["t"]][1] > 1, 0, NA),
  V=ifelse(p[["t"]][1] > 2, 1/(1-2/p[["t"]][1]), NA)
  )

ev[["unif"]] <- c(
  E=mean(p[["unif"]]),
  V=diff(p[["unif"]])^2/12
  )

ev[["weibull"]] <- c(
  E=p[["weibull"]][2] * gamma(1 + 1/p[["weibull"]][1]),
  V=p[["weibull"]][2]^2 * (gamma(1 + 2/p[["weibull"]][1]) - (gamma(1 + 1/p[["weibull"]][1]))^2)
  )

ev[["wilcox"]] <- c(
  E=prod(p[["wilcox"]])/2,
  V=prod(p[["wilcox"]])*(sum(p[["wilcox"]])+1)/12
  )
ev <- lapply(ev, "*", c(1, 1/n))

# 分布からの推定と理論値の確認
m <- rowMeans(res)
v <- apply(res, 1, var)
a <- cbind(do.call(rbind, ev), m=m, v=v)

cols <- c("distribution"="orange", "theorem"="blue")
par(mfrow=c(length(f), 2), cex.main=2)
for(i in seq(f)){
  ml <- list(m[i], ev[[i]]["E"])
  vl <- list(v[i], ev[[i]]["V"])
  de <- mapply(function(z1, z2) mapply(qnorm, xq, z1, z2), ml, vl, SIMPLIFY=FALSE)
  y <- mapply(function(z0, z1, z2) mapply(dnorm, z0, z1, z2), de, ml, vl, SIMPLIFY=FALSE)
  yl <- c(0, max(unlist(y), na.rm=TRUE))
  hist(res[i,], probability=TRUE, xlab="mean", ylab="density", main=f[i], ylim=yl, cex.main=2)
  for(j in seq(y)){
    lines(de[[j]], y[[j]], col=cols[j], lwd=5)
  }
  legend("topright", legend=names(cols), col=cols, pch=15, cex=2)
  plot(dens[[i]]$q, dens[[i]]$d, type="l", lwd=5, xlab="", ylab="density", main=sprintf("%s distribution", f[i]))
}