循環とは、酸素運搬と組織灌流

読んだ。

COI:評判よいらしいから買った。著者は全く知らない。

循環とは、よく見るDO2 の式の酸素運搬と、組織灌流\Delta P、すなわち脳還流などでよく見る脳血流=MAP - ICP のもの。
その他、筆者がものすごい好きなのだろうが、改定スターリングの法則とかグリコカリックスによる血管透過性などが語られていた。
この本で読むべきなのは、後半の血管作動薬、利尿薬の項で、こちらの方が実践的で明日からすぐ使える知識感が高かった。

間質性肺炎、完全に理解した

読んだ。

COI:なし

間質性肺炎急性増悪からICU入室が連続していたので読んだ。
間質性肺炎については分類が多すぎてほとんど理解していなかったが、とりあえず分かるようになった。
胸部CTの所見から市中肺炎、間質性肺炎など一般的な呼吸器内科診療が載っていたのでそれなりによかった。
肺癌については専門でないのか、記述は薄かった。

R shiny で出力される図のサイズを入力パラメータでうまくしたい

R shiny で図を出力するのだが、図の大きさをinputで入力したとき、renderPlotの引数(widthheight)にそれぞれ持たせようとすると、エラーになる。
状況としてはこんな感じである。

library(shiny)
server <- function(input, output, session) {
  height <- as.numeric(input$imageHeight)
  width  <- as.numeric(input$imageWidth)
  output$myplot <- renderPlot({
    plot(iris)
  }, height=height, width=width)
}

ui <- fluidPage(
  numericInput("imageWidth", "図の幅", 480, min = 1, max = 2000),
  numericInput("imageHeight", "図の幅", 560, min = 1, max = 2000),
  mainPanel(
    plotOutput(outputId="myplot", height="auto")
  )
)
shinyApp(ui, server)
 警告:  Error in : Can't access reactive value 'imageHeight' outside of reactive consumer.
ℹ Do you need to wrap inside reactive() or observe()?
  55: <Anonymous>
Error : Can't access reactive value 'imageHeight' outside of reactive consumer.
ℹ Do you need to wrap inside reactive() or observe()?

outputrenderPlotのオブジェクトを代入しようとしているが、renderPlot自体にはそれより先に作成されているはずのオブジェクト(height)はうまくなんやかんや出来ないらしい。
reactiveobserveを使え、と言われている。
というわけで探したらあった。
Problem getting height of renderPlot to be reactive to number of questions plotted - #6 by WillP - shiny - RStudio Community

library(shiny)
server <- function(input, output, session) {
  height <- reactive(input$imageHeight)   # reactive
  width  <- reactive(input$imageWidth)
  observe({                                     # observe
    output$myplot <- renderPlot({ 
      plot(iris)
    }, height=height(), width=width())    # () と関数化されている
  })
}

ui <- fluidPage(
  numericInput("imageWidth", "図の幅", 480, min = 1, max = 2000),
  numericInput("imageHeight", "図の幅", 560, min = 1, max = 2000),
  mainPanel(
    plotOutput(outputId="myplot", height="auto")
  )
)
shinyApp(ui, server)

出力されるhtml環境に合わせて図を拡大縮小しようとも思ったが、ブラウザのサイズ(というか図のサイズ)を定量的に保存しておきたいときに、ドラッグでブラウザのサイズを変えるのはよろしくないので採用できなかった。
Shiny でプロットの高さをブラウザ画面のサイズに合わせて変更する | Atusy's blog

そもそも図を描出するときに、Rだとdev.new(width=7, height=7)でプロットのためのキャンバスが登場するが、これをドラッグして拡大縮小したあとで、定量的にサイズ指定する方法がわからない。
拡大縮小したあとは、par()$dinでキャンバスの大きさを取得はできるが、par()$dinは読み取り専用で、任意の値の指定が出来ないので困っている。
https://www.rstudio.com/wp-content/uploads/2016/10/how-big-is-your-graph.pdf

日付を曜日に変換するのに日本語環境と英語環境で罠にハマった話

