Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2018-10-10

MikuHatsune2018-10-10

ようやくShinyを使ってインタラクティブにパラメータをいじって遊べるシミュレーターをweb上にあげることができた

Shinyというhtml 上でいろいろパラメータをいじったものを即座に計算して描出するインタラクティブでなんかすごいアプリがあったのだが、rmarkdown と同様に手をつけることなく放置していた。

せっかちな人はこちら → https://yfujii08.shinyapps.io/pkpdsimulator/

一番出来のいいエクセルは → https://home.hiroshima-u.ac.jp/r-nacamura/

 

話は変わるのだが、手術の最中に鎮痛(痛み止め)としてフェンタニルという麻薬を使う。この薬は適当に投与していてもなんとかなるが、本当にいい感じに使いたかったら、痛みがとれる濃度に収めつつ、かつ、副作用として呼吸抑制(息をしてくれない)を起こさない濃度にしなければならない。

そして、手術の最中は適当に投与していてもいい(手術中は人工呼吸器という機械で息をしてくれているのは保証されているから)が、手術が終わったあとも、上記のいい感じの濃度に保って鎮痛を継続して図りたいのである。

ここで、濃度のシミュレーションアプリとして、FentaSim(無料)、iTIVA anesthesia(無料といいながらアプリ内課金でぜんぜん無料で使えない)、AnestAssist PK/PD(潔く有料)というものがいくつか出ているのだが、これらのアプリがやってくれているのは、結局、

・濃度シミュレーション

UI として可視化

というわけである。濃度シミュレーションAnesthesiology. 1990 Dec;73(6):1091-102.あたりから必要な微分方程式パラメータを持ってきている。要は3コンパートメントモデルと効果部位を想定して、コンパートメントi からj への移動をk_{ij} と表記する。ここで効果部位はe, メイン(というか静脈注射されて薬が一気に流入してくるコンパートメントのことで、たいてい1番目)コンパートメントから他のコンパートメント以外への排泄を0 として、使われるパラメータk_{10}, k_{12}, k_{13}, k_{21}, k_{31}, k_{e0} となる。

というわけでフェンタニルパラメータを使って開始時と1時間のときに適当に投与した結果がこちら。

f:id:MikuHatsune:20181010204339p:image

# patient setting
weight <- 50
height <- 160
gender <- "M"
LBM <- ifelse(gender == "M", 1.1*weight-128*(weight/height)^2,1.07*weight-148*(weight/height)^2)
timestep <- 120

# parameter setting
k_1N <- c(0.0827, 0.471, 0.225)
k_N1 <- c(0.102, 0.006)
k_e0 <- 0.114

k <- c(k_1N, k_N1, k_e0)/(60/timestep)
names(k) <- c(sprintf("k1%d", c(0, 2, 3)), sprintf("k%d1", 2:3), "ke0")
k1 <- k["k10"] + k["k12"] + k["k13"]

V1 <- 0.105*weight
CLs <- k[c("k10", "k12", "k13")]*(60/timestep)*V1
V2 <- CLs[2]/(k["k21"]*(60/timestep))
V3 <- CLs[3]/(k["k31"]*(60/timestep))
# 適当な投与
X[1, "Bolus"] <- 123000*(60/timestep)
X[60, "Bolus"] <- 100000*(60/timestep)

X <- rbind(0, X)

for(i in 1:(nrow(X)-1)){
  X$dA1_1[i+1] <- k["k21"]*X$A2[i]+k["k31"]*X$A3[i]-k1*X$A1[i]+X$Cont[i+1]/(60/timestep)
  X$dA2_1[i+1] <- k["k12"]*X$A1[i] - k["k21"]*X$A2[i]
  X$dA3_1[i+1] <- k["k13"]*X$A1[i] - k["k31"]*X$A3[i]
  X$dAe_1[i+1] <- k["ke0"]*(X$A1[i] - X$Ae[i])
  X$dA1_2[i+1] <- k["k21"]*(X$A2[i]+X$dA2_1[i+1]/2)+k["k31"]*(X$A3[i]+X$dA3_1[i+1]/2)-k1*(X$A1[i]+X$dA1_1[i+1]/2)+X$Cont[i+1]/(60/timestep)
  X$dA2_2[i+1] <- k["k12"]*(X$A1[i]+X$dA1_1[i+1]/2) - k["k21"]*(X$A2[i]+X$dA2_1[i+1]/2)
  X$dA3_2[i+1] <- k["k13"]*(X$A1[i]+X$dA1_1[i+1]/2) - k["k31"]*(X$A3[i]+X$dA3_1[i+1]/2)
  X$dAe_2[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_1[i+1]/2 - (X$Ae[i]+X$dAe_1[i+1]/2))
  X$dA1_3[i+1] <- k["k21"]*(X$A2[i]+X$dA2_2[i+1]/2)+k["k31"]*(X$A3[i]+X$dA3_2[i+1]/2)-k1*(X$A1[i]+X$dA1_2[i+1]/2)+X$Cont[i+1]/(60/timestep)
  X$dA2_3[i+1] <- k["k12"]*(X$A1[i]+X$dA1_2[i+1]/2) - k["k21"]*(X$A2[i]+X$dA2_2[i+1]/2)
  X$dA3_3[i+1] <- k["k13"]*(X$A1[i]+X$dA1_2[i+1]/2) - k["k31"]*(X$A3[i]+X$dA3_2[i+1]/2)
  X$dAe_3[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_2[i+1]/2 - (X$Ae[i]+X$dAe_2[i+1]/2))
  X$dA1_4[i+1] <- k["k21"]*(X$A2[i]+X$dA2_3[i+1])+k["k31"]*(X$A3[i]+X$dA3_3[i+1])-k1*(X$A1[i]+X$dA1_3[i+1])+X$Cont[i+1]/(60/timestep)
  X$dA2_4[i+1] <- k["k12"]*(X$A1[i]+X$dA1_3[i+1]) - k["k21"]*(X$A2[i]+X$dA2_3[i+1])
  X$dA3_4[i+1] <- k["k13"]*(X$A1[i]+X$dA1_3[i+1]) - k["k31"]*(X$A3[i]+X$dA3_3[i+1])
  X$dAe_4[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_3[i+1] - (X$Ae[i]+X$dAe_3[i+1]))
  X$A1[i+1]    <- X$A1[i]+(X$dA1_1[i+1]+2*X$dA1_2[i+1]+2*X$dA1_3[i+1]+X$dA1_4[i+1])/6+X$Bolus[i+1]/(60/timestep)
  X$A2[i+1]    <- X$A2[i]+(X$dA2_1[i+1]+2*X$dA2_2[i+1]+2*X$dA2_3[i+1]+X$dA2_4[i+1])/6
  X$A3[i+1]    <- X$A3[i]+(X$dA3_1[i+1]+2*X$dA3_2[i+1]+2*X$dA3_3[i+1]+X$dA3_4[i+1])/6
  X$Ae[i+1]    <- X$Ae[i]+(X$dAe_1[i+1]+2*X$dAe_2[i+1]+2*X$dAe_3[i+1]+X$dAe_4[i+1])/6
  X$Cp[i+1]    <- X$A1[i+1]/V1/1000
  X$Ce[i+1]    <- X$Ae[i+1]/V1/1000
  X$C1[i+1]    <- X$A1[i+1]/V1/1000
  X$C2[i+1]    <- X$A2[i+1]/V2/1000
  X$C3[i+1]    <- X$A3[i+1]/V3/1000
}
X <- X[-1,]


