2023-09-27

臨床外科：術後急変

麻酔外科学

読んだ。

臨床外科 2023年 7月号特集　術後急変！　予知・早期発見のベストプラクティス

医学書院

Amazon

臨床外科 (78巻7号) | 医書.jp
だいたい数日後にくる。

2023-09-21

プロット上のテキストの文字を一部分だけ色を変えたい

R Rを使いこなす

ABCDみたいなことをしたいが、文字をそのままtext関数で配置すると、一色しか使えない。
というわけで一文字ずつ着色するのだが、そうすると文字の配置がずれる。
なのでsubstituteにphantomを使って文字と同じ幅の空白を作って無理やり文字を配置する。

txt <- "RAINBOW"
w <- strsplit(txt, "")[[1]]
cols <- rainbow(length(w))
# 地道に書いて配置してみる
cex <- 2
par(mar=c(1, 1, 1, 1))
plot(1, type="n", frame=FALSE, xaxt="n", yaxt="n", xlab="", ylab="")
pa <- par()$usr
x <- mean(pa[1:2])
y <- seq(pa[4], pa[3], length=length(w)+2)
text(x, y[1], txt, cex=cex, xpd=TRUE)
text(x, y[2], substitute(phantom("")*"R"*phantom("AINBOW")), cex=cex, col=cols[1])
text(x, y[3], substitute(phantom("R")*"A"*phantom("INBOW")), cex=cex, col=cols[2])
text(x, y[4], substitute(phantom("RA")*"I"*phantom("NBOW")), cex=cex, col=cols[3])
text(x, y[5], substitute(phantom("RAI")*"N"*phantom("BOW")), cex=cex, col=cols[4])
text(x, y[6], substitute(phantom("RAIN")*"B"*phantom("OW")), cex=cex, col=cols[5])
text(x, y[7], substitute(phantom("RAINB")*"O"*phantom("W")), cex=cex, col=cols[6])
text(x, y[8], substitute(phantom("RAINBO")*"W"*phantom("")), cex=cex, col=cols[7])

# いっぺんに
for(i in seq(w)){
  w0 <- c("", w, "")
  hoge <- substitute(phantom(w1)*w2*phantom(w3), list(w1=paste0(w0[1:i], collapse=""), w2=w0[i+1], w3=paste0(w0[(i+2):(length(w)+2)], collapse="")))
  text(x, y[9], hoge, cex=cex, col=cols[i], xpd=TRUE)
}

2023-09-14

せん妄関係で読んだ雑誌

医学麻酔高度情報化

読んだ。

精神医学 2020年 4月号特集　ベンゾジアゼピン受容体作動薬の問題点と適正使用

医学書院

Amazon

精神医学 (62巻4号) | 医書.jp
ベンゾジアゼピンがなぜだめか、どう適正使用するべきかいろいろ書いてあった。

不眠に対してベンゾジアゼピン受容体作動薬を投与することは，精神薬理学的にはアルコールを服用させて入眠をさせているのと大きな違いがない

というのには最初は笑ってしまったが、GABA受容体と作用機序を考えると確かにそうだと納得した。

精神医学 2018年 3月号特集　せん妄をめぐる最近の動向

医学書院

Amazon

精神医学 (60巻3号) | 医書.jp
せん妄の対応について、一般病棟レベルの対応の話が多いがICUでも参考になる。
上記のベンゾジアゼピンのことを勉強するとようやくアルコール離脱せん妄でセルシンを内服させるのが治療になる機序が理解出来た。

精神医学 2012年 03月号オピニオン認知症の終末期医療の対応:現状と課題−尊厳をどう守るか

医学書院

Amazon

精神医学 (54巻3号) | 医書.jp
重症認知症は終末期ということが理解されていない。

2023-09-10

帰無仮説検定が正しいのか正しくないのか

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