Rの日付に関する関数で、

format(Sys.Date(), "%a")

とすると、日本語環境ならば

[1] "月"

これをR Shiny で使おうと思っていたら、shiny のサーバーは英語環境なので、デフォルトでは日月火水木金土が出てこず、英語表記のSun Mon... になってしまって罠にハマった。
環境を英語にしようと

Sys.setenv(LANG="en")

にしてみるが、これではうまく行かない。

困っていたら弟子に教えてもらったが

temp <- Sys.getlocale(category="LC_TIME")
Sys.setlocale(category="LC_TIME", "C")
format(Sys.Date(), "%a")
Sys.setlocale(category="LC_TIME", temp)
[1] "C"
[1] "Mon"
[1] "ja_JP.UTF-8"

となり、英語化はうまくいった。

英語から日本語化は環境設定的にうまくいくのか確かめにくいので、面倒ではあるが、曜日番号を%u(月火水木金土日が1-7)もしくは%w(月火水木金土が1-6で、日が0)にして、月火水木金土日のオブジェクトを作って変換したほうが無難かもしれない。

COVID-19と無症候性低酸素血症

書いた。
Asymptomatic Hypoxemia as a Characteristic Symptom of Coronavirus Disease: A Narrative Review of Its Pathophysiology. COVID 2022, 2(1), 47-59.
COVID-19において、低酸素血症なのにぜんぜん苦しそうにない、いわゆるhappy hypoxia もしくはsilent hypoxia と言われている現象の総説である。
COVID-19罹患時に限らず、いわゆる低酸素応答とは、みたいなことも概説してあるので興味があれば。

査読で「日本の感染者数だけではなく世界について語って」みたいなことを言われたのでデータを探したら、Google のトップページでもこれを基に図を作っている、という元データがあった。
Coronavirus (COVID-19) Testing - Statistics and Research - Our World in Data
GitHub - owid/covid-19-data: Data on COVID-19 (coronavirus) cases, deaths, hospitalizations, tests • All countries • Updated daily by Our World in Data
国別のPCR陽性者数、検査数、死亡者数などいろいろあって、かつ国が大陸ごとにラベル付けられていたのでこれを使った。
f:id:MikuHatsune:20220109215619p:plain

library(xts)
url <- "https://covid.ourworldindata.org/data/owid-covid-data.csv"
dat <- read.csv(url, check.names=FALSE, stringsAsFactors=FALSE)

dat$date <-  as.Date(dat$date)
Month <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

