読んだ。
COI：著者に知り合いがいたが買ってから気づいた。安いので初学者は買ってもいいかもしれない。
読んだ。
COI：安かったので買った。知り合いはまったくいない。新型肺炎で2回PCRしたけど陰性で、でも感染症の専門家がやっぱり疑わしいからダメ押しでもう一回PCRやって、というのでやったら陽性だった
と言われた。
検査には感度と特異度という検査特性が存在するので、検査が陽性だからといって本当にその疾患が陽性であるかは限らないし、検査が陰性だからといって本当に陰性であるとは限らない。
ではどうやって陽性もしくは陰性の判定を下すのか、というと、「その判定を下すに足るとみなされる閾値を超える/下回る場合に判定を下す」か、「その閾値になるまで検査を追加する」のが戦略である。
いま、新型肺炎のPCR検査の感度を、特異度を
、新型肺炎であろう事前確率を
とする。一回の検査で陰性であるとき、事後確率
は、陰性尤度比とオッズを用いて
とかける。
検査が連続して 回陰性であるとき、
回後の事後確率はオッズのままの表記を利用すると（
とした）
回までの反復検査において、事前確率と事後確率の関係は図の通りになる。
検査特性として特異度が高い検査は、陽性である場合に、本当に陽性である確率が高くなる。これはSpIn と呼ばれる。この原理によると、陽性者を補足することには非常に有用である。
一方で、感度が高い検査は、陰性である場合に、本当に陰性である確率が高くなる。これはSnOut と呼ばれる。この原理によると、陽性ではないことを主張するには有用である。ということで、感度70% のPCR検査ではこちらの主張はできないため、「（疾患の）陰性証明書」なるものが意味をなさない、というのが通説である。
を大きくすればいいだろう、と思うかもしれないが、たしかに
回すべてで陰性であれば、本当に陰性と言っていいかもしれないが、そのうち1回でも陽性になってしまった場合に真陽性なのか偽陽性なのか判断しにくくなるし、現実的な運用や費用の問題がある。2類感染症という社会的要因のため、1度でも陽性なら判定としては陽性、という運用になる。
例えば結核の3連痰においては、検査の感度と特異度は以下のとおりである。
3連痰の活動性肺結核に対する感度は約70%、特異度は90%以上
亀田感染症ガイドライン：結核を疑う時とその対応（version 2） - 亀田総合病院 感染症科
結核は空気感染するため公衆衛生上、診断をつけて適切な対応をすることが重要であり、3連痰のうち1度でも陽性なら、陽性と判定することになっている。
postp <- function(x, Sn=0.7, Sp=0.95, n=1){ z <- ((1-Sn)/Sp)^n z*x/(1-x+z*x) } x <- seq(0, 1, length=1000) Y <- mapply(function(z) postp(x, n=z), 1:5) cols <- jet.colors(ncol(Y)) par(mar=c(5, 5, 2, 2), cex.axis=1.5, cex.lab=1.5) matplot(x, Y, type="l", lty=1, col=cols, lwd=3, las=1, xlab="Prior probability", ylab="Posterior probability") legend("topleft", legend=sprintf("%d回陰性", seq(cols)), col=cols, pch=15, cex=2) abline(0, 1, lty=3)
ちょっとしたアプリケーションを作って、自分以外の、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年５月５日18時時点
感染者数の推移はNHKからデータが取得できる。
新型コロナウイルス 都道府県別の感染者数・感染者マップ｜NHK特設サイト
厚生労働省のデータは1周間ごとくらいに重症者、死亡者数を年齢階級ごとにまとめているので、これを前処理して使用する。
重症度の定義が（真面目に下調べしてないから）曖昧であり、かつ一時欠損しているので、厳密にはECMOに依存している患者数ではないが、死亡数のデータはさすがに間違えないし欠損なく存在しているので、死亡数を扱うことである時点で死亡率が増減しているかどうかを検討する。
モデルとして死亡率 を推定する。年齢区分
、集計週
において
・年齢区分が大きくなるほど死亡率が増加する：
・週ごとの変化は2階差分モデル：
として、 から死亡数
は感染者数
から二項分布で発生するとするが、データは検査陽性になった時点でのカウントで、そこから酸素療法を要したり挿管になったりするのはだいたい1～2週間くらいの時差があり、挿管からECMOになったりするのはまた時差があるが、挿管からいきなりECMOだったり10日くらい挿管で粘ったりするので、このあたりは西浦先生法を使ったりして確率分布を使ったモデリングが必要かもしれないが、面倒くさいので「感染した週から2週間後に死亡する」とする。つまり、
での感染者は
での死亡数と関連する、とすると
となる。
（もっとまじめにモデリングするなら、死亡確率のとある確率分布を仮定して
での死亡数を何週か先まで足し合わせて
やるのがいいと思う。
ここからスクリプトを流用してrstanで記述する。
mikuhatsune.hatenadiary.com
結果としては、70代および80代以上が高い死亡率を示し、25%近くに及ぶ瞬間もある。
（ひとつのグラフに2つの軸を載せるのはあんまり好ましくないが）感染者数のピークに遅れること3-4週くらいで死亡率が急増している。
それに比べて20～50歳代はというと、グラフが潰れていて見えない。つまりそんなに死亡率は高くない。
変異株が日本に入ってきた瞬間もしくは、「若年の重症化率が高い」と言われ始めたのがいつの瞬間なのかあやふやなので、死亡率が低い2020年11月11日の週を比較として、2021年以降の死亡率 の増減を検討した。さすがに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か月後の時点で重症化率や死亡数（率）を指標にしていては、このときに数が最大化されているので医療現場の実感としては一番しんどいときだろう。
もうちょっと真面目にモデルを考えたり、Interrupted time series analysis を使ったりしたほうがいいのではないかとか考えたがrstanがそれなりに動いたのでよしとする。
読んだ。
法律素人だが非常に面白かった。麻酔関連の訴訟でなにかあるかと思って読み始めたが、医療水準の変化や法的解釈、民事と刑事裁判の違いなどいろいろ勉強になった。
5冊もあるが半世紀近く経っているので、内容は様々に変遷しており収録されている判例もかなり変わっている。特に1989年と1976年の判例集はもはや現代ではありえないような医療が行われているので必ずしも読む必要はないが、医療の歴史という意味では読む価値はあるかもしれない。特に国民皆保険制度が始まっての判例は正直なところ法律素人なので解釈とか判例の意味がわからなかったが、興味があれば読んでみるといいかもしれない。
次の判例百選が出たら普通に買ってしまうと思う。
Interrupted time series analysis は介入の前後でlevel（増減、つまり切片）とtrend もしくはslope（単調な時系列増減、つまり傾き）を断片的に回帰するモデル（式は雑）
を推定することが多いが、時系列データなので普通は自己相関とか季節性とかそういうことを考慮しないといけないはずで、ARIMAを使ったという話で、実装上もRのforecast パッケージだけでできる、という話。
サンプルデータは2014年に処方規制がされたクエチアピン25mg錠の処方数をもとに解析していて、本文中にはstep（level のこと）とramp（trend のこと）のほかにpulse（介入が起きた直後にスパイク状に突然増える/減る効果がある場合があると思うが、これはすぐ次の時刻でもとに戻る）というものも想定している。また、介入後の増減の仕方はtransfer function という を用いてもっと複雑にできるが、これは
と
のパラメータを設定する必要があり、これはこれでまた複雑な話なので本文中ではさらっと流され、しかもどちらも
、
とされている。
ARIMAのモデルの決め方は、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 も推定できる。この場合は赤線になる。