t0 <- Sys.time()
tx <- t0 + (seq(nrow(X))-1)*60
yl <- c(0, max(X$Ce))
par(mar=c(5, 5, 1, 1), las=1)
plot(tx, X$Ce, type="l", ylim=yl)


さて、自分のPCが手元にあるときは自由パラメータを設定して任意にシミュレーションできるが、PCがないときやRを使わない他の人に使ってもらいたいときに困る。このようなときにアプリ化しておくと便利である。これはShiny でできる。

Shiny のサンプルはRstudio の新規作成 > shiny にすれば、ヒストグラムを描いてbin をインタラクティブに操作できるものが作成される。出力はoutput$hoge に入っていて、サンプルではoutput$distPlot でヒストグラムが翔ようになっている。ここをいじる。

output にはinput$bins から必要な変数が取得できるので、ここをいじる。生値を入力するのか、スライドバーで設定するのか、選択肢を作るのかはhogehogeInput 関数で設定できる。

アプリを作るにはui.R とapp.R を作る必要があるようだが、app.R 内にui とserver をごちゃごちゃ書きこめばなんとかなるっぽい。アプリを公開するには手っ取り早いのはshinyserver を使う。rmd をRpubs に公開するのと同じように、Rstudio から run app したあとに出力されるhtml の右上に公開、というボタンがあるのでそれを押すと、「アカウントある?」と聞かれるのでとりあえずアカウントを作る。アカウントを作成したあとは、アプリ公開のためのキーを取得して紐づけする必要があるので、SECRET を直して他人に見られないようにRにコピペする。これはツイッターAPI などを使う人はやったことがあるはず。

そのあと、Rstudio から公開しようとするとなぜかできなかった。なのでRコンソールから

rsconnect::deployApp('path/to/your/app')

とするとなぜかうまくいった。このとき、app はapp.R ではなく、app.R があるディレクトリを参照する。

f:id:MikuHatsune:20181010204641p:image

### app.R
#
# This is a Shiny web application. You can run the application by clicking
# the 'Run App' button above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#
# http://shiny.rstudio.com/gallery/
# https://shiny.qhmqk.me/shinywidgets/

library(shiny)
library(shinyTime)
# rsconnect::deployApp('path/to/your/app')

# Define UI for application that draws a histogram
ui <- fluidPage(
   # Application title
   titlePanel("Fentanyl PK/PD"),
   # Sidebar with a slider input for number of bins 
   sidebarLayout(
      sidebarPanel(
         sliderInput("hour", "Hour", min = 1, max = 12, value = 4, step=1),
         #dateInput("date", label = h3("Date input"), value = Sys.Date()),
         timeInput("starttime", label = "Simulation start time", value = strptime("09:00:00", "%T"), seconds=FALSE),
         hr(),
         radioButtons("gender", "Gender:", c("Male" = "M", "Female" = "F")),
         sliderInput("age", label="Age", value = 40, min = 0, max = 120, step=1),
         numericInput("height", label = "Height [cm]", min=1, max=250, value = 150),
         numericInput("weight", label = "Weight [kg]", min=1, max=150, value = 50),
         hr(),
         timeInput("ftime1", label = "Fentanyl time 1", value = strptime("09:00:00", "%T"), seconds=FALSE),
         numericInput("iv1", label = "Fentanyl 1 [ng]", min=0, max=1000, value = 50),
         timeInput("ftime2", label = "Fentanyl time 2", value = strptime("10:00:00", "%T"), seconds=FALSE),
         numericInput("iv2", label = "Fentanyl 2 [ng]", min=0, max=1000, value = 0),
         timeInput("ftime3", label = "Fentanyl time 3", value = strptime("11:00:00", "%T"), seconds=FALSE),
         numericInput("iv3", label = "Fentanyl 3 [ng]", min=0, max=1000, value = 0),
         timeInput("ftime4", label = "Fentanyl time 4", value = strptime("12:00:00", "%T"), seconds=FALSE),
         numericInput("iv4", label = "Fentanyl 4 [ng]", min=0, max=1000, value = 0),
         timeInput("ftime5", label = "Fentanyl time 5", value = strptime("13:00:00", "%T"), seconds=FALSE),
         numericInput("iv5", label = "Fentanyl 5 [ng]", min=0, max=1000, value = 0),
         hr(),
         #numericInput("ylim", label = h3("ylim [ng/ml]"), min=0, max=50, value=3),
         sliderInput("ylim", label="ylim [ng/ml]", value = 0, min = 0, max = 5, step=0.5),
         sliderInput("Cont_conc", label="ivPCA [ng/hr]", value = 0, min = 0, max = 50, step=5),
         timeInput("ivPCAstart", label = "ivPCA start time", value = strptime("12:00:00", "%T"), seconds=FALSE),
         hr(),
         fluidRow(column(3, verbatimTextOutput("value")))
         #submitButton("Submit")
      ),
      # Show a plot of the generated distribution
      mainPanel(
         plotOutput("distPlot", click="plot_click"),
         verbatimTextOutput("info"),
         p("やりたいこと", style = "font-family: 'times'; font-size: 16pt"),
         p("フェンタニルの投与回数が5回までなのでsidebarPanel をどうやれば任意に追加できるか"),
         p("フェンタニル以外の薬剤の並列表記、切り替え"),
         p("プロット上にカーソルを持ってくると時刻と濃度のインタラクティブな表示"),
         strong("Strong() makes bold text."),
         em("And em() makes italicized (i.e, emphasized) text."),
         br(),
         code("code displays your text like computer code"),
         div("Use span and div to create segments of text that all have a similar style. For example, this division of text is all blue because I passed the argument 'style = color:blue' to div", style = "color:blue"),
         br(),
         p("span does the same thing, but it works with",
           span("groups of words", style = "color:blue"),
           "that appear inside a paragraph."),
         br(),
         p("https://home.hiroshima-u.ac.jp/r-nacamura/")
      )
   )
)