検定とp値がいつまで経ってもわからないので久々にシミュレーションした。
対照群 $C$ と治療群 $T$ があって、なんらかの指標となる値に変化があるかどうか調べたい。
比較するのは、母集団（真の分布）である（ここでは真の分布は正規分布であると自分は分かっているので、 $N(\mu, \sigma)$ は自分で決められるしシミュレーションでは正確に比較検討できる）の平均 $\mu$ である。
ここではシミュレーションなので $\mu$ は分かっているが、実際のデータ解析においては母集団（真の分布）は分布の形もわからないしもちろん真のパラメータの値もわからないのが、実験して標本（サンプル）は分かる。この標本のデータを使って
帰無仮説 $H_0:\mu_C=\mu_T$ と対立仮説 $\mu_C\not=\mu_T$ を考えるのが一般的な帰無仮説検定とp値の扱いである。
$H_0$ が真と仮定してp値を計算したらp<0.05 なのでこれは棄却（本当は有意水準 $\alpha$ ）．．．というのがよくある考え方である。
単純に総当りとして $H_0$ が真か偽で、それぞれp>0.05 で棄却しないもしくはp<0.05 で棄却するパターンの2x2 通りをシミュレーションした。
B：帰無仮説 $H_0$ が真（本当に差がない）のにも関わらず、p<0.05 となり棄却する、というのは、p値の定義からは確率論的にありうるので、統計学的には正しい。しかし実践医学的には正しくはない。
C：帰無仮説 $H_0$ が偽（差があることが神視点では分かっている）なのに、p<0.05 にならずに帰無仮説が棄却されない、というのも、統計学的にはありえる正しい結果である。これは検出力による。だがこれも実践医学的には正しくない、というか好ましくない。

n <- 25 # サンプルサイズ
g <- rep(1:2, each=n)
alpha <- 0.05
dist_q <- 0.01 # プロット用の真の分布の上下限の分位点%

# 真の分布を正規分布から
ms <- c(10, 15)
sds <- 3

# シミュレーションデータをつくる
set.seed(1234)
flag <- TRUE
while(flag){
  y <- list(
    mapply(rnorm, n, ms[c(1, 1)], sds, SIMPLIFY=FALSE),
    mapply(rnorm, n, ms[c(1, 1)], sds, SIMPLIFY=FALSE),
    mapply(rnorm, n, ms[c(1, 2)], sds, SIMPLIFY=FALSE),
    mapply(rnorm, n, ms[c(1, 2)], sds, SIMPLIFY=FALSE)
  )
  pv <- mapply(function(z) t.test(z ~ g)$p.value, lapply(y, unlist))
  flag <- !(pv[1] > alpha & pv[2] < alpha & pv[3] > alpha & pv[4] < alpha)
}
for(i in seq(y)) names(y[[i]]) <- c("C", "T")
# 一気に見つけるにはよいが

# 4パターン同時に見つけるには対立仮説が真のときに帰無仮説が棄却されないパターンがなかなか見つからなくなるので個別にシミュレーションデータを作る
flag <- TRUE
while(flag){
  tmp0 <- mapply(rnorm, n, ms[c(1, 1)], sds, SIMPLIFY=FALSE)
  p0 <- t.test(unlist(tmp0) ~ g)$p.value
  flag <- !(p0 > alpha)
}
flag <- TRUE
while(flag){
  tmp1 <- mapply(rnorm, n, ms[c(1, 1)], sds, SIMPLIFY=FALSE)
  p1 <- t.test(unlist(tmp1) ~ g)$p.value
  flag <- !(p1 < alpha)
}
flag <- TRUE
while(flag){
  tmp2 <- mapply(rnorm, n, ms[c(1, 2)], sds, SIMPLIFY=FALSE)
  p2 <- t.test(unlist(tmp2) ~ g)$p.value
  flag <- !(p2 > alpha)
}
flag <- TRUE
while(flag){
  tmp3 <- mapply(rnorm, n, ms[c(1, 2)], sds, SIMPLIFY=FALSE)
  p3 <- t.test(unlist(tmp3) ~ g)$p.value
  flag <- !(p3 < alpha)
}
y <- list(tmp0, tmp1, tmp2, tmp3) 
for(i in seq(y)) names(y[[i]]) <- c("C", "T")
pv <- c(p0, p1, p2, p3)



cols <- c("orange", "green2")
colsa <- c(rgb(255, 165, 0, 50, maxColorValue=256), rgb(0, 165, 0, 125, maxColorValue=256))
cols <- gsub(".{2}$", "", colsa)
dx <- 0.3
xl <- c(-1, 3)
yl <- range(unlist(y))
yl <- range(mapply(function(z) mapply(qnorm, z, ms, sds), c(dist_q, 1-dist_q)), unlist(y))
x <- seq(yl[1], yl[2], length=3000)

