読んだ。
COI：評判よいらしいから買った。著者は全く知らない。
循環とは、よく見るDO2 の式の酸素運搬と、組織灌流、すなわち脳還流などでよく見る脳血流=MAP - ICP のもの。
その他、筆者がものすごい好きなのだろうが、改定スターリングの法則とかグリコカリックスによる血管透過性などが語られていた。
この本で読むべきなのは、後半の血管作動薬、利尿薬の項で、こちらの方が実践的で明日からすぐ使える知識感が高かった。
R shiny で出力される図のサイズを入力パラメータでうまくしたい
R shiny で図を出力するのだが、図の大きさをinputで入力したとき、renderPlotの引数（widthとheight）にそれぞれ持たせようとすると、エラーになる。
状況としてはこんな感じである。
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()?
outputにrenderPlotのオブジェクトを代入しようとしているが、renderPlot自体にはそれより先に作成されているはずのオブジェクト（height）はうまくなんやかんや出来ないらしい。
reactiveかobserveを使え、と言われている。
というわけで探したらあった。
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陽性者数、検査数、死亡者数などいろいろあって、かつ国が大陸ごとにラベル付けられていたのでこれを使った。
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%以内の能力 が昇進時に加算される、とする。この
の選び方は明示されていなかったので正規分布でも一様分布でもいいが、一様分布にした。
昇進のさせ方は3通りあり、ある役職が空いたときに、下の役職の社員たちの中から、能力値が（1）最大の者、（2）最低の者、（3）に関係なくランダム、を選ぶという作戦である。
シミュレーションの細かい設定としては、全社員人
nの会社において、役職が6つある。である。各役職の人数は
niであり、,
である。
社員の能力はcompetenceというある実数で定義される。これはいろいろな能力が合算された結果、なんらかのひとつの値になったものとする。社員は年齢ageが定義されており、competenceは正規分布 かつ
の範囲、
ageは正規分布かつ
の範囲である。両端が切られているのでなんとなく切断正規分布
truncnorm パッケージを使用した。各社員agentには、competence、age、役職の階層levelの属性がある。
さて、昇進させるには役職が空かないといけないが、役職を空けるには解雇しなければならない。解雇のルールは、competenceについて、恣意的に4以下の場合は、次の時刻ステップで解雇される。また、各時刻ステップでageは1つ加算され、60歳以上は解雇される、とする。役職が空いたときに、下の階層から、上の3つの昇進ルールで昇進させ、新たにcompetenceとageをランダムに割り振ったlevelが1の新入社員を採用し、全体は 人が保たれるようにする。
会社全体の価値は、各社員のcompetenceを、各役職の重要度, (
)
riで重み付け和を取ったものとする。会社の価値global efficiency は、ある階層での
competence の和 を用いて
となる。これは の範囲になる。
これを1000時間単位、50回（こちらの環境では30回）行うと、figure 2 となる。
社員の能力は昇進先でも活かされるはずだ、という一般常識的な仮説が真であれば、昇進させる人物は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以下、くらいがわかればいいのでとりあえずやってみる。
クレアチニンクリアランスは、クレアチニン値、年齢
、体重
として
で、女性は0.85倍される。
体重60kgとすると、<30 となる領域は80歳以上でCre >1.5 であるので、そんなに高齢でなければCre 2 くらいまではなんにも腎調節のことは考えずに初回投与してもたぶん大丈夫そう。
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) }