# Define server logic required to draw a histogram
server <- function(input, output) {
   test <- 1
   output$distPlot <- renderPlot({
      # generate bins based on input$bins from ui.R
      weight <- input$weight
      comp <- c(1:3, "e")
      cname <- c("Bolus", "Cont", sprintf("A%s", comp), mapply(function(z) sprintf("dA%s_%d", comp, z), 1:4), "Cp", sprintf("C%s", comp))
      X <- matrix(0, input$hour*60, length(cname), dimnames=list(NULL, cname))
      X <- as.data.frame(X)
      timestep <- 60
      k_1N <- c(0.0827, 0.471, 0.225)
      k_N1 <- c(0.102, 0.006)
      k_e0 <- 0.114
      
      k <- c(k_1N, k_N1, k_e0)/(60/timestep)
      names(k) <- c(sprintf("k1%d", c(0, 2, 3)), sprintf("k%d1", 2:3), "ke0")
      k1 <- k["k10"] + k["k12"] + k["k13"]
      
      V1 <- 0.105*weight
      CLs <- k[c("k10", "k12", "k13")]*(60/timestep)*V1
      V2 <- CLs[2]/(k["k21"]*(60/timestep))
      V3 <- CLs[3]/(k["k31"]*(60/timestep))
      t0 <- input$starttime - 61*0
      tx <- t0 + (seq(nrow(X))-1)*60
      X[tx > input$ivPCAstart, "Cont"] <- input$Cont_conc/60*1000
      X <- na.omit(X)
      #X[1, "Bolus"] <- input$iv1*1000*(60/timestep)
      #X[100, "Bolus"] <- 150000*(60/timestep)
      X <- rbind(0, X)
      ### together
      #ivs <- unlist(mapply(function(z) eval(parse(text=z)), sprintf("input$iv%d", 1:5)))
      #ft <- unlist(mapply(function(z) eval(parse(text=z)), sprintf("input$ftime%d", 1:5)))
      #idx <- mapply(function(z) tail(which(tx < z), 1), ft)
      #X[idx, "Bolus"] <- ivs*(60/timestep)
      ### each
      X[tail(which(tx < input$ftime1+61), 1), "Bolus"] <- input$iv1*1000*(60/timestep)
      X[tail(which(tx < input$ftime2), 1), "Bolus"] <- input$iv2*1000*(60/timestep)
      X[tail(which(tx < input$ftime3), 1), "Bolus"] <- input$iv3*1000*(60/timestep)
      X[tail(which(tx < input$ftime4), 1), "Bolus"] <- input$iv4*1000*(60/timestep)
      X[tail(which(tx < input$ftime5), 1), "Bolus"] <- input$iv5*1000*(60/timestep)
      for(i in 1:(nrow(X)-1)){
        X$dA1_1[i+1] <- k["k21"]*X$A2[i]+k["k31"]*X$A3[i]-k1*X$A1[i]+X$Cont[i+1]/(60/timestep)
        X$dA2_1[i+1] <- k["k12"]*X$A1[i] - k["k21"]*X$A2[i]
        X$dA3_1[i+1] <- k["k13"]*X$A1[i] - k["k31"]*X$A3[i]
        X$dAe_1[i+1] <- k["ke0"]*(X$A1[i] - X$Ae[i])
        X$dA1_2[i+1] <- k["k21"]*(X$A2[i]+X$dA2_1[i+1]/2)+k["k31"]*(X$A3[i]+X$dA3_1[i+1]/2)-k1*(X$A1[i]+X$dA1_1[i+1]/2)+X$Cont[i+1]/(60/timestep)
        X$dA2_2[i+1] <- k["k12"]*(X$A1[i]+X$dA1_1[i+1]/2) - k["k21"]*(X$A2[i]+X$dA2_1[i+1]/2)
        X$dA3_2[i+1] <- k["k13"]*(X$A1[i]+X$dA1_1[i+1]/2) - k["k31"]*(X$A3[i]+X$dA3_1[i+1]/2)
        X$dAe_2[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_1[i+1]/2 - (X$Ae[i]+X$dAe_1[i+1]/2))
        X$dA1_3[i+1] <- k["k21"]*(X$A2[i]+X$dA2_2[i+1]/2)+k["k31"]*(X$A3[i]+X$dA3_2[i+1]/2)-k1*(X$A1[i]+X$dA1_2[i+1]/2)+X$Cont[i+1]/(60/timestep)
        X$dA2_3[i+1] <- k["k12"]*(X$A1[i]+X$dA1_2[i+1]/2) - k["k21"]*(X$A2[i]+X$dA2_2[i+1]/2)
        X$dA3_3[i+1] <- k["k13"]*(X$A1[i]+X$dA1_2[i+1]/2) - k["k31"]*(X$A3[i]+X$dA3_2[i+1]/2)
        X$dAe_3[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_2[i+1]/2 - (X$Ae[i]+X$dAe_2[i+1]/2))
        X$dA1_4[i+1] <- k["k21"]*(X$A2[i]+X$dA2_3[i+1])+k["k31"]*(X$A3[i]+X$dA3_3[i+1])-k1*(X$A1[i]+X$dA1_3[i+1])+X$Cont[i+1]/(60/timestep)
        X$dA2_4[i+1] <- k["k12"]*(X$A1[i]+X$dA1_3[i+1]) - k["k21"]*(X$A2[i]+X$dA2_3[i+1])
        X$dA3_4[i+1] <- k["k13"]*(X$A1[i]+X$dA1_3[i+1]) - k["k31"]*(X$A3[i]+X$dA3_3[i+1])
        X$dAe_4[i+1] <- k["ke0"]*(X$A1[i]+X$dA1_3[i+1] - (X$Ae[i]+X$dAe_3[i+1]))
        X$A1[i+1]    <- X$A1[i]+(X$dA1_1[i+1]+2*X$dA1_2[i+1]+2*X$dA1_3[i+1]+X$dA1_4[i+1])/6+X$Bolus[i+1]/(60/timestep)
        X$A2[i+1]    <- X$A2[i]+(X$dA2_1[i+1]+2*X$dA2_2[i+1]+2*X$dA2_3[i+1]+X$dA2_4[i+1])/6
        X$A3[i+1]    <- X$A3[i]+(X$dA3_1[i+1]+2*X$dA3_2[i+1]+2*X$dA3_3[i+1]+X$dA3_4[i+1])/6
        X$Ae[i+1]    <- X$Ae[i]+(X$dAe_1[i+1]+2*X$dAe_2[i+1]+2*X$dAe_3[i+1]+X$dAe_4[i+1])/6
        X$Cp[i+1]    <- X$A1[i+1]/V1/1000
        X$Ce[i+1]    <- X$Ae[i+1]/V1/1000
        X$C1[i+1]    <- X$A1[i+1]/V1/1000
        X$C2[i+1]    <- X$A2[i+1]/V2/1000
        X$C3[i+1]    <- X$A3[i+1]/V3/1000
      }
      X <- X[-1,]
      
      yl <- c(0, ifelse(input$ylim==0, max(X$Ce), ifelse(input$ylim<max(X$Ce), input$ylim, max(X$Ce))))
      par(mar=c(5, 5, 2, 1), las=1, cex.lab=1.5)
      plot(tx, X$Ce, type="n", ylim=yl, xlab="Time", ylab="Effect concentratoin [ng/ml]")
      lines(tx, X$Ce, lwd=7)
      abline(h=1.0, lty=3)
      abline(h=2.0, lty=3, col="red", lwd=2)
      abline(v=seq.POSIXt(tx[1], tail(tx, 1)+100, by=15*60), lty=3)
      if(input$Cont_conc > 0){
        x_ivPCA <- as.numeric(input$ivPCAstart)
        abline(v=x_ivPCA, lty=3, col="yellow")
        text(x_ivPCA, par()$usr[4], "ivPCA", pos=3, xpd=TRUE, cex=1.6)
      }
   })
   output$info <- renderText({
     xy_str <- function(e) {
       if(is.null(e)) return("Click figure\n")
       click_time <- format(as.POSIXlt(as.numeric(e$x), origin="1970-01-01"), "%H:%M")
       sprintf("Ce is %.4f at %s", e$y, click_time)
       #paste0("x=", round(e$x, 1), " y=", round(e$y, 1), "\n")
     }
     paste0(
       xy_str(input$plot_click)
     )
   })
}

