Bayesian Estimation for the GreenLab Plant Growth Model with Deterministic Organogenesis

Abstract

Plant growth modeling has attracted a lot of attention due to its potential applications. Many scientific disciplines are involved, and a lot of research effort and intensive computer methods were needed to understand better the complex mechanisms underlying plant evolution. Among the numerous challenges, one can cite mathematical modeling, parameterization, estimation and prediction. One of the most promising models that have been proposed in the literature is the GreenLab functional–structural plant growth model. In this study, we focus only on one of its versions, named GreenLab-1, particularly adapted to a certain class of plants with known organogenesis, such as sugar beet, maize, rapeseed and other crop plants. The parameters of the model are related to plant functioning, and the vector of observations consists of organ masses measured only once at a given observation time. Previous efforts for parameter estimation in GreenLab-1 include Kalman-type filters, stochastic variants of EM and/or ECM algorithms, and hybrid sequential importance sampling algorithms with Bayesian estimation only for the functional parameters of the model. In this paper, the first purely Bayesian approach for parameter estimation of the GreenLab-1 model is proposed. This approach has much more flexibility in handling complex structures, thus providing a useful tool for analyzing such types of models. In order to sample from the posterior distribution an MCMC algorithm is used and its implementation issues are also discussed. The performance of this method is illustrated on a simulated and a real dataset from the sugar beet plant, and a comparison is made with the MLE approach.

Appendix

The posterior density $\pi (\theta , q_{0:N}| y_{0:N})$ of the model under the following prior distributions

$$\begin{aligned} \alpha _b\sim & {} {\mathcal {N}} \left( m_{\alpha _b},\ \sigma _{\alpha _b}^2\right) \ {\mathbf {I}}_{(\alpha _b>1)},\ \alpha _p \sim {\mathcal {N}}\left( m_{\alpha _p},\ \sigma _{\alpha _p}^2\right) {\mathbf {I}}_{(\alpha _p>1)},\ \sigma _Q^2 \sim {{\mathcal {I}}}{{\mathcal {G}}}(\alpha , \beta ),\\ \mu ^{-1}\sim & {} {\mathcal {N}}_+\left( m_{\mu ^{-1}},\ \sigma _{\mu ^{-1}}^2\right) ,\ \ p_p \sim {\mathcal {N}}_+\left( m_{p_p},\ \sigma _{p_p}^2\right) ,\ \Sigma _Y \sim {{\mathcal {I}}}{{\mathcal {W}}}(\Psi ^{-1}, v). \end{aligned}$$

can be expressed in the form:

$$\begin{aligned}&\pi (\theta , q_{0:N}| y_{0:N})\\&\quad \propto \prod _{i=0}^{N-1}\frac{1}{\sigma _Q F_{i}\left( Q_{(i-T+1)^+:i}, E_{i}, \phi \right) }\cdot \exp \left\{ -\frac{\left( q_{i+1}-F_{i}\left( Q_{(i-T+1)^+:i}, E_{i}, \phi \right) \right) ^2}{2\sigma _Q^2F_{i}^2\left( Q_{(i-T+1)^+:i}, E_{i}, \phi \right) }\right\} \\&\qquad \prod _{i=0}^{N} |\Sigma _Y|^{-1/2}\exp \left\{ (y_i-G_{i}(Q_{i:(i+T-1)},p_{al}))'\Sigma _Y^{-1}(y_i-G_{i}(Q_{i:(i+T-1)},p_{al}))\right\} \\&\quad \exp \left\{ -\frac{(\alpha _b-m_{\alpha _b})^2}{2\sigma _{\alpha _b}^{{ 2}}}-\frac{(\alpha _p-m_{\alpha _p})^2}{2\sigma _{\alpha _p}^{{ 2}}} -\frac{(p_p-m_{p_p})^2}{2\sigma _{p_p}^{{ 2}}}-\frac{(\mu ^{-1}-m_{\mu ^{-1}})^2}{2\sigma ^{{ 2}}_{\mu ^{-1}}}\right\} \\&\quad I(\alpha _b>1) \cdot I(\alpha _p>1) \cdot I(p_p>0) \cdot I(\mu ^{-1}>0)\\&\quad {\prod _{i=0}^{N-1}}\cdot \frac{1}{(\sigma _Q^2)^{\alpha +1}}\exp \left\{ -\frac{\beta }{\sigma ^2_{Q}}\right\} |\Sigma _Y|^{(v+3)/2}\exp \left\{ -\frac{1}{2}tr(\Psi \Sigma _Y^{-1})\right\} . \end{aligned}$$

A hybrid Gibbs sampler is used to sample from the aforementioned posterior distribution. A hybrid Gibbs sampler is a Gibbs sampler (Geman and Geman 1984, Gelfand and Smith 1990) where at least one of the direct simulations from the full conditional distributions is replaced by a Metropolis–Hastings step (Metropolis et al. 1953, Hastings 1970).

Let $\pi (\theta , q_{1:N}|y_{0:N})$ be the density of the posterior (target) distribution with respect to some underlying measure $\mu$ and $\pi _i(\cdot |\cdot ), \ i= 1,\dots ,N+8$ be the full conditional distribution of the i component of ${\mathbf {w}} = (\theta , q_{1:N})$. The hybrid Gibbs sampler simulates an observation of ${\mathbf {w}} =(\theta , q_{1:N})$ into $N+8$ sequential steps by sampling from the full conditional distributions. In the case where direct sampling from $\pi _i(\cdot |\cdot )$ is very difficult or even infeasible a MH step is incorporated into Gibbs sampler. The MH algorithm with proposal distribution $p_i(\cdot |\cdot )$ generates a Markov chain with stationary distribution $\pi _i(\cdot |\cdot )$. Thus, the Hybrid Gibbs sampler generates a Markov Chain through the following transitions.

It starts from the state ${\mathbf {w}}^0$ in the support of $\pi$ and at time $t+1$ gives ${\mathbf {w}}^t$ for $i=1,\ldots , N+8$ proposes $z_i \sim p_i(z_i|w_i^t)$ and set $w_i^{t+1}=z_i$ with probability $\alpha (w_i^t, z_i)$; otherwise, set $w_i^{t+1}= w_i^{t}$ where

$a(w_i^{t}, z_i)=\min \left\{ 1, \frac{\pi _i(z_i|w_1^{t+1}\dots w_{i-1}^{t+1}, w_{i+1}^{t} \dots w_{N+8}^{t})}{\pi _i(w_i^{t}|x_1^{t+1}\dots w_{i-1}^{t+1}, w_{i+1}^{t} \dots x_{N+8}^{t})}\times \right. \left. \frac{p_i(w_i^{t}z_i)}{p_i(z_i| w_i^{t})} \right\}$

More specifically, the unobserved biomasses $q_1,q_2,\ldots , q_N$ and the parameters $\alpha _b$, $\alpha _p$, $p_p$ and $\mu ^{-1}$ are simulated using Metropolis–Hasting steps. The remaining parameters of the model, $\sigma _Q^2$ and $\Sigma _Y$, are updated by using the corresponding full conditional distributions, since they have a known form. The specific sweep strategy for each scan (iteration) of this hybrid Gibbs sampler is described below in detail.

• Step 1: $q_n,\ n=1,\ldots , N \ \hookrightarrow$ MH with target distribution the full conditional distribution of $q_n$

  $$\begin{aligned}&\pi _n(q_n|\cdots ) \\&\quad \propto \prod _{i=n+1}^{(n+T) \wedge N}\big (1-\exp \left\{ -\delta \, s_{i-1}^{act}(q_n)\right\} \big )^{-1} \nonumber \\&\qquad \times \exp \left\{ -\frac{1}{2{\sigma _Q^2}}\left( \left( \frac{q_n}{F_{n-1}(q_n)}-1\right) ^2+ \sum _{i=n+1}^{(n+T)\wedge N}\left( \frac{q_i}{F_{i-1}(q_n)}-1\right) ^2\right) \right\} \nonumber \\&\qquad \times \exp \left\{ -\frac{1}{2}\sum _{i=(n-T+1)^{+}}^{n}(y_i-G_i(q_n))^{\top }\Sigma _{Y}^{-1}(y_i-G_i(q_n))\right\} , \end{aligned}$$

  and proposal distribution the one resulting from the prior transition kernel of the hidden chain under the current parameters’ values,

  $$\begin{aligned} Q_{n+1}\sim {\mathcal {N}}_+(F_{n}\left( Q_{(n-T+1)^+:n}, E_{n}, \phi \right) , \sigma _Q^2 F_{n}^2\left( Q_{(n-T+1)^+:n}, E_{n}, \phi \right) ). \end{aligned}$$

• Step 2: $\mu ^{-1}\ \hookrightarrow$ Random Walk MH with target distribution the full conditional distribution of $\mu ^{-1}$

  $$\begin{aligned} \pi (\mu ^{-1}|\cdots )\propto & {} ({\mu ^{-1}})^N\exp \left\{ -\frac{1}{2 \sigma _{\mu 0}^2}\left( {\mu ^{-1}}-\frac{\mu _{\mu 0}}{\sigma _{\mu ^{-1}}^2}\right) ^2\right\} \end{aligned}$$

  where

  $$\begin{aligned} \sigma _{\mu 0}^2= & {} \left( \frac{1}{\sigma ^2_{\mu ^{-1}}}\sum _{i=0}^{N-1}\frac{q_{i+1}^2}{c(i)^2}+\frac{1}{\sigma ^2_{\mu ^{-1}}}\right) ^{-1}\\ \mu _{\mu 0}= & {} \sigma _{\mu 0}^2 \cdot \left( \frac{1}{\sigma ^2_{\mu ^{-1}}}\sum _{i=0}^{N-1}\frac{q_{i+1}}{c(i)}+\frac{\mu _{\mu ^{-1}}}{\sigma _{\mu ^{-1}}^2}\right) \\ c(i)= & {} E_{i} \mu S_p\left\{ 1-\exp \left( -\,\frac{k}{e S_p}M_{i,b}^{act}(q_{0:i}, \phi )\right) \right\} \end{aligned}$$

  and proposal distribution of ${\mathcal {N}}_+(\mu ^{-1}_{old}, \sigma ^{2}_{\mu ^{-1}MH})$.

• Step 3:     $\alpha _b \ \hookrightarrow$ Random Walk MH with target distribution the full conditional distribution of $\alpha _b$ and proposal distribution ${\mathcal {N}}(\alpha _{b-old}, \sigma ^{2}_{\alpha _bMH})\cdot {\mathbf {I}}_{(\alpha _b>1)}$.

• Step 4:     $\alpha _p \ \hookrightarrow$ Random Walk MH with target distribution the full conditional distribution of $\alpha _p$ and proposal distribution ${\mathcal {N}}(\alpha _{p-old}, \sigma ^{2}_{\alpha _pMH})\cdot {\mathbf {I}}_{(\alpha _p>1)}$.

• Step 5:     $p_p \ \hookrightarrow$ Random Walk MH with target distribution the full conditional distribution of $p_p$ and proposal distribution ${\mathcal {N}}_{+}(p_{b-old}, \sigma ^{2}_{p_pMH})$.

Since the functional parameters $\alpha _b, \alpha _p, p_p$ appear in all densities functions of both the observed and hidden data, their functional form is complex and it is omitted here. The remaining parameters of the model are simulated directly from their full conditional distributions,

• Step 6:     $\sigma _Q^2\sim {{\mathcal {I}}}{{\mathcal {G}}}\left( \alpha +\frac{N}{2},\ \beta + \frac{1}{2}\sum _{i=0}^{N-1}\frac{\left( q_{i+1}-F_{i}\left( Q_{(i-T+1)^+:i}, E_{i}, \phi \right) \right) ^2}{F_{i}\left( Q_{(i-T+1)^+:i}, E_{i}, \phi \right) ^2}\right)$.

• Step 7:     $\Sigma _Y\sim {{\mathcal {I}}}{{\mathcal {W}}}(A +\Psi ,\ v+N+1)$

  where $A= Y^* {Y^{*}}'$, and $y_i^* = y_i - G_{i}(Q_{i:(i+T-1)},p_{al})$.

Since all the steps (steps 2–5) which involve Random Walk MH proposals concern unidimensional quantities, all the variances of the proposal distributions are selected in such a way that the corresponding acceptance probabilities are as near as possible to the optimal one which, in the univariate case, is 0.44 (see, e.g., Roberts and Rosenthal 2001). The initial variances that we used in the simulated and the real data cases in steps 2–5 of the proposed algorithm, common to all priors, are given below for reference:

$$\begin{aligned} \sigma ^2_{MH(\mu _{-1})} =3, \quad \sigma ^2_{MH(a_b)} =0.06, \quad \sigma ^2_{MH(a_p)} = 0.08, \quad \sigma ^2_{MH(p_p)} =0.035. \end{aligned}$$