2016年4月19日火曜日

R で推移行列モデル ~ Shiny の活用~

前回、推移行列モデルを R で行うコードを紹介しました。
各齢の生残率の状態が保持されると仮定すれば、将来の個体群の齢構成がどうなるかということが計算できました。

ですが、現実の問題として、保全施策や個体群の維持管理などを考えたとき、「どの齢の生残率を高く保っておかなければいけないか?」といったことを考える必要がもあると思います。

これを推移行列モデルで計算するなら、$p_1, p_2, p_3, …, p_{w-1}$ を変更し、再計算を繰り返さなければなりません。例えば、1 歳の個体を保護対象とせず全滅させてしまったらどうなるかということを計算するには、$p_1 = 0$ で再計算する必要があります。さらに、2 歳の個体を半分は保護するとか、1 歳から 4 歳までを満遍なく守るとか、1 歳から 4 歳にかけて守る量を減らしてみるとかいう計算をするとなると、その度に $p$ を打ち直して計算するという、探索的な計算をしないといけないわけです。





どう考えても面倒です。
そこで、R のパッケージ Shiny を活用します。


Shiny に関する説明は、他で説明されている方がいますので、そちらにお任せします。Shiny を実行して得られるものは web アプリです。このアプリは、こちらの入力に対する再計算結果やグラフ表示を返してくれます。こちらからの値の入力方法はいくつも用意されており、直接打ち込むタイプのテキスト入力や数値入力がある他、マウスで操作するスライダーバーやチェックボックス・ラジオボタン・セレクトボックスなどがあり、日付を入力する際はカレンダーを利用することもできます。数値や条件を入力すると即座に結果を返してくれるので、面倒で探索的な計算には非常に有効です。



使い方 --------------------------
まず、RStudio の File → New Project → New Directory → Shiny Web Application を選択します。そして、Directory name に任意の名前を入れて、この名前のフォルダを作る場所を Create project as subdirectory of の Browse ボタンで選びます。すると、RStudio に ui.R と server.R というタブができます。(この時点で Run App をクリックすると、ヒストグラムのバーの数を変更して表示するアプリが出てくるので、Shiny がどんなものなのか触ってみるといいと思います。)

次に、ui.R と server.R を書き換えます。
以下が推移行列モデルの ui.R と server.R ですので、それぞれをコピペして書き換えて下さい。
この ui.R はスライダーバーなどの入力パネルを作るコードで、server.R は計算するコードです。
server.R はほとんど R コードと一緒です。

最後に、両方とも書き換えたら Run App をクリックすると推移行列モデルのアプリが表示されます。


「よし!自分で作ろう!」という方はチュートリアルをご一読下さい。


ui.R

# ui.R

shinyUI(fluidPage(
  titlePanel(h1("推移行列モデル")),

  sidebarLayout(
    sidebarPanel(sliderInput("p1", "p1", min = 0, max = 1, value = 0.505, step = 0.001),
                 sliderInput("p2", "p2", min = 0, max = 1, value = 0.723, step = 0.001),
                 sliderInput("p3", "p3", min = 0, max = 1, value = 0.553, step = 0.001),
                 sliderInput("p4", "p4", min = 0, max = 1, value = 0.229, step = 0.001),
                 sliderInput("range", "グラフの範囲", min = 0, max = 100, value = 10)
                 ),
    mainPanel("加入個体数",
              plotOutput("recruit"))
  )
))


server.R

#作成日時 : 15:00 2015/10/20
#作成者名:T. Nakano
rm(list=ls())
gc();gc()

# ///////////////////////////////////////////////////////////////////

library(shiny)

shinyServer(function(input, output) {
  output$recruit <- renderPlot({

    # データ入力 ----------------------------------------------------
    LM <- c(0, 0.699, 3.754, 6.747, 4.414,
            input$p1, 0, 0, 0, 0,
            0, input$p2, 0, 0, 0,
            0, 0, input$p3, 0, 0,
            0, 0, 0, input$p4, 0)

    LM <- matrix(LM, ncol=sqrt(length(LM)), byrow=T)

    # 計算開始 ------------------------------------------------------
    DT <- N <- rep(10, ncol(LM))

    for(i in 1:input$range){
      act1N <- rep(N, nrow(LM))
      act1N <- matrix(act1N, ncol=length(N), byrow=T)

      N <- apply(LM*act1N, 1, sum)
      DT <- c(DT, N)
    }

    WR <- matrix(DT, ncol=length(N), byrow=T)

    # 加入量のグラフ --------------------------------------------------
    par(family="Meiryo")
    maintxt <- paste("固有値 =", round(abs(eigen(LM)$values[1]), 3))
    k<- 1:input$range
    barplot(WR[k,1], ylab="個体数", xlab="時刻 (年)", names.arg=k, main=maintxt)

  })

})





出力されるアプリの画面





グラフの範囲のスライダーバーを動かすと、棒グラフの x 軸が変更されます。
つまり、計算される世代数が変わるということです。

p1~p4 は推移行列モデルの $p_1 ~ p_4$ に対応しています。
このスライダーバーを動かすと、各齢の生残率を変更し、再計算されます。
つまり、p1~p4 を動かすだけで、探索的な計算ができるというわけです。


Shiny を使うことで解析の全てが完了するわけではありませんが、値の当たりをつけたり、解析の方向性を検討するということには有効なんじゃないかなぁ、と感じました。