# Run the application 
shinyApp(ui = ui, server = server)

2018-01-17

MikuHatsune2018-01-17

(サッカー解説)2点差は危険なスコアですね ← ???

高校サッカーを見ていた。2017年度は前橋育英が初優勝で幕を閉じた。

どの試合だったか忘れてしまったが、2点差がついたときに解説が「2点差は危険」ということを言っていた。

 

調べてみると、やはりよく言われていることのようだが、実際にデータをとってみると、プレミアリーグでは2点差をひっくり返して勝つ確率は1.71%、Jリーグでは2点差からドローが5%、2点差から敗北5%だったらしい。

 

自分はユース年代のファンなので、せっかくJFA が公式に試合記録を出してくれるということもあって、冬の高校選手権の得点時間を抽出して、2点差が危険なのかどうかを解析したい。

 

JFA から公式記録PDF を取得するが、2009年(88回大会)から2017年(96回大会)まで存在していて、各大会47試合ある。ただし、2009年はPDF の都合でデータをパースできなかったので全部で379試合が対象である。

試合記録から、得点時間と得点チームを抽出する。というのも、得点時間の分布と、ある時刻でのリードが最終的に勝率にどのくらい影響するかをみたいからである。

ただし、冬の選手権は決勝までは80分+即PK戦、決勝は90分+15分ハーフの延長戦なので、いい感じに分類したかったが、90分のときのロスタイムでうまく処理できてない気もするが無視。

また、2点差が危険の定義だが、N点差がついていたにもかかわらず、追いつかれて点差が0 になったときとする。というのも、追いつかれたとしてもトーナメント方式なのでPK戦までして勝敗を決めるが、得点時間を取得した都合で、80分内での点差の推移を考えることにする。なので、80分時点での点差が勝ち、引き分け、負けになる。

90分の試合を補正して80分にしてもよかったが、していない。ロスタイムの得点で最大80+6分というのがある。

 

結論から言うと、2点差がつくとたいてい90% 以上の確率でリードを保ったまま試合を終えることができる。

ある試合時刻t のとき、2点差がついている試合数N_{t,2} があり、残りの80-t 時間、リードを保ったまま(点差が0 より大)試合を終えた試合数N_{t,2}^{’} を単純に数えて割ったものが各点であり、指数関数モデルで曲線フィッティングしたのが黒線である。

試合経過時間が短いときに得点しても、追いつくチャンスがまだまだあるので、リードを保ち続けられる確率は低めだが、時間が経てばそのチャンスが減っていくので、リードを保てる確率は増える。

1点差のリードのとき、80-85分のあたりで確率が低くなっており、このせいで推定曲線が全体的に下がったものが推定されているが、ここを除外して推定曲線を引くと、80分時点で1点リードしていればほぼ勝率が100% になる。

f:id:MikuHatsune:20180117130322p:image

 

前半、後半について得点時間の分布を見てみると、特にどの時間帯でピークがある、という感じではなく、一様分布のようである。なのでどの時間帯にトイレにいっても得点シーンを見逃す確率は同じ程度と思われる。

f:id:MikuHatsune:20180117130349p:image

 

一方で、得点が入ってから次の得点がはいるまでの時間の分布は、一様分布ではない。横軸が経過時間、縦軸が確率密度とすれば、ガンマ分布のようである。また、1試合あたりの得点数も、ポアソン分布っぽい。

というわけで、次の得点までの時間分布をガンマ分布で、1試合あたりの両チーム合計の得点数をポアソン分布でモデル化してみると、こんな感じになる。

ガンマ分布のパラメータは(1.33, 0.06) で、最頻5分、中央16分程度、平均20分で点がはいる。

ポアソン分布のパラメータは2.9であり、1試合あたり3点くらいの点がはいる。

f:id:MikuHatsune:20180117130440p:image

 

というわけで、「2点差は危険」というのは、1回のゴールで1点しか点が動かないサッカーのルールの特性上、2点差から追いつくもしくは逆転されるという試合が伝説のように語られる思い出バイアスなので、解説の人は勉強してください()

 

解析

### analysis
library(stringr)
dat <- read.table("soccer.txt", header=TRUE, stringsAsFactors=FALSE)

tmp <- split(dat[,-1], dat$year)
tmp <- mapply(function(z) split(z[,-1], z$no), tmp)

res <- vector("list", length(tmp))
k <- c(1, -1)
for(i in seq(tmp)){
  out <- vector("list", length(tmp[[i]]))
  for(j in seq(tmp[[i]])){
    hoge <- tmp[[i]][[j]]
    foo <- lapply(str_extract_all(hoge$time, "\\d+"), as.numeric)
    time <- mapply(function(z) append(z, rep(0, 2-length(z))), foo)
    team <- k[c(factor(hoge$team, levels=unique(hoge$team)))]
    bar <- rbind(time, team)
    rownames(bar) <- c("time", "extra", "team")
    out[[j]] <- bar
  }
  res[[i]] <- out
}


