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

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

$\frac{dw}{dt} = rw$

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

$w_t = w_0 e^{rt}$

しかし、植物は無限に大きくなるわけではないので、一般的には相対成長速度は時間とともに減少する。ここに例えば適当な定数 $a$ と $b$ を用意して、 $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}$ と仮定したので、定数 $a$ と $b$ は、対数リンク関数を設定したglmで推定できる。

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

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で $a$ と $b$ を推定する。 $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