読者です 読者をやめる 読者になる 読者になる

RGRからゴンペルツ曲線を推定する

植物の重量wについて、瞬間的な成長速度は植物の重量に比例すると仮定する。

\frac{dw}{dt} = rw

このときrを相対成長速度(RGR: Relative Growth Rate)と呼ぶ。もしrが変化しないのであれば、植物は指数関数的に増大を続けてしまい、地球は程なくして滅ぶだろう。

w_t = w_0 e^{rt}

しかし、植物は無限に大きくなるわけではないので、一般的には相対成長速度は時間とともに減少する。ここに例えば適当な定数abを用意して、rは次のように減少すると仮定してみよう。

r= e^{a-bt}

先の式に代入すれば次のようになる。

\frac{dw}{dt} = we^{a-bt}

そして面倒なのでWolfram Alphaにまかせて微分方程式を問いてもらうとこうなる。

w_t = c_1 e^{-\frac{e^{a-b t}}{b}}

さらに定数項を少し整理して、c_2をこのように定義する。

c_2 = e^{-\frac{e^a}{b}}

そして得られた式をゴンペルツ関数(Gompertz function)と呼ぶ。

w_t = c_1c_2^{e^{-bt}}

瞬間的なRGRの推定

最初の仮定と若干違ってしまうのだが、短い時間\Delta tで植物の体重を測定したためにrが変化しなかったとしよう。すると、最初の式を\Delta tの間で積分することで次の式が得られる。

r = \frac{\Delta \ln w}{\Delta t}

実際の測定で推定できそうな雰囲気がしてきただろう。これで植物体重の変化からRGRの変化が推定できる。そして、RGRの変化についてはr=e^{a-bt}と仮定したので、定数abは、対数リンク関数を設定したglmで推定できる。

RGRからゴンペルツ曲線を推定する

まず、abが残ったバージョンでゴンペルツ関数を定義し、適当なパラメータでゴンペルツ曲線を描いてみよう。

gompertz <- function(t, c1, a, b){
  c1 * exp(-(exp(a-b*t)/b))
}

t <- seq(1, 20, 0.5)
a <- 0.3
b <- 0.25
c1 <- 10
plot(t, gompertz(t,c1, a, b), type = "l")

f:id:Rion778:20161020222308p:plain

次に重量を対数変換し、その変化量としてRGRの近似値を得てみよう。

g <- gompertz(t, c1, a, b)
r <- diff(log(g))/0.5
plot(t[-1], [r)

f:id:Rion778:20161020224446p:plain]

本当の植物の葉っぱなどを実際に測定をしてみる(重量測定は破壊が伴うので、例えば葉面積で代替する)と、RGRは実際にこのような減衰をすることが多い。

あとはglmでabを推定する。c_1についてはoptimize()で頑張って求める。

newt <- t[1:length(r)]+0.5
result <- glm(r ~ newt, family = gaussian(link="log"), start = c(1, 0))
result_coef <- coef(result)
plot(t, g)
f <- function(c1) (sum(gompertz(newt, c1, result_coef[1], -result_coef[2]) - g[-1]))^2
c1 <- optimize(f, interval = c(0, 100))$minimum
c1 <- 10
points(newt, gompertz(newt, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

するとこうだ。
f:id:Rion778:20161020225004p:plain
測定間隔の間はRGRが変化しないと仮定しているので、少しずれる。測定間隔が広いほどずれは大きくなる。

重量データに直接当てはめる

ちなみに、もし本当にゴンペルツ曲線のパラメータがほしいのであれば、重量の測定データに直接当てはめるべきで、このような回り道をすべきではない。
例えばoptim()を使うとこうなる。

f <- function(c) (sum((gompertz(t, c[1], c[2], c[3]) - g)^2))
res <- optim(c(10, 0.5, 1), f)$par
res # > [1] 9.9993114 0.3007776 0.2500798
plot(t, g)
points(t, gompertz(t, res[1], res[2], res[3]), col = "red", type = "l")

f:id:Rion778:20161020233012p:plain