y0 <- lapply(res, do.call, what=cbind)
y1 <- do.call(cbind, y0)
half <- c("前半", "後半")
idx <- list(y1[1,] <= 40, y1[1,] >  40)
goaltime <- mapply(function(z) colSums(y1[1:2, z]), idx)
yl <- c(0, 4)
par(mfcol=c(2, 2), las=1)
for(i in seq(goaltime)){
  pl1 <- table(goaltime[[i]])/sum(unlist(idx), na.rm=TRUE)
  pl2 <- table(goaltime[[i]])/length(goaltime[[i]])
  barplot(pl1*100, xlab="経過時間", ylab="percentage", ylim=yl, col=i+1, main=sprintf("%sの得点時間分布 (%s時間)", half[i], "全試合"))
  barplot(pl2*100, xlab="経過時間", ylab="percentage", ylim=yl, col=i+1, main=sprintf("%sの得点時間分布 (%s時間)", half[i], half[i]))
}


score <- vector("list", sum(sapply(tmp, length)))
k <- 1
for(i in seq(tmp)){
  for(j in seq(tmp[[i]])){
    hoge <- res[[i]][[j]]
    score[[k]] <- cbind(c(0, colSums(hoge[1:2,,drop=FALSE])), cumsum(c(0, hoge[3,,drop=FALSE])))
    k <- k + 1
  }
}

mat <- NULL
for(i in seq(res)){
  for(j in seq(res[[i]])){
    k <- res[[i]][[j]]
    if(!is.na(k[1])){
      l <- append(rep(NA, max(k[1,])), rep(NA, tail(k[2,], 1)+1))
      elp <- c(1, colSums(k[1:2,,drop=FALSE])+1) # 得点時間のindex
      elp <- append(elp, ifelse(max(k[1,])<=80, max(80, elp), ifelse(max(k[1,])<=90, max(90, elp), 120)))
      s <- c(0, cumsum(k[3,]))
      for(m in 1:(length(elp)-1)){
        l[ elp[m]:elp[m+1] ] <- s[m]
      }
    } else {
      l[1:80] <- 0
      l <- rep(0, 80)
    }
    l <- append(l, rep(NA, 120-length(l)))
    mat <- rbind(mat, l)
  }
}
library(rstan)
# 同時に
code <- "
data{
  int Ngame;
  int M;
  int Ngoal[Ngame];
  vector<lower=-1>[M] G[Ngame];
}
parameters{
  real<lower=0, upper=10> lambda;
  real<lower=0, upper=5> p[2];
}
model{
  Ngoal ~ poisson(lambda);
  for(i in 1:Ngame){
    G[i, 1:Ngoal[i] ] ~ gamma(p[1], p[2]);
  }
}
"

m <- stan_model(model_code=code)
Ng <- sapply(goaltime, length)
game <- t(mapply(function(z) append(z, rep(-1, max(Ng)-length(z))), goaltime))
dat <- list(Ngame=length(Ng), Ngoal=Ng, M=max(Ng), G=game)
fit <- sampling(m, data=dat, iter=3000, warmup=1500, seed=1234)
ex <- extract(fit)


par(mfrow=c(1, 2), las=1, mar=c(5, 4.5, 2, 2), cex.lab=1.5)
goaltime <- mapply(function(z) head(rle(as.numeric(na.omit(mat[z,])))$length, -1), seq(nrow(mat)))
g <- unlist(goaltime)
x <- 1:100
tab <- table(factor(g, x))

dg <- dgamma(x, median(ex$p[,1]), median(ex$p[,2]))
plot(tab/length(g)*100, xlab="直前のゴールからの経過時間(分)", ylab="頻度 [%]")
lines(x, dg*100, col=2, lwd=3)

x <- 0:max(Ng)
Ngd <- table(factor(Ng, x))
dp <- dpois(x, median(ex$lambda))
plot(Ngd/length(Ng)*100, xlab="1試合の総得点数", ylab="頻度 [%]")
lines(x, dp*100, col=2, lwd=3)


# ある時刻での勝利確率を求める
mat80 <- mat[is.na(mat[,90]),]
s <- replicate(2, matrix(0, 2, 90), simplify=FALSE)
for(k in 1:2){
  for(i in 1:nrow(mat80)){
    a <- abs(mat80[i,])
    M <- sum(!is.na(a))
    for(j in 1:M){
      if(a[j] >= k){
        s[[k]][2, j] <- s[[k]][2, j] + 1
        if(all(a[j:M] > 0)){
          s[[k]][1, j] <- s[[k]][1, j] + 1
        }
      }
    }
  }
}

xl <- c(0, 90)
yl <- c(0.5, 1)
cols <- c("blue", "red")
plot(0, type="n", xlim=xl, ylim=yl, xlab="試合経過時間(分)", ylab="試合終了時までリードし続ける確率", cex.lab=1.2, cex.axis=1.2, las=1)
for(i in seq(s)){
  v <- s[[i]][1,]/s[[i]][2,]
  points(v, col=cols[i], pch=16)
  d <- drm(v ~ seq(v), fct=AR.3())
  lines(d$data[,1], d$predres[,1], lwd=3)
}
legend("bottomright", legend=sprintf("%d 点差", 1:2), col=cols, pch=16, cex=2)

 

データ取得と前処理

# 2017-2015 年までの96-93 大会
wd <- "/soccer/"
dir.create(wd)
y <- 2017:2014
m <- 1:47

for(i in y){
  urls <- sprintf("http://www.jfa.jp/match/alljapan_highschool_%d/match_report/m%d.pdf", i, m)
  dir.create(sprintf("%s%d", wd, i))
  for(j in urls){
    command <- sprintf("wget %s -P %s%d", j, wd, i)
    system(command)
  }
}

for(i in 2013){
  urls <- sprintf("https://www.jfa.or.jp/match/matches/2014/0113koukou_soccer/match_report/m%d.pdf", m)
  dir.create(sprintf("%s%d", wd, i))
  for(j in urls){
    command <- sprintf("wget %s -P %s%d", j, wd, i)
    system(command)
  }
}

y <- 2012:2008
d1 <- c("0114", "0110", "0110", "0110", "")
for(i in seq(y)){
  urls <- sprintf("http://www.jfa.or.jp/match/matches/%d/%skoukou_soccer/schedule_result/pdf/m%02d.pdf", y[i]+1, d1[i], m)
  dir.create(sprintf("%s%d", wd, y[i]))
  for(j in urls){
    command <- sprintf("wget %s -P %s%d", j, wd, y[i])
    system(command)
  }
}
system(sprintf("mv %s%d/* %s%d/; rm -r %s%d", wd, 2008, wd, 2009, wd, 2008))

