読んだ。
BeyondER (2巻1号) | 医書.jp
ヘリには乗らないのでドクターヘリやドクターカーの運用やエビデンスを勉強した。
院外心停止の蘇生も復習がてら読んだ。
ショック不応のアドレナリン投与までの時間とROSC確率で、スプラインを使った解析をしているので投与10分までは生存確率が横ばい、ないしは推定中央値でみるとやや上昇しているようにも見えなくもないが、普通に考えて投与までに時間がかかると生存率は低下するはずである。
となれば時間がかかると低下する、という制約をいれた解析をするべきである。
Time to Epinephrine Administration and Survival From Nonshockable Out-of-Hospital Cardiac Arrest Among Children and Adults
Circulation. 2018 May 8;137(19):2032-2040. のFigure 2 を使う。
https://pubmed.ncbi.nlm.nih.gov/29511001/ library(tiff) library(drc) a <- readTIFF("2032fig02.tif", as.is=TRUE) image(seq(nrow(a)), seq(ncol(a)), a) lab <- c(mapply(function(z) c(paste(head(z, -1), tail(z, -1), sep="-"), paste0("<", tail(z, 1))), list(c(1, seq(2, 30, by=2))))) xy$x <- rep(312, 16) points(xy) y_0 <- 312 # y=0 のピクセル x_0 <- 88 # x=0 の軸の値があるピクセル xy <- list(x=rep(312, 16), y=c(107, 145, 187, 222, 263, 302, 338, 379, 417, 457, 493, 536, 570, 614, 650, 688) ) Y <- y_0-mapply(function(i) tail(which(a[1:y_0, i] == 255), 1), xy$y) P <- mapply(function(i) range(which(10 < a[1:y_0, i] & a[1:y_0 ,i] < 51)), xy$y) v <- which(100 < a[1:y_0, x_0] & a[1:y_0, x_0] < 230) r <- rle(diff(v)) r0 <- r$lengths[r$values < 2] + 1 r1 <- tapply(head(v, -1), unlist(mapply(rep, seq(r0), r0)), c) p_unit <- 5/(y_0 - mean(r1[[1]]))/100 n_unit <- 6000/(y_0 - mean(r1[[1]])) p <- (y_0 - colMeans(P))*p_unit N <- floor(Y*n_unit) pN <- ceiling(Y*n_unit*p) library(rstan) library(vioplot) rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) code <- " data { int N; int Y[N]; int D[N]; } parameters{ ordered[N] theta; real<lower=0> s; } model{ for(i in 3:N){ theta[i] ~ normal(2*theta[i-1] - theta[i-2], s); } for(i in 1:N){ D[i] ~ binomial(Y[i], inv_logit(theta[i])); } } " m0 <- stan_model(model_code=code) standata <- list(N=length(N), Y=N, D=pN) standata <- lapply(standata, rev) fit <- sampling(m0, standata, iter=10000, warmup=5000, chain=40) ex <- extract(fit, pars=head(fit@model_pars, -1)) esp <- 1/(1+exp(-ex$theta[, fit@par_dims$theta:1])) vioplot(esp) code0 <- " data { int N; int Y[N]; int D[N]; } parameters{ ordered[N] theta; } model{ //for(i in 1:N){ // D[i] ~ binomial(Y[i], theta[i]); //} D ~ binomial(Y, inv_logit(theta)); } " code1 <- " data { int N; int Y[N]; int D[N]; } parameters{ real<lower=0, upper=1> d; real<lower=0> e; } model{ for(i in 1:N){ D[i] ~ binomial(Y[i], d*exp(-i/(1/e))); } } " ms <- mapply(stan_model, model_code=list(code0, code1)) standata <- list(N=length(N), Y=N, D=pN) fit <- mapply(sampling, ms, list(lapply(standata, rev), standata), iter=1500, warmup=1000, chain=16) ex <- mapply(extract, fit) esp <- list( 1/(1+exp(-ex[[1]]$theta[, fit[[1]]@par_dims$theta:1]))*100, c(ex[[2]]$d)*exp(outer(ex[[2]]$e, -seq(standata$N), "*"))*100 ) xl <- range(seq(lab)) yl <- c(0, max(mapply(max, esp))) cols <- c("Monotonic decrease"="skyblue", "Exponential decrease"="green2", "EXD.2 model"="yellow3", "Raw data"="black") par(mar=c(5, 5, 1, 1), las=1, cex.lab=2) plot(0, type="n", xlim=xl, ylim=yl, xlab="Minutes to Epinephrine Administration", ylab="Percent Survival", xaxt="n") abline(v=1:length(lab), lty=3) legend("topright", legend=names(cols), col="black", pt.bg=cols, pch=22, cex=1.5) axis(1, at=1:length(lab), labels=NA) for(i in seq(esp)){ vioplot(esp[[i]], add=TRUE, frame.plot=FALSE, side=c("left", "right")[i], col=cols[i]) } # exponential decrease の推定線 xd <- seq(0, length(lab)+2, length=3000) yd <- mean(ex[[2]]$d) * exp(-xd*mean(ex[[2]]$e))*100 lines(xd, yd, col=cols[2], lwd=3) lines(xd, yd, lwd=1, lty=3) # Monotonic decrease の推定線 xd <- seq(lab) yd <- colMeans(esp[[1]]) lines(xd, yd, col=cols[1], lwd=3) lines(xd, yd, lwd=1, lty=3) # drc を使ったEXD モデル drc_model <- drm(p ~ xd, fct=EXD.2(c(NA, NA))) xd <- seq(0, length(lab)+2, length=3000) yd <- predict(drc_model, as.data.frame(xd))*100 lines(xd, yd, col=cols[3], lwd=3) lines(xd, yd, lwd=1, lty=3) # 実のデータ mapply(text, 1:length(lab), par()$usr[3]-0.6, lab, srt=30, xpd=TRUE) lines(p*100, lwd=3, col=cols[4])
