R で推移行列モデル

今回は R で推移行列モデルを計算するコードを紹介します。
推移行列モデルの詳細は書籍：海洋ベントスの生態学を参照して下さい。



早速ですが、推移行列モデルとは以下のような数式で表されます。

$$\begin{align*} A = \begin{pmatrix} f_1 & f_2 & f_3 & ... & f_{w-1} & f_w \\ p_1 & 0 & 0 & ... & 0 & 0 \\ 0 & p_2 & 0 & ... & 0 & 0 \\ 0 & 0 & p_3 & ... & 0 & 0 \\ . & . & . & ... & . & . \\ . & . & . & ... & . & . \\ . & . & . & ... & . & . \\ 0 & 0 & 0 & ... & p_{w-1} & 0 \end{pmatrix} \end{align*}$$

ここで、$w$ はコホートの最高齢を表して、第 1 行目は各齢における繁殖量を表しているとのことです。$p_x$ は齢別生存率を表しており、例えば $p_1$ は 1 歳から 2 歳までの生残率を意味します。この値の算出には、予めコホート解析 (混合正規分布の分解) で得られたコホートの生残率が使えそうです。


次に、時刻 $t$ における個体群の年齢構成を以下のように記述します。

$$\begin{align*} n_t = \begin{pmatrix} n_1 \\ n_2 \\ n_3 \\ . \\ . \\ . \\ n_w \end{pmatrix} \end{align*}$$

ここでは、$n_1$ は 1 歳の個体数、$n_2$ は 2 歳の個体数...、というように表されています。
そして、これらの行列を使って、以下のように時刻 $t+1$ における年齢構成を計算します。

$n_{t+1}=A\times n_t$

そして、この $n_{t+1}$ を $n_t$ とし、繰り返し計算することで未来の年齢構成を予測することができるわけです。どうして掛け算するだけで次の世代が計算できるのかということに関しては、行列の計算方法を考えると分かると思います。





プログラムについて --------------------------------------------
基本的には推移行列と齢構成のベクトルを繰返し計算するコードです。

推移行列から個体群増加率 λ が算出でき、これが λ > 1 なら個体群は増加するそうです。この λ は推移行列の最大固有値だということなので、以下のコードでは eigen() 関数で算出しています。



自分のデータで解析するときは、


1. 計算したい世代数を入力する
2. 以下のコードの変数 LM に推移行列を入力する
3. 個体群の年齢構成の初期値 N を入力する (デフォルトでは各齢 10 個体)
4. 実行


という手順で行って下さい。

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

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

# 変数の設定 ---------------------------------------------------
# 計算したい世代数 ----------
f <- 100

# データ入力 ----------
LM <- c(0, 0.699, 3.754, 6.747, 4.414,
              0.505, 0, 0, 0, 0,
              0, 0.723, 0, 0, 0,
              0, 0, 0.553, 0, 0,
              0, 0, 0, 0.229, 0)

# 個体群の年齢構成の初期値 -----
RN <- sqrt(length(LM))
N <- rep(10, RN)
DT <- N


# 計算開始 ---------------------------------------------------
LM <- matrix(LM, ncol=RN, byrow=T)

# 行列の計算 -----
for(i in 1:(f-1)){
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)

# 固有値 λ --------------------------
eigen(LM)$values[1]

# 加入量のグラフ -------------------
k <- 1:10    # 最大 1:nrow(WR)
barplot(WR[k,1], xlab="世代数", ylab="加入量", names.arg=k)

​