###
library(pdftools)
library(stringr)
year <- list.files(wd)
year <- as.character(2017:2013)
# 2013-2017
output <- NULL
for(j in seq(year)){
  print(sprintf("%s", year[j]))
  files <- list.files(sprintf("%s%s", wd, year[j]), pattern="pdf", recursive=TRUE, full.names=TRUE)
  for(f in seq(files)){
  text <- pdf_text(files[f])
  text <- strsplit(text, "\n")[[1]]
  i <- 1
  while(i < length(text)){
    if(any(strsplit(text[i], " ")[[1]] == "得点時間")){
      i <- i + 1
      res <- NULL
      while(length(grep("PK戦の経過", text[i])) < 1){
        tmp <- strsplit(gsub(" +", " ", text[i]), " ")[[1]][2:3] if(length(grep("分", tmp)) > 0){
          res <- rbind(res, tmp)
        }
        i <- i + 1
      }
      tmp <- table(res[,2])
      if(length(tmp) > 1){
        if(tmp[1] == tmp[2]){
          i <- i + 1
          a1 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
          a1 <- a1[a1 != " "]
          i <- i + 1
          a2 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
          a2 <- a2[a2 != " "]
          a3 <- c(a1[1], a2[1])[which.max(c(length(grep("○", a1)), length(grep("○", a2))))]
          res <- cbind(res, (res[,2] == a3)+0)
        } else {
          i <- i + 100
          a3 <- names(tmp)[which.max(tmp)]
          res <- cbind(res, (res[,2] == a3)+0)
        }
      } else {
        res <- cbind(res, 1)
      }
    }
    i <- i + 1
  }
  res <- switch(1 + (length(res)==1), res, matrix(c(NA, NA, res), 1))
  res <- cbind(year[j], f, res)
  output <- rbind(output, res)
  }
}

year <- as.character(2012:2009)
for(j in seq(year)){
  files <- list.files(sprintf("%s%s", wd, year[j]), pattern="pdf", recursive=TRUE, full.names=TRUE)
  for(f in seq(files)){
  text <- pdf_text(files[f])
  text <- strsplit(text, "\n")[[1]]
  if(length(text) > 0){
    i <- 1
    while(i < length(text)){
      pktext <- str_extract(text[i], "先.*先")
      if(!is.na(pktext)){
        pk <- as.numeric(str_extract_all(pktext, "\\d+")[[1]])
      }
      if(any(strsplit(text[i], " ")[[1]] == "得点チーム")){
        i <- i + 1
        res <- NULL
        while(length(grep("備考欄|戦評者氏名", text[i])) < 1){
          tmp <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
          tmp <- tmp[grep("分", tmp)]
          if(length(grep("分", tmp)) > 0){
            tmp <- strsplit(gsub("分", "分 ", tmp), " ")[[1]]
            res <- rbind(res, tmp)
          }
          i <- i + 1
        }
        tmp <- table(res[,2])
        if(length(tmp) > 1){
          if(tmp[1] == tmp[2]){
            i <- i + 1
            a1 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
            a1 <- a1[a1 != " "]
            i <- i + 1
            a2 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
            a2 <- a2[a2 != " "]
            a3 <- c(a1[1], a2[1])[which.max(c(length(grep("○", a1)), length(grep("○", a2))))]
            res <- cbind(res, (res[,2] == a3)+0)
          } else {
            i <- i + 100
            a3 <- names(tmp)[which.max(tmp)]
            res <- cbind(res, (res[,2] == a3)+0)
          }
        } else {
          res <- cbind(res, 1)
        }
      }
      i <- i + 1
    }
    res <- switch(1 + (length(res)==1), res, matrix(c(NA, NA, res), 1))
    res <- cbind(year[j], f, res)
    output <- rbind(output, res)
    } else {
      print(files[f])
    }
  }
}

2017-08-11

杏仁豆腐「あんかけどうふ」の待機列が牛歩だったので待ち行列モデルで考えてみる

C92 1日目の「あんかけどうふ」というサークルは、デレマスのみくにゃんなどで有名な杏仁豆腐氏の作品が出るサークルということで、叶姉妹のサークル、けものフレンズのたつき監督と並んで3大徹夜サクチケ転売屋の巣窟、という下馬評があったりなかったり。

 

結論から言うと始発初手で新刊だけは買えた。

詳細はこちらにある。

 

当時の状況を整理すると

あんかけどうふのサークルは東6ホールの端

東京ビックサイトによると、東ホールは90m*90m

待機列整理は3列。

0957 東4 を折り返したところに並んだとのツイートがある。

1016 始発の自分が列に並び始める。このとき2列目の最後尾は、東6ホールの真ん中くらい。

1105 2限から1限になる。

1206 1時間くらい牛歩だったがようやく折り返し地点。このとき有明病院の真正面だということに気づく。

1234 なんとか購入できた。このとき、目の前の集団は12行あって、列がさばけるのを測ってみたら58秒だった。しかしこのときは1限で新刊のみ、そして500円玉払いがほとんどなのでたぶん秒速。自分も3秒くらいだった。

1254 杏仁豆腐氏により完売宣言

 

これを待ち行列モデルで考えてみる。

まず、単位長さの行列に何人参加者がいるかだが、これによると専有面積0.3m2 で混んでいる、くらいの印象のようなので、肩幅0.4m として1m あたり0.75人いることになる。

ここで、東ホール1つあたり90m で、東4ホールまで伸びたけど有明病院の目の前なので、2ホール分+ 20mと簡単にしておくと、1列あたり200m なので1列あたり266人、これが3列と折り返し(2列分) なので自分が並び始めた1016の時点で1600 人いた。

また、0957 時点で1列分とちょっとの折り返しなので、220m としておくと、1000 の開場時点で880人、徹夜とサクチケ組がいた事になる。

さて、自分が並び始めた1016 までの16分間で、720 人増えたことになるので、新規に参加者が並ぶのは¥lambda=45 となる。

さて、列の捌け具合はどうだったかというと、2限だった間は牛歩だったのでよくわからない。自分がようやく買える順番に来た時、目の前の人たちは12行分いたので、3列で36人である。これは58秒で捌けたので、窓口の対応速度は¥mu=¥frac{36}{58/60}=12.4 だった。

 

あとは完売のタイミングから用意していた冊数など考えようと思ったけど余力があれば続くを書く。