y1 <- as.numeric(format(dat$date[1], "%Y"))
y1 <- as.Date(sprintf("%s-01-01", format(dat$date[1], "%Y")))
s <- as.numeric(mapply(format, dat$date[1], c("%Y", "%m"))

t0 <- c(
  c(mapply(function(z) tapply(z$new_cases, z$date, sum, na.rm=TRUE), split(dat, dat$continent))),
  Japan=mapply(function(x) tapply(x$new_cases, x$date, sum), list(subset(dat, location=="Japan")), SIMPLIFY=FALSE)
)

d <- as.Date(names(t0[[which.max(sapply(t0, length))]]))
t1 <- mapply(function(z) replace(rep(NA, length(d)), match(names(z), as.character(d)), z), t0)
rownames(t1) <- d
t1[, c(colnames(t1)=="")] <- rowSums(t1[, colnames(t1) != "" & colnames(t1) != "Japan"], na.rm=TRUE)
colnames(t1)[c(colnames(t1)=="")] <- "World"

t2 <- ts(t1, start=c(s[1], dat$date[1]-y1), frequency=365)
t7 <- apply(t1, 2, filter, rep(1, 7))/7
emg_state <- lapply(list(c("2020-04-07", "2020-05-25"),
                         c("2021-01-08", "2021-03-07"),
                         c("2021-04-25", "2021-06-20"),
                         c("2021-07-12", "2021-09-30")), as.Date)

m <- as.numeric(format(d, "%m"))
m0 <- rle(m)
m1 <- c(1, cumsum(m0$lengths))

Y <- as.numeric(format(d, "%Y"))
Y0 <- rle(Y)
Y1 <- c(1, cumsum(Y0$lengths))

# svg("figure.svg", 1200/96, 800/96)
cols <- c("black", "#0000FF", "#00FFFF", "green2", "yellow3", "orange", "#FF0000", "white")
names(cols) <- colnames(t2)
xl <- range(d) + c(1, -1) * 20

w7 <- 100000
par(mfrow=c(2, 1), mar=c(5, 6, 1, 1), las=1)
plot(d, t7[,1]/w7, type="n", lwd=10, las=1, xlab="", ylab="", cex.axis=1.2,
     xaxt="n", yaxt="n", xlim=xl, ylim=c(0, 10), frame=FALSE)
pa <- par()$usr
# d1 <- c(pa[1], d, pa[2])
d1 <- as.numeric(c(d, pa[2]))
x1 <- (head(d1[m1], -1) + tail(d1[m1], -1))/2
# axis(1, at=d1[m1[-1]], labels=NA)
axis(1, at=d1[m1], labels=NA)
ats <- seq(0, 10, 2); axis(2, at=ats, cex.axis=1.5)
abline(h=ats, lty=3)
for(i in seq(ncol(t7)-1)){
  lines(d, t7[,i]/w7, col=cols[i], lwd=7)
  lines(d, t2[,i]/w7, col=cols[i], lwd=2)
}
for(i in seq(x1)){axis(1, at=x1[i], labels=Month[m0$values][i], tick=FALSE, cex.axis=1.2)}
y1 <- c(pa[1], head(d[cumsum(Y0$lengths)], -1), pa[2])
xd <- 5
yd <- pa[3]-1.75
for(i in 1:(length(y1)-1)){
  segments(y1[i]+xd, yd, y1[i+1]-xd, xpd=TRUE, lwd=2)
  axis(1, at=mean(y1[c(i, i+1)]), labels=Y0$values[i], line=2, tick=FALSE, cex.axis=1.2)
  #text(mean(y1[c(i, i+1)]), yd, Y0$values[i], adj=c(0.5, 2), xpd=TRUE)
}
axis(1, at=pa[1], labels="Month", las=1, cex.axis=1, tick=FALSE, line=0, hadj=1, cex.axis=1.2)
axis(1, at=pa[1], labels="Year", las=1, cex.axis=1, tick=FALSE, line=2, hadj=1, cex.axis=1.2)
legend("topleft", legend=head(names(cols), -1), pch=15, col=cols, cex=1.5)
box()
#mtext("Time", 1, line=5, cex=1.5)
mtext(substitute("Positive cases ["*10^x*"]", list(x=log10(w7))), 2, line=3, cex=2, las=3)

w2 <- 10000
plot(d, t2[,"Japan"]/w2, type="n", lwd=5, las=1, xlab="", ylab="", cex.axis=1.2, xaxt="n", yaxt="n", xlim=xl, ylim=c(0, 3), frame=FALSE, cex.axis=1.5)
xi <- c(1, 2, 2, 1); yi <- c(3, 3, 4, 4)
pa <- par()$usr
d1 <- as.numeric(c(d, pa[2]))
axis(1, at=d1[m1], labels=NA)
ats <- seq(0, 3, 1); axis(2, at=ats, cex.axis=1.5)
abline(h=ats, lty=3)
for(i in seq(emg_state)){polygon(emg_state[[i]][xi], pa[yi], border=NA, col=grey(0.9))}
lines(d, t2[,"Japan"]/w2, col="black", lwd=2)
lines(d, t7[,"Japan"]/w2, col="black", lwd=7)
x1 <- (head(d1[m1], -1) + tail(d1[m1], -1))/2
for(i in seq(x1)){axis(1, at=x1[i], labels=Month[m0$values][i], tick=FALSE, cex.axis=1.2)}
y1 <- c(pa[1], head(d[cumsum(Y0$lengths)], -1), pa[2])
xd <- 5
yd <- pa[3]-0.5
for(i in 1:(length(y1)-1)){
  segments(y1[i]+xd, yd, y1[i+1]-xd, xpd=TRUE, lwd=2)
  axis(1, at=mean(y1[c(i, i+1)]), labels=Y0$values[i], line=2, tick=FALSE, cex.axis=1.2)
  #text(mean(y1[c(i, i+1)]), yd, Y0$values[i], adj=c(0.5, 2), xpd=TRUE)
}
axis(1, at=pa[1], labels="Month", las=1, cex.axis=1, tick=FALSE, line=0, hadj=1, cex.axis=1.2)
axis(1, at=pa[1], labels="Year", las=1, cex.axis=1, tick=FALSE, line=2, hadj=1, cex.axis=1.2)
mtext(substitute("Positive cases ["*10^x*"]", list(x=log10(w2))), 2, line=3, cex=2, las=3)
box()
#  dev.off()

ランダムに昇進を決めることが全体の価値を最も高める、というPeterの法則とイグノーベル賞の話

読んだ。
[0907.0455] The Peter Principle Revisited: A Computational Study

こんな動画を見た。

www.youtube.com
2010年のイグノーベル賞で、Peterの法則というものをシミュレーションで示した、という話。

社員地位構造が階層的な場合、昇進するときに能力に応じて昇進させるより、ランダムに昇進させたほうが、全体としての価値が高まる、という結論である。

論文では、Peterの法則と一般常識的な仮説のふたつと、昇進させるときのルールみっつを想定して、各々かけあわせて6つの結果をシミュレーションしている。
Peterの法則は、各階層において社員が昇進すると、その新たな部署においてはいままでの能力は引き継がれず、再度(ランダムに)能力を開発する。例えば、SE業だけやっていた平社員が突然管理職になったときに、新たに必要となる管理職能力に対して、これまで培ってきたSE能力は一切加味されない。これはシミュレーションでは、能力値を昇進時に新たにランダムサンプリングすることで表現される。

一般常識的(Common sense strategy)な仮説では、昇進時にこれまでの能力は引き継がれる、という設定である。これは先の例でいくと、管理職業務にSE能力を流用できる、ということになる。これはシミュレーションでは、昇進前の能力値に対して、(個人ではなく、すべての社員が持ちうる理論上)最大能力値の±10%以内の能力\delta が昇進時に加算される、とする。この\delta の選び方は明示されていなかったので正規分布でも一様分布でもいいが、一様分布にした。

昇進のさせ方は3通りあり、ある役職が空いたときに、下の役職の社員たちの中から、能力値が(1)最大の者、(2)最低の者、(3)に関係なくランダム、を選ぶという作戦である。

シミュレーションの細かい設定としては、全社員n=160nの会社において、役職が6つある。i=\{1,2,3,4,5,6\}である。各役職の人数はn_i niであり、n_i > n_{i+1}, n=\sum_i n_i である。
社員の能力はcompetenceというある実数で定義される。これはいろいろな能力が合算された結果、なんらかのひとつの値になったものとする。社員は年齢ageが定義されており、competence正規分布N(7, 2) かつ[1, 10]の範囲、age正規分布N(25, 5)かつ[18, 60] の範囲である。両端が切られているのでなんとなく切断正規分布truncnorm パッケージを使用した。各社員agentには、competenceage、役職の階層levelの属性がある。

さて、昇進させるには役職が空かないといけないが、役職を空けるには解雇しなければならない。解雇のルールは、competenceについて、恣意的に4以下の場合は、次の時刻ステップで解雇される。また、各時刻ステップでageは1つ加算され、60歳以上は解雇される、とする。役職が空いたときに、下の階層から、上の3つの昇進ルールで昇進させ、新たにcompetenceageをランダムに割り振ったlevelが1の新入社員を採用し、全体はn 人が保たれるようにする。

会社全体の価値は、各社員のcompetenceを、各役職の重要度r_i, (r_i < r_{i+1})riで重み付け和を取ったものとする。会社の価値global efficiency E は、ある階層でのcompetence の和C_i を用いて
E=10\frac{\sum_i C_i r_i}{\sum_i n_i r_i}
となる。これは[0, 100] の範囲になる。

これを1000時間単位、50回(こちらの環境では30回)行うと、figure 2 となる。
f:id:MikuHatsune:20211201235743p:plain

社員の能力は昇進先でも活かされるはずだ、という一般常識的な仮説が真であれば、昇進させる人物はcompetenceが高い者を選ぶと、会社全体の価値は上がる。一方で、昇進先ではこれまでの能力は活かされず、一新される(新たにランダムサンプリングした能力で働かなければならない)というPeterの法則が真であれば、competenceが低いものを昇進させるほうが、会社全体の価値は上がる。というより、後者においては、配置換えに近いイメージである。
一般常識的な仮説もしくはPeterの法則と、最大能力/最低能力者を昇進させる、という戦略があべこべの場合、会社全体の価値は下がる。

一般常識的な仮説が真なのか、Peterの法則が真なのか、はたまたその中間が真なのかは誰にもわからないが、ランダムに昇進させる社員を選べば、少なくとも最初の状態よりかは会社全体の価値は上がる、ということになる。

library(truncnorm)
ni <- c(81, 41, 21, 11, 5, 1)
n <- sum(ni)
ri <- c(0.2, 0.4, 0.6, 0.8, 0.9, 1)
delta <- 0.1

# competence parameter
comp_p <- c(7, 2)
# age parameter
age_p <- c(25, 5)
# range
trunc <- list(competence=c(1, 10), age=c(18, 60))

thres <- c(competence=4, age=60)


niter <- 30
Time <- 1000
promote_method <- c("The best", "The worst", "Random")
hypothesis <- c("Common", "Peter")
Eff <- array(0, c(niter, Time, length(promote_method), length(hypothesis)))
# the best: 1
# the worst: 2
# random: 3
# Peter: P
# Common: C
for(hy in seq(hypothesis)){
for(pm in seq(promote_method)){
  for(iter in seq(niter)){
    agent <- cbind.data.frame(
      competence=rtruncnorm(n, trunc$competence[1], trunc$competence[2], comp_p[1], comp_p[2]),
      age=floor(rtruncnorm(n, trunc$age[1], trunc$age[2], age_p[1], age_p[2])),
      level=unlist(mapply(rep, seq(ni), ni))
      )
    Eff[iter, 1, pm, hy] <- sum(agent$competence*ri[agent$level])/(sum(ni*ri))*10
    for(ti in seq(Time)[-1]){
      agent$age <- agent$age + 1
      # retirement
      ret_idx <- agent$competence < thres["competence"] | agent$age > thres["age"]
      ret_agent <- agent[ret_idx,]
      ret_levels <- unique(ret_agent$level)
      n_ret <- c(table(factor(ret_agent$level, factor(seq(ni)))), 0)
      nl <- sum(ret_idx)
      if(nl > 0){ # no retirement
        new_agent <- cbind.data.frame(
          competence=rtruncnorm(nl, trunc$competence[1], trunc$competence[2], comp_p[1], comp_p[2]),
          age=floor(rtruncnorm(nl, trunc$age[1], trunc$age[2], age_p[1], age_p[2])),
          level=1
        )
        # promotion
        if(promote_method[pm] == "Random"){ # random
          promote_idx <- unlist(mapply(function(l) sample(which(agent$level==l-1 & !ret_idx), sum(n_ret[l:length(n_ret)])), 2:length(ni)))
        } else if (promote_method[pm] == "The best"){ # best
          promote_idx <- unlist(
            mapply(function(l){
              a <- which(agent$level==l-1 & !ret_idx)
              tail(a[order(agent$competence[a])], sum(n_ret[l:length(n_ret)]))
            }, 2:length(ni))
          )
        } else { # worst
          promote_idx <- unlist(
            mapply(function(l){
              a <- which(agent$level==l-1 & !ret_idx)
              head(a[order(agent$competence[a])], sum(n_ret[l:length(n_ret)]))
            }, 2:length(ni))
          )
        }
        agent[ret_idx,] <- new_agent
        agent$level[promote_idx] <- agent$level[promote_idx] + 1
        if(hypothesis[hy] == "Peter"){
          # Peter hypothesis
          agent$competence[promote_idx] <- rtruncnorm(length(promote_idx), trunc$competence[1], trunc$competence[2], comp_p[1], comp_p[2])
        } else {
          # Common sense hypothesis
          d <- trunc$competence[2]*delta
          inherits <- runif(length(promote_idx), min=-d, max=d)
          agent$competence[promote_idx] <- agent$competence[promote_idx]+inherits
        }
      }
    Eff[iter, ti, pm, hy] <- sum(agent$competence*ri[agent$level])/(sum(ni*ri))*10
    }
  }
  # boxplot(Eff[, , pm, hy], main=sprintf("%s hypothesis, %s strategy", hypothesis[hy], promote_method[pm]))
}
}

yl <- range(Eff)
xl <- c(1, Time)
h <- mean(Eff[,1,,])
cols <- matrix(c("green", "skyblue", "cyan", "red", "orange", "yellow3"), length(promote_method), length(hypothesis))
txt <- outer(promote_method, hypothesis, paste, sep="; ")
par(mar=c(4, 4.5, 2, 10), cex=2)
plot(0, xlim=xl, ylim=yl, type="n", xlab="Time", ylab="Efficiency [%]", las=1)
for(hy in seq(hypothesis)){
for(pm in seq(promote_method)){
  v <- mapply(function(z) z$conf.int, apply(Eff[,,pm,hy], 2, t.test))
  for(i in seq(Time)){
    segments(i, v[1, i], y1=v[2, i], col=cols[pm, hy])
  }
  text(par()$usr[2], mean(v[, Time]), txt[pm, hy], xpd=TRUE, pos=4, font=2, col=cols[pm, hy]
)
}}
abline(h=h, lty=3, lwd=3)
text(par()$usr[2], h, "Averaged Initial Efficiency", xpd=TRUE, pos=4, font=2, cex=0.7)

クレアチニンクリアランス

薬剤の投与量を考えるときに、クレアチニンリアランスというものを考慮しないといけないのだが、クレアチニンリアランスにはクレアチニンと男女と体重が必要である。
そんなものは電子カルテで自動計算してくれればいいのだが、デフォルトではしてくれない仕様である。
クレアチニンリアランスで、だいたい50以下、30以下、10以下、くらいがわかればいいのでとりあえずやってみる。
クレアチニンリアランスは、クレアチニンC、年齢a、体重b として
f(C)=\frac{(140-a)b}{72C}
で、女性は0.85倍される。

体重60kgとすると、<30 となる領域は80歳以上でCre >1.5 であるので、そんなに高齢でなければCre 2 くらいまではなんにも腎調節のことは考えずに初回投与してもたぶん大丈夫そう。
f:id:MikuHatsune:20211127180554p:plain

age <- 16:100
cre <- seq(0, 10, by=0.1)[-1]
bw <- 60

d <- c(0, 10, 30, 50, 10000)
CCR <- mapply(function(z) matrix(c(cut(outer((140-age)*bw/72, cre, "/")*z, d)), length(age), length(cre)), c("male"=1, "female"=0.85), SIMPLIFY=FALSE)

v0 <- seq(0, 100, by=5)
v1 <- seq(0, 10, by=0.5)
txt <- sprintf("%d - %s", head(d, -1), c(d[-c(1, length(d))], ""))
cols <- rev(c("lightgreen", "yellow", "orange", "red"))
par(mfrow=c(1, 2), mar=c(5, 5, 2, 2), las=1, cex.lab=1.5, cex.axis=2)
for(i in seq(CCR)){
  image(age, cre, CCR[[i]], xlab="Age", ylab="Creatinine clearance [mg/dL]", col=cols, main=names(CCR)[[i]])
  abline(v=v0, h=v1, lty=3)
  legend("topright", legend=txt, bg="white", title="CCr", col=cols, pch=15, cex=1.5)
}