Rのshiny をweb公開したときに日本語フォントが表示されないのがlibrary(showtext) で解決しなかった

ちょっとしたアプリケーションを作って、自分以外の、Rが使えない人でもデータを入力したら出力を得られるような仕様にしたい。
そこで、R のshiny を使ってweb にアクセスしたら使えるようにするわけだが、ローカル環境で構築したshiny アプリではプロット上の日本語フォントがきちんと表示されるのに、shinyapps.io にアップロードしてserver 環境下で使うと、日本語フォントが表示されない。
これはshiny server がlinux ベースで、日本語フォントを標準装備していないかららしいが、library(showtext) とすると日本語フォントが対応するらしい。
shinyapps.io - RjpWiki

しなかった。

下記を参考に、shiny server にR からコマンドを実行してフォントを読み込ませる。
リンクではsh ファイルをbash で実行しているが、system 関数でR側からコマンドを実行しているので、全部R側のスクリプトに納めておいた。
shinyapps.ioで任意の日本語フォントを使う - mana.bi

library(shiny)
#download.file("https://raw.githubusercontent.com/ltl-manabi/shinyapps.io_japanese_font/master/use_ipaex_font.sh", destfile = "use_ipaex_font.sh")

bashcode <- "
mkdir ~/tmp;
cd ~/tmp;
curl -O -L https://moji.or.jp/wp-content/ipafont/IPAexfont/IPAexfont00401.zip;
unzip IPAexfont00401.zip;
mkdir -p ~/.fonts/ipa;
cp ./IPAexfont00401/*.ttf ~/.fonts/ipa;
fc-cache -f ~/.fonts;
"
system(bashcode)

# Define UI for data upload app ----
ui <- fluidPage(
  # App title ----
  titlePanel("Uploading Files"),
  # Sidebar layout with input and output definitions ----
  sidebarLayout(
    # Sidebar panel for inputs ----
    sidebarPanel(
      # Input: Select a file ----
      fileInput("file1", "Choose CSV File",
                multiple = FALSE,
                accept = c("text/csv",
                         "text/comma-separated-values,text/plain",
                         ".csv")),
      # Horizontal line ----
      tags$hr(),
      # Input: Checkbox if file has header ----
      checkboxInput("header", "Header", TRUE),
      # Input: Select separator ----
      radioButtons("sep", "Separator",
                   choices = c(Comma = ",",
                               Semicolon = ";",
                               Tab = "\t"),
                   selected = ","),
      # Input: Select quotes ----
      radioButtons("quote", "Quote",
                   choices = c(None = "",
                               "Double Quote" = '"',
                               "Single Quote" = "'"),
                   selected = '"'),
      # Horizontal line ----
      tags$hr(),
      # Input: Select number of rows to display ----
      radioButtons("disp", "Display",
                   choices = c(Head = "head",
                               All = "all"),
                   selected = "head")
    ),
    # Main panel for displaying outputs ----
    mainPanel(
      # Output: Data file ----
      tableOutput("contents"),
      plotOutput("plot"),
    )
  )
)

# Define server logic to read selected file ----
server <- function(input, output) {
  hoge <- iris
  dat <- reactive({
    if(is.null(input$file1$datapath)){
      return(NULL)
    } else {
      req(input$file1)
      huga <- read.csv(input$file1$datapath)
      return(huga)
    }
  })
  output$contents <- renderTable({
    # input$file1 will be NULL initially. After the user selects
    # and uploads a file, head of that data file by default,
    # or all rows if selected, will be shown.
    req(input$file1)
    # when reading semicolon separated files,
    # having a comma separator causes `read.csv` to error
    tryCatch(
      {
        df <- read.csv(input$file1$datapath,
                 header = input$header,
                 sep = input$sep,
                 quote = input$quote)
      },
      error = function(e) {
        # return a safeError if a parsing error occurs
        stop(safeError(e))
      }
    )
    if(input$disp == "head") {
      return(head(df, input$is_holiday))
    }
    else {
      #return(input$is_holiday)
      return(df)
    }
  })
  output$plot <- renderPlot({
    plot(1, type="n")
    title("日本語")
    text(1, 1, "ほげほげ")
  })
}

# Create Shiny app ----
shinyApp(ui, server)

新型肺炎の変異株で若年者は本当に重症化しやすいのか

結論から言うと、0.うんぬん%以内の増加は40-50歳代の年齢層であるようだが、それよりも80歳以上の高齢者での死亡率増加が大きく、かつ感染者数が増えており重症のままベッドを占拠しておりながらもすぐに死亡してベッドを空ける、というわけではないので、医療逼迫感から「若年の重症化が~」というのが広まっている、のではないかと思う(素人の感想レベル

この1-2ヶ月くらいでイギリス型変異株が蔓延し、インド型の変異株も蔓延しつつあるので、よくTVやネットで耳にするのが
「40(20代も?)から50代までの若くて、しかも基礎疾患のない人が重症化しており、変異株で重症化リスクが増加しているおそれがある」
という言説である。
大阪で勤務している人のインタビューや市長がそんなことを言っているのでそうなのかもしれないが、少なくとも(大阪ではない場末の)地方勤務の自分の感触では、重症化しているのは高齢者もしくは肥満、糖尿病などの高リスク50代以上である。

変異株は「若い方が重症化しやすい」 小池都知事、若者の重症化や後遺症に警鐘 | ENCOUNT
こんな記事を見かけたが、記事いわく

従来型より感染力が高く、都内で感染が置き換わりつつある変異株「N501Y」について、小池知事は「若い方が重症化しやすい。都の感染者は20代~30代が半数以上を占めており、重傷者数も20代~50代で倍増している。ECMO(人工心肺装置)は都内に7台あるが、現在その7台すべてが60代以下の方に使われている。だからこそ若い方にも注意していただきたい」と警戒を呼び掛けた。

確かにこの7人はecmonetの東京都でのECMO人数と一致している。
ではこのECMO 7人が全員若年だから若年の重症度化率が高まっているかというと、そもそもECMOを70歳以上の高齢者にしたところで、費用対効果として意味があるのかどうか、というのはそもそも新型肺炎が大流行する直前の2020年2月くらいで集中治療をする人たちにはぼんやり頭のなかにはあって、某意識高い系救急病院においては「65歳以上のVVECMO導入は見送る」なんて書いたマニュアルが裏では流れていた。

ecmonetの「COVID-19 ECMO症例(関西地方)の年齢変化」の図では

暫定的なデータですが、統計学的にも第四波では関西地区においてECMO実施症例の若年化が目立っております。全国的にはこの傾向はみられないので変異株の影響が強く疑われます。感染者も20才台30才台が多いとの報道もありますのでその影響も考えられます。いずれにしろ関西から全国にこの傾向が広がっていくことが懸念されます。2021/4/21

これはecmonetの関西におけるECMO年齢層の解析において、1-3波で103人平均67歳だったのが、4波で22人平均53歳になっていることを鑑みるに、VVECMOを行われる件数が増えたことと、感染者数の増加に伴って人工呼吸患者も増えて、それにより潜在的にVVECMOに移行しうる患者群が増えたので、高齢者へのVVECMOの差し控えが(裏では)行われているものと思われる。
つまりは単純に感染者が増えて確率的に40-50代の(若年と言われる)人でもECMOになる人が出てきたのでECMOが行われているだけであって、若年の重症化率がどうとかはあまり関係がないように思う。しかしながら高齢者へのECMOがなんとなくwithdrawal な風潮になっていそう(要検証)なのはそれはそれで残念ではあるが。

個人の感想を述べていてもエビデンスがないので、データから考えてみる。
週あたりの死亡者数は厚生労働省発表データおよび東洋経済のコロナページから取得できる。
GitHub - kaz-ogiwara/covid19: 新型コロナウイルス感染症(COVID-19)の国内における状況を厚生労働省の報道発表資料からビジュアルにまとめた。
国内の発生状況など|厚生労働省
新型コロナウイルス感染症の国内発生動向:2021年5月5日18時時点
感染者数の推移はNHKからデータが取得できる。
新型コロナウイルス 都道府県別の感染者数・感染者マップ|NHK特設サイト

厚生労働省のデータは1周間ごとくらいに重症者、死亡者数を年齢階級ごとにまとめているので、これを前処理して使用する。
重症度の定義が(真面目に下調べしてないから)曖昧であり、かつ一時欠損しているので、厳密にはECMOに依存している患者数ではないが、死亡数のデータはさすがに間違えないし欠損なく存在しているので、死亡数を扱うことである時点で死亡率が増減しているかどうかを検討する。
モデルとして死亡率\theta_{a, t} を推定する。年齢区分a=\{1,2,\dots,9\}、集計週t=\{0, 1, \dots,T\}において
・年齢区分が大きくなるほど死亡率が増加する:\theta_{a, t} \leq \theta_{a+1, t}
・週ごとの変化は2階差分モデル:\theta_{a, t+2} - \theta_{a, t+1} \sim normal(\theta_{a, t+1} - \theta_{a, t}, \sigma)
として、\theta_{a, t} から死亡数D_{a,t}は感染者数I_{a,t}から二項分布で発生するとするが、データは検査陽性になった時点でのカウントで、そこから酸素療法を要したり挿管になったりするのはだいたい1~2週間くらいの時差があり、挿管からECMOになったりするのはまた時差があるが、挿管からいきなりECMOだったり10日くらい挿管で粘ったりするので、このあたりは西浦先生法を使ったりして確率分布を使ったモデリングが必要かもしれないが、面倒くさいので「感染した週から2週間後に死亡する」とする。つまり、t での感染者はt+2 での死亡数と関連する、とすると
D_{a,t+2} \sim binomial(I_{a, t}, \theta_{a, t})
となる。
(もっとまじめにモデリングするなら、死亡確率のとある確率分布f_{a,t}を仮定してt での死亡数を何週か先まで足し合わせて\displaystyle\sum_t^TD_{a,t}f_{a,t}やるのがいいと思う。
ここからスクリプトを流用してrstanで記述する。
mikuhatsune.hatenadiary.com

結果としては、70代および80代以上が高い死亡率を示し、25%近くに及ぶ瞬間もある。
(ひとつのグラフに2つの軸を載せるのはあんまり好ましくないが)感染者数のピークに遅れること3-4週くらいで死亡率が急増している。
それに比べて20~50歳代はというと、グラフが潰れていて見えない。つまりそんなに死亡率は高くない。
f:id:MikuHatsune:20210511152230p:plain

変異株が日本に入ってきた瞬間もしくは、「若年の重症化率が高い」と言われ始めたのがいつの瞬間なのかあやふやなので、死亡率が低い2020年11月11日の週を比較として、2021年以降の死亡率\theta_{a,t} の増減を検討した。さすがに2021年になればイギリス株がどうたらこうたら言っていたような気がするので適当である。
差分をとったのは若年の死亡率はもともとものすごく低い値で、推定の都合でばらついた値が出たときに比が100倍(0.0001が0.01)となるとこれは誇張になりすぎる、と思ったので差分にしている。

2021年2月17日の週で死亡率の増加が最大となっている。70代においては10%、80代以上においては20%近く死亡率が増加している。60代でも5%近く増加している。
50代以下でも増加はしていることはしているが、0以上なだけであって0...% の増加である。それでも、感染者数が膨大に増えれば数としては影響が大きく見えるだろう。1000人/週が一桁増えて10000人/週になったとすれば、0.1%の増加であっても10人死亡するのでインパクトは大きく見える。
2回目の緊急事態宣言で感染者数はすぐさま減少しているが、重症化や死亡に関しては、すでに感染して現在治療中の人の転帰によるので、これは4週くらいの時差がある。緊急事態宣言をして1か月で解除します、といっても、1か月後の時点で重症化率や死亡数(率)を指標にしていては、このときに数が最大化されているので医療現場の実感としては一番しんどいときだろう。
f:id:MikuHatsune:20210511152243p:plain

もうちょっと真面目にモデルを考えたり、Interrupted time series analysis を使ったりしたほうがいいのではないかとか考えたがrstanがそれなりに動いたのでよしとする。

医療裁判の判例

読んだ。

医事法判例百選 第2版 (別冊ジュリスト 219)

医事法判例百選 第2版 (別冊ジュリスト 219)

  • 発売日: 2014/03/31
  • メディア: ムック
これに1989年と1976年に出た別冊ジュリスト医療過誤判例百選と医事判例百選の5冊がある。

法律素人だが非常に面白かった。麻酔関連の訴訟でなにかあるかと思って読み始めたが、医療水準の変化や法的解釈、民事と刑事裁判の違いなどいろいろ勉強になった。
5冊もあるが半世紀近く経っているので、内容は様々に変遷しており収録されている判例もかなり変わっている。特に1989年と1976年の判例集はもはや現代ではありえないような医療が行われているので必ずしも読む必要はないが、医療の歴史という意味では読む価値はあるかもしれない。特に国民皆保険制度が始まっての判例は正直なところ法律素人なので解釈とか判例の意味がわからなかったが、興味があれば読んでみるといいかもしれない。
次の判例百選が出たら普通に買ってしまうと思う。

Interrupted time series analysis でARIMAを使って真面目にやる

読んだ。
Interrupted time series analysis using autoregressive integrated moving average (ARIMA) models: a guide for evaluating large-scale health interventions

Interrupted time series analysis は介入の前後でlevel(増減、つまり切片)とtrend もしくはslope(単調な時系列増減、つまり傾き)を断片的に回帰するモデル(式は雑)
Y_t=\beta_0+\beta_t t +\beta_{intervention}X_{intervention}+\beta_{level}t+\beta_{trend}t
を推定することが多いが、時系列データなので普通は自己相関とか季節性とかそういうことを考慮しないといけないはずで、ARIMAを使ったという話で、実装上もRのforecast パッケージだけでできる、という話。

サンプルデータは2014年に処方規制がされたクエチアピン25mg錠の処方数をもとに解析していて、本文中にはstep(level のこと)とramp(trend のこと)のほかにpulse(介入が起きた直後にスパイク状に突然増える/減る効果がある場合があると思うが、これはすぐ次の時刻でもとに戻る)というものも想定している。また、介入後の増減の仕方はtransfer function という\omega_0 を用いてもっと複雑にできるが、これはhr のパラメータを設定する必要があり、これはこれでまた複雑な話なので本文中ではさらっと流され、しかもどちらもh=0r=0 とされている。
ARIMAのモデル(p, d, q)の決め方は、forecast パッケージのauto.arima関数を使ってAICもしくはBICで最適化されたものを勝手に採用してくれる。

解析結果としては、2014年の介入で処方数は3285、trendも1397/月減少したことがわかった。confint関数を使えば95%信頼区間も計算できるし、lmtest パッケージのcoeftest関数もしくはベタに計算すれば、点推定のp値も計算できる。

library(lmtest)
coeftest(model1)
z test of coefficients:

        Estimate  Std. Error  z value  Pr(>|z|)    
ar1     -0.87301     0.12396  -7.0427 1.885e-12 ***
ar2     -0.67314     0.12587  -5.3480 8.894e-08 ***
sma1    -0.60694     0.38723  -1.5674     0.117    
step -3284.77920   602.33653  -5.4534 4.942e-08 ***
ramp -1396.65226   106.63300 -13.0977 < 2.2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
(1-pnorm(abs(model1$coef)/sqrt(diag(model1$var.coef))))*2
         ar1          ar2         sma1         step         ramp 
1.884937e-12 8.894107e-08 1.170208e-01 4.941704e-08 0.000000e+00 

ARIMAモデルを使っているので、もし介入がなかった場合の仮想現実counterfactrual も推定できる。この場合は赤線になる。

医療事故の舞台裏

読んだ。

COI:中古で安かったので買った。

別に判例集を読んでいるが、一般的な判例集とは違い、臨床経過があって、場合によってはいくらか検査値データなどがあって、転機があって、それに対する患者・関係者の心情があって、訴えの内容、判決内容(無罪、賠償金など)が書いてあるので、法律素人向けにはわかりやすかった。
内容はたいはんがIC関連でこじれたようなものだが、薬剤投与や手技に関連した合併症のことなど多く、参考になった。特に、判例集などで超有名なもの以外の判例が挙げられているので、ほとんどが新鮮な気持ちで読める。

化合物の構造式を描きたい

という相談を受けた。
化合物なんてまったく扱ってないし構造式も大学入試以来描いたことがないが、こういうのはRで出来るだろうと思って調べたらあった。
ChemmineR: a compound mining framework for R | Bioinformatics | Oxford Academic
ChemmineR: Cheminformatics Toolkit for R
ChemmineRを使ってみよう【1】 - アメリエフの技術ブログ

SDFでファイルをもらったのでSDFでやってみる。
化合物の情報はpubchemというデータベースからいろいろ取ってこれるので、ここからフェンタニルを取ってきてプロットしてみる。
f:id:MikuHatsune:20210329122939p:plain

sdfset <- read.SDFset("https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/CID/3345/record/SDF/?record_type=3d&response_type=display")
plot(sdfset, print_cid=sdfset@SDF[[1]]@header["Molecule_Name"])

メタアナリシスで、サンプルサイズが最も大きいわけではないのに、weightが最大になるのはおかしくないですか?

という質問を受けた。
結論から言うとおかしくない。メタアナリシスのweight w_iは各研究内の分散v_i と研究間の分散\tauにより
w_i=\frac{1}{v_i}(fixed model の場合)もしくはw_i=\frac{1}{v_i+\tau}(random effect model の場合)
で決まるから、分散が小さい、すなわち推定精度の高い研究はweight が大きくなる。サンプルサイズが大きいと推定精度が増すので、サンプルサイズが大きいことがweight を大きくする要因だが、分散がそもそも小さいことも重要である。

とあるメタアナリシスで、サンプルサイズが114+115のKienast という研究(29%)より、サンプルサイズが33+23のBaudoという研究(57%)のほうがweight が大きい。
f:id:MikuHatsune:20210318132522p:plain

Baudo の研究だけ取り出してみると、分割表は
\begin{matrix}23&10\\20&3\end{matrix}
であるが、周辺度数を固定したときに、この分割表の取りうる運命は、\Delta を変動させて
\begin{matrix}23+\Delta&10-\Delta\\20-\Delta&3+\Delta\end{matrix}
となる。
\Deltaを変動させて、そのときの分割表のMH検定を行い、信頼区間の幅を見てみると、\Delta=-2のときに信頼区間の幅が最小になる。すなわち、推定の精度がよいことになる。
f:id:MikuHatsune:20210318132535p:plain

すべての\Delta の取りうる範囲が、各々メタアナリシスに組み入れられた研究だと仮定して、それら全体をメタアナリシスしてみる。確かに、信頼区間の幅の狭い\Delta=2などの研究が、weight が大きい、ということになっている。
f:id:MikuHatsune:20210318132730p:plain

library(meta)
library(epiR)
library(latex2exp)
# 元のデータ
dat <- cbind.data.frame(
   author=c("Baudo", "Fourrier", "Gando", "Inthom", "Kienast"),
   year=c(1998, 1993, 2013, 1997, 2006),
   ev.exp=c(23, 4, 3, 3, 29),
   n.exp=c(33, 17, 30, 7, 114),
   ev.cont=c(20, 9, 4, 7, 46),
   n.cont=c(23, 18, 30 ,7, 115)
   )

cols <- replace(rep("black", nrow(dat)), dat$author=="Baudo", "blue")

m <- metabin(ev.exp, n.exp, ev.cont, n.cont,
             method="Inverse",
             data=dat, studlab = paste(author, year),
             label.e="Antithrombin", label.c="Control",
             label.right="Favours Control", col.label.right="red",
             label.left = "Favours Antithrombin", col.label.left="blue")

forest(m, col.study=cols, layout="RevMan5", comb.fixed=FALSE)

# Gando に着目する
i <- 1
b <- rbind(c(dat[i, "ev.exp"], dat[i, "n.exp"]-dat[i, "ev.exp"]),
           c(dat[i, "ev.cont"], dat[i, "n.cont"]-dat[i, "ev.cont"]))

# MH test
A <- matrix(c(1, -1, -1, 1), 2)
r <- epi.2by2(b, method="cohort.count")
r$massoc$RR.strata.taylor

# Delta をすべて動かして検定する
x <- (-sort(b)[1]):(sort(b)[2])
bs <- mapply(function(z) b + A*z, x, SIMPLIFY=F)
r <- mapply(function(z) try(unlist(epi.2by2(z)$massoc$RR.strata.taylor)), bs, SIMPLIFY=FALSE)
g <- t(mapply(function(z) switch(any((is(z)%in%"try-error"))+1, z, rep(NA, 3)), r))
g <- matrix(unlist(g), nrow(g))
colnames(g) <- c("RR", "CI lower", "CI upper")

cols <- c("black", "blue", "red")
txt <- mapply(function(z, w) substitute(x~y~Delta, list(x=z, y=w)), c(b), c("+", "-", "-", "+"))
matplot(x, g, type="o", pch=15, lty=1, xlab=expression(Delta), ylab="RR", col=cols, lwd=3, log="y")
legend("topleft", legend=colnames(g), col=cols, lty=1, lwd=3)
legend("bottomright", legend=sapply(txt, as.expression), ncol=2)

# すべてのDelta に対応する分割表が各研究だとみなしてメタアナリシスする
bs1 <- sapply(bs, c)
bs0 <- cbind(bs1[1,], bs1[1,]+bs1[3,], bs1[2,], bs1[2,]+bs1[4,])
colnames(bs0) <- tail(colnames(dat), 4)
bs0 <- cbind.data.frame(year=x, bs0)
m <- metabin(ev.exp, n.exp, ev.cont, n.cont,
             method="Inverse",
             data=bs0, studlab=year,
             label.e="Antithrombin", label.c="Control",
             label.right="Favours Control", col.label.right="red",
             label.left = "Favours Antithrombin", col.label.left="blue",
             )

forest(m, layout="RevMan5", comb.fixed=FALSE)