東ホール1つ分の1列で120くらいが並んでいる、というのがたぶん目安になる、かもしれない(適当

2017-06-14

MikuHatsune2017-06-14

single cell RNA-seq のdropout

読んだ。

MAGIC: A diffusion-based imputation method reveals gene-gene interactions in single-cell RNA-sequencing data

コードはPython で書かれている。

computational flowcytometry のDana Peer ラボ。RNA-seq のデータ行列が取るであろう高次元空間の形からデータ点(細胞)を復元しようという雰囲気。

single cell RNA-seq ではmRNA の読みの精度の悪さのせいで、全体の80-90% くらいは0 という結果が返ってくるdropout という現象に悩まされる。有効なデータが10% 程度しかなく、他が0 のときに元のデータがどのようなものだったかを推定するのが必要な作業だが、MAGIC(Markov Affinity-based Graph Imputation of Cells) という、マルコフ連鎖と近傍のデータ点からもともと読まれるはずだったmRNA のデータと細胞を再構成してデータを復元する方法があるらしい。

 

f:id:MikuHatsune:20170613192643j:image

図はNat Methods. 2014 Jul; 11(7): 740–742.のFig.1a

 

n 遺伝子 m 細胞の行列D がある。このなかのほとんどはdropout のせいで0 であり、データを復元したD_{imputed}がほしい。細胞で眺めれば遺伝子発現のベクトルになるので、このベクトルで表現型が分類でき、ベクトルが近いものはその近傍に細胞が固まるはずである(diffusion)。グラフ理論的にある表現型(ベクトルパターン)をもつ細胞の近傍k 個を考える。

D からは簡単に距離行列D_{dist}が計算できる(適当にユークリッド)。そして、D_{dist}を適切なadaptive ガウスカーネル関数で変換するのだが、

A=e^{-(¥frac{D_{dist}}{¥sigma})^2}

というように変換してaffinity matrix を作る。このとき、¥sigma はkernel width と呼ばれ、このアルゴリズムでは非常に重要である。¥sigma 以下のD_{dist}は速やかに0 affinity になるため、相当近い細胞同士でないとneighborhood が構成できない。

A は対称行列でないのでM=A+A^T として強制的に対称行列化して、マルコフ過程を満たすために行和を1 になるようにする。つまり確率推移行列になるようにする。

ここで、t というべき乗パラメータを使って、D_{imputed}=M^tD とすると、目的だったD_{imputed}が得られる。ここでD_{imputed}は、dropout して80-90% が0 だったRNA-seq データが、もともと得られるはずだったRNA-seq リードカウントデータになっている、ということらしい。

(M^t は行列積ではなく、ただ単に要素のべき乗である)

 

scImpute: Accurate And Robust Imputation For Single Cell RNA-Seq Data

MAGIC とは異なり、発現の多寡とdropout するかしないかを確率分布からモデル化して、それを説明するパラメータを推定する。

データの取得のされ方について以下のようなモデルを考える。

f_{X_{i}}(x)=¥lambda_i ¥textrm{Gamma}(x;a_i,b_i)+(1-¥lambda_i)¥textrm{Normal}(x;¥mu_i,¥sigma_i)

¥lambda_i:dropout rate

a_i,b_i:ガンマ分布のパラメータ。ガンマ分布はdropout をモデル化しており、遺伝子発現としては高くてばらつきの少ないものだが、これが0になるのはdropout によるもの、という感じ。

¥mu_i,¥sigma_i:正規分布のパラメータ。正規分布は通常の遺伝子発現をモデル化しており、低めから中程度で、ばらつきの大きな遺伝子発現をしており、これが0 になるということはまあ普通に考えて遺伝子発現がないから0 である、という感じ。

これをEM アルゴリズムでゴリゴリ解く。ベイズでやってもいいと思う。

実際にi 遺伝子がj 細胞でdropout する確率は

d_{ij}=¥frac{¥lambda_i ¥textrm{Gamma}(x;a_i,b_i)}{¥lambda_i ¥textrm{Gamma}(x;a_i,b_i)+(1-¥lambda_i)¥textrm{Normal}(x;¥mu_i,¥sigma_i)}

となる(各パラメータはEM アルゴリズムで推定されたあとのハット付きである)

 

dropout 確率と、細胞遺伝子発現の隠れモデルが推定されたので、実際にdropout しているデータをimputation する。ある適当なdropout 確率閾値t により、細胞j についてimputation されないといけない遺伝子発現セットA_j=¥{i:d_{ij}¥geq t¥} と、dropout 確率閾値を下回っているため確実にデータとして信頼できる遺伝子発現セットB_j=¥{i:d_{ij}< t¥} に分けて、B_j を用いてimputation する。これはL_1 正則化で行い、

¥hat{¥beta}^{(j)}=¥underset{¥beta^{(j)}}{¥textrm{argmin}}¥frac{1}{|B_j|}||w_{B_j}^T(X_{B_j,j}-X_{B_j,-j}¥beta^{(j)})||_2^2+¥theta||¥beta^{(j)}||_1

となる発現量を復元する。

dropout するかしないかというのは、例えば行動定量学で鳥の生態を観察するような、鳥は生きているが観察されないのか、鳥は死んでいるので観察されないのか、というのをモデル化するzero inflated model 的な感じがした。

 

Nat Methods. 2014 Jul; 11(7): 740–742.

もともとはこれがいっぱい引用されている。

ベイズ的にやるのだが、モデル式としては、dropout はポアソン分布、amplification は負の二項分布からサンプリングされるとしている。負の二項分布を使うと裾が伸びて分散が大きい値を表現できるから、これが高発現というかamplification を意図しているのだと思う。

あるsubpopulation で平均してx の発現量がある遺伝子の事後確率を

p_{S}=E¥[¥Pi_{c¥in B}p(x|r_c,¥Omega_c)¥]

とする。ただし、BS のブートストラップサンプル、r はRPM (RNA-seq で得られるリードカウントを適当に補正したようなもので、なんらかの定量値と思えばいい)、¥Omega はエラー項であり、c の細胞について、p(x|r_c,¥Omega_c) は事後確率であり、dropout を観測する確率p_d を用いて

p(x|r_c,¥Omega_c)=¥p_d(x)p_{¥textrm{Poisson}}(x)+(1-p_d(x))p_{¥textrm{NB}}(x|r_c)

と書く。

というようなことをやっているので、結局はdropout で0 になる、その0 になるのが通常発現パターンか、高発現パターンなのかを別のそれっぽい分布から出てきたようにモデル化するっぽい。

2017-05-22

MikuHatsune2017-05-22

糖質制限ダイエットを始めたらたった1日で体重が2kg減った

親戚が集まったときの会で、ふいにそんなことを言われた。その人としては減量をがんばったということを言いたかったのだろうが、自分自身、糖質制限ダイエットに対してエビデンスを持っていないことと、1日で2kg 減るっていうのは、偶然「2日間連続で体重を測ったら2kg 減っていた」というポジティブな観測を強烈に覚えているだけの観測バイアスなので「水分だけでも1〜2kg 程度は体重の日内変動はある。1ヶ月間毎日体重測ってプロットしてから出直せ」と返してしまったのだが、本当に糖質制限ダイエットをするだけで1日で2kg 減ったらどうしよう、と思ったので、データで考えてみる。

 

ある程度の期間、それなりに管理された測定データが欲しいのだが、実は某氏の体重データがある。某氏はデータ取得と解析の知識があるので、データのとり方は「毎朝、起床時に自宅の体重計で測定する」というものらしい。出張などで測れないときは欠測である。170cm 半ばの某氏は、さすがにラーメンを食べ過ぎたので減量しようと思ったようである。ただし、食生活を変えたり、運動習慣を強化した、というようなダイエット行動がいまいちわからなかったのでよくわからない。

 

ここで、体重を測っていなかったとしても、決して体重が存在しなかった、というわけではない。観測していなかっただけである。そういう場合に、状態空間モデルを使って、「真の体重」というのもは常にあるが、それを観測するときに誤差が生じるし、場合によっては欠測にもなる、というのをstan でやる。

 

真の体重w は、ふたつ前までの状態の差分で決まるとする。このとき、差分に対して適当なパラメータが決まり、その分散は固定されたb という実数とすると

w_t-w_{t-1} ¥sim ¥textrm{normal}(w_{t-1}-w_{t-2}, b)

つまり

w_t ¥sim ¥textrm{normal}(2w_{t-1}-w_{t-2}, b)

である。

実際に観測した体重W (大文字) は、測定時に着ている服とか、体内の水分バランスとか、体重計のキャリブレーション具合とかで、なんらかの観測誤差を含むと考えられるので、適当な分散s を用いて

W_t¥sim ¥textrm{normal}(w_t, s)

でサンプリングされる、とする。毎回の観測の条件を揃えると、s は小さくて済む。

欠測は、rstan ではNA を受け付けないので、欠測の日を先にindex で作成しておいて、rstan のmodel ブロック内でfor loop で回避する。

 

b, s をcauchy 分布からサンプリングしているが、なるべくばらついて大げさな値が出て欲しかったのでやってみた。適当である。

 

¥hat{R}<1.1 でうまく収束したようだ。計算時間は1000 iter で1 chain 3分程度。

f:id:MikuHatsune:20170522195344p:image

1日での体重減少の事後分布を見てみると、中央値で0.03 kg 程度減少する。95% 信用区間では、実は有意に体重が減少するという結果ではないようだ。しかし、200日程度のデータがあるので、分布の正の領域と負の領域が半々でない以上、負の方向(痩せる方向)へ傾いていくのは妥当だろう。

quantile(c(apply(ex$w, 1, diff)), c(a/2, 0.5, 1-a/2))
       2.5%         50%       97.5% 
-0.08474003 -0.03705018  0.01160505 

何も考えず、データ観測期間内での体重減少から日割りすると、-0.04 kg 程度だった。

diff(range(dat$weight, na.rm=TRUE))/nrow(dat)
[1] 0.04056604

 

さて、ここで気になるのが、「1日で2kg 痩せる」というのがどれくらいの出来事なのだろうか、ということである。体重の減少は、「真の体重減少」に加えて、そのときの偶然の誤差を加味した「体重計の指し示す値」で決まる。ここでいう偶然の誤差は、先に言った通り、体重計のキャリブレーションや、そのときの服装、直前に食事をしたり排泄をしたかどうか、が含まれる。

この「誤差」の作用は、上のモデルのs によって表現される。ここから、2kg 減ったというのが95% 信頼区間(ここでは頻度主義的に) でどうか、というのを調べるには

qnorm(1-a/2, sd=s)
[1] 0.8670893

となり、0.87kg 以上の観測誤差があるとき、2kg の1日減量やばい、ということになる。よって、糖質制限ダイエットの真の1日体重減少量は1.2 kg 程度となる。

(体重の減少、が興味の対象なので、片側では? という話もあるが、このときは0.72 kg である)

単純に、糖質制限ダイエットの真の減量効果が1kg/day としよう。1ヶ月で30kg 痩せる。1ヶ月で30kg 痩せたら自分としては悪性腫瘍を疑って医療機関の受診を勧める。

 

というのはまあ笑い話なのだが、糖質制限はするがたんぱく質、脂質はどんだけ摂取してもよいらしい。糖質がないと脳内でエネルギーにはならないけど? と聞くと、午前中はツライらしく、脂質のエネルギーって9kcal/g だけど? と聞くと、脂質は肌の恒常性の維持に必要だからむしろたくさん摂取しろという。謎()

 

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
dat <- read.csv("weight.csv")


code <- "
data{
  int Day;
  real<lower=0, upper=200> W[Day];
  int Na;
  int na[Na];
}
parameters{
  real<lower=0, upper=200> w[Day]; # weight
  real<lower=0> b; # sigma for hidden model
  real<lower=0> s; # sigma for observation
}
model{
  for(i in 3:Day){
    w[i] ~ normal(2*w[i-1] - w[i-2], b);
  }
  for(i in 1:Na){
    W[ na[i] ] ~ normal(w[ na[i] ], s);
  }
  b ~ cauchy(0, 0.5);
  s ~ cauchy(0, 0.5);
}
"

m <- stan_model(model_code=code)

W = replace(dat$weight, is.na(dat$weight), 1)
standata <- list(Day=nrow(dat), W=W, Na=sum(!is.na(dat$weight)), na=which(!is.na(dat$weight)))

fit <- sampling(m, standata)

ex <- extract(fit)
a <- 0.05
w <- apply(ex$w, 2, median)
b <- median(ex$b)
s <- median(ex$s)
ci <- apply(ex$w, 2, quantile, c(a/2, 1-a/2))

par(mar=c(5, 4.5, 2, 2), cex.lab=1.5, cex.axis=1.5)
xi <- seq(nrow(dat))
yl <- range(dat$weight, ci, na.rm=TRUE)
plot(w, type="l", ylim=yl, xlab="Day", ylab="Weight [kg]")
polygon(c(xi, rev(xi)), c(ci[1,], rev(ci[2,])), col=grey(0.8), border=NA)
lines(w)
points(dat$weight, pch=16, col=4, cex=0.8)
points(which(is.na(dat$weight)), w[which(is.na(dat$weight))], pch=16, col=2)
legend("topright", legend=c("実測値", "欠測時"), col=c(4, 2), pch=16, cex=3)