# 真の分布のプロット用
Y <- list(
  mapply(dnorm, list(x), ms[c(1, 1)], sds, SIMPLIFY=FALSE),
  mapply(dnorm, list(x), ms[c(1, 1)], sds, SIMPLIFY=FALSE),
  mapply(dnorm, list(x), ms[c(1, 2)], sds, SIMPLIFY=FALSE),
  mapply(dnorm, list(x), ms[c(1, 2)], sds, SIMPLIFY=FALSE)
)
Y <- rapply(Y, function(z) -z*1/max(unlist(Y)), how="replace")
dY <- -10
# Y <- rapply(Y, function(z) z*dY, how="replace")


dx <- 0.5 # 帰無仮説が真のときにプロットを少しずらす
marginY <- 7
marginy <- 2
text <- c("差がなくて，有意にならない",
          "差がないのに，有意になる",
          "差があるのに，有意にならない",
          "差があって，有意になる")
text <- paste(LETTERS[seq(text)], text, sep=": ")
bcols <- c("green3", "red2")[c(1, 2, 2, 1)]

# png("significant.png", 640, 640)
par(mfrow=c(2, 2), mar=c(3, 3, 2, 1))
for(j in seq(y)){
  plot(0, type="n", xlim=xl, ylim=yl+c(0, marginY), frame=FALSE, xaxt="n", yaxt="n", xlab="", ylab="", main=text[j], cex.main=1.5)
  boxplot(y[[j]], at=1:2, add=TRUE, outline=FALSE, boxwex=0.5, col=colsa, frame=FALSE, xaxt="n", yaxt="n")
  axis(1, at=1:2, labels=names(y[[j]]), lwd=0, lwd.ticks=1, cex.axis=1.5)
  axis(1, at=(par()$usr[1]+0)/2, labels="真の分布", lwd.ticks=0, cex.axis=1.5)
  axis(2, lwd=0, lwd.ticks=1)
  box(col=bcols[j], lwd=3)
  stripchart(y[[j]], method="jitter", vertical=TRUE, add=TRUE, col=cols, pch=15, cex=1.5)
  for(i in 1:2) lines(Y[[j]][[i]], x+dx*(i-1), col=cols[i], lwd=5)
  abline(v=0, lty=3)
  yy <- sapply(y[[j]], max)
  # y2 <- (par()$usr[4] + max(yy))/2
  # y2 <- (par()$usr[4] + max(unlist(y)))/2
  y2 <- sum(c(par()$usr[4], max(unlist(y))) * c(1, 2)/3)
  for(i in 1:2) segments(i, y0=yy[i]+marginy, y1=y2, lwd=2)
  segments(x0=1, x1=2, y0=y2, lwd=2)
  if(pv[j] > alpha){
    txt <- substitute(atop("NS", "("*italic(p)~"="~x*")"), list(x=sprintf("%.3f", pv[j])))
  } else {
    txt <- substitute(atop(italic(p)~"<"~y, "("*italic(p)~"="~x*")"), list(x=sprintf("%.3f", pv[j]), y=alpha))
  }
  text(1.5, y2, txt, pos=3, cex=1.5)
}
# dev.off()

2023-09-09

みんなの呼吸器で読んだやつ

医学麻酔

みんなの呼吸器 Respica(レスピカ) 2022年5号(第20巻5号)特集:きほん～応用まで　人工呼吸器グラフィッククイズ50＋深おもしろい豆知識Q

メディカ出版

Amazon

いろいろなパターンが出ており勉強になる。

みんなの呼吸器 Respica(レスピカ) 2022年3号(第20巻3号)特集:症例で理解が深まる! “脱・初級"のための人工呼吸器設定とアセスメント

メディカ出版

Amazon

症例ベースで学ぶが大半が呼吸器非同調の原因探しである。
Oxygenation index (OI) というものがあるらしく OI=Paw * FiO2 / PaO2 で決まるようだが、平均気道内圧 Paw = PEEP * (PIP - PEEP) * Ti * RR/60 という人工呼吸器のグラフィックモニターの圧波形について動的コンプライアンス部分が1分間の呼吸に占める割合を、PF値で割れば求まるようで、下記のように評価するそうだ。