Using a Bayesian Hierarchical Linear Mixing Model to Estimate Botanical Mixtures

Abstract

In grazing systems, estimating the dietary choices of animals is challenging but can be achieved using plant-wax markers, natural compounds that provide a signature of individual plants. If sufficiently distinct, these signatures can be used to characterize the makeup of a botanical mixture or diet. Bayesian hierarchical models for linear unmixing (BHLU) have been widely used for hyperspectral image analysis and geochemistry, but not diet mixtures. The aim of this study was to assess the efficiency of BHLU to estimate botanical mixtures. Plant-wax marker concentrations from eight forages found in Nebraska were used for simulating combinations of two, three, five and eight species. Also, actual forage mixtures were constructed in laboratory and evaluated. Analyses were performed using BHLU with 2 prior choices for forage proportions (uniform and Gaussian), 2 covariance structures (independent and correlated markers), stable isotope mixing models (SIMM), and nonnegative least squares (NNLS). Accounting for correlations between markers increased efficiency. Estimation error increased when Gaussian priors were used to model forage proportions. Performance of BHLU, SIMM, and NNLS was reduced with the more complex botanical mixtures and the limited number of markers. For simple mixtures, BHLU is a reliable alternative to NNLS for estimation of forage proportions. Supplementary materials accompanying this paper appear online.

Appendix

1.1 Forage Proportions

Although the Dirichlet distribution is not conjugate to the Gaussian likelihood, due to the choice of $\varvec{\upalpha }$ the posterior distributions resulting from combining (4) with (5), (6), or (7) have closed forms (Dobigeon et al. 2008) and can be efficiently sampled. Under (5) or (6) the full conditional distribution of the forage proportions is a truncated Gaussian on the simplex

$$\begin{aligned} \begin{aligned} f\big (\varvec{\uptheta }_{i,-p}|\varvec{\Sigma }, \mathbf {y}_{i}\big )&\propto |\varvec{\Psi }| \mathrm {exp}\bigg [-\dfrac{1}{2}\big (\varvec{\uptheta }_{i,-p} -\varvec{\upmu }_i \big )^{\prime } \varvec{\Psi }^{-1}\big (\varvec{\uptheta }_{i,-p} -\varvec{\upmu }_i \big )\bigg ]\\&\quad \times I\big (\varvec{\uptheta }_{i,-p}\in \mathbb {S}\big ),\\&\sim \phi _{\mathbb {S}}\big (\varvec{\upmu }_i,\; \varvec{\Psi }\big ), \end{aligned} \end{aligned}$$

(12)

where the mean vector is $\varvec{\upmu }_i=\varvec{\Psi }\Big [\big (\mathbf {M}_{-p}-\mathbf {m}_{p} \mathbf {1}^{\prime }\big )^{\prime }\varvec{\Sigma }^{-1}\big (\mathbf {y}_i- \mathbf {m}_{p}\big )\Big ]$ and the variance matrix is $\varvec{\Psi }=\Big [\big (\mathbf {M}_{-p}-\mathbf {m}_{p} \mathbf {1}^{\prime }\big )^{\prime }\varvec{\Sigma }^{-1}\big (\mathbf {M}_{-p}-\mathbf {m}_{p} \mathbf {1}^{\prime }\big )\Big ]^{-1}$. In addition, $\mathbf {M}_{-p}=\big [\mathbf {m}_1,\ldots , \mathbf {m}_{p-1}\big ]$, $\mathbf {m}_p$ is the $p\mathrm{th}$ column of M, and $\mathbf {1}=\big [1,\ldots ,1\big ]$ is a unit vector with $p-1$ elements. Under (7) the posterior distribution of forage proportions is also multivariate Gaussian truncated on a simplex,

$$\begin{aligned} \begin{aligned} f\big (\varvec{\uptheta }_{i,-p}|\varvec{\Sigma }, \mathbf {y}_{i}\big )&\propto |\varvec{\Psi }| \mathrm {exp}\bigg [-\dfrac{1}{2}\big (\varvec{\uptheta }_{i,-p} -\varvec{\widetilde{\upmu }}_i \big )^{\prime } \varvec{\widetilde{\Psi }}\big (\varvec{\uptheta }_{i,-p} -\varvec{\widetilde{\upmu }}_i \big )\bigg ]\\&\quad \times I\big (\varvec{\uptheta }_{i,-p}\in \mathbb {S}\big ),\\&\sim \phi _{\mathbb {S}}\big (\varvec{\widetilde{\upmu }}_{i},\; \varvec{\widetilde{\Psi }}\big ), \end{aligned} \end{aligned}$$

(13)

with updated mean vector $\overset{\sim }{\varvec{\upmu }}_i=\overset{\sim }{\varvec{\Psi }}\Big (\varvec{\Psi }^{-1} \varvec{\upmu }_i+\mathbf {V}^{-1}\mathbf {q}\Big )$, and updated variance matrix $\overset{\sim }{\varvec{\Psi }} =\big (\varvec{\Psi }^{-1}+\mathbf {V}^{-1}\big )^{-1}$.

1.2 Variance Matrix $\Sigma$

Given that the likelihood is Gaussian, and the choice of priors, in the case of $\varvec{\Sigma }=\mathrm {diag}\big (\sigma _{1}^{2},\ldots ,\sigma _{r}^{2}\big )$, the posterior distribution for the marker variances $\sigma _{k}^{2}$ is a scaled inverse Chi-squared as shown by Kazianka et al. (2011)

$$\begin{aligned} \begin{aligned} f\big (\sigma _{k}^{2}|\varvec{\uptheta }_{i},\lambda ,\nu ,\mathbf {y}_i\big )&\propto \big (\sigma _{k}^{2}\big )^{\frac{n+\nu }{2}-1}\\&\quad \times \mathrm {exp}\Bigg \lbrace -\dfrac{1}{2\sigma _{k}^{2}}\Bigg (\nu \lambda + \sum _{i=1}^{n} \big [y_{ik} - \mathbf {m}_k\varvec{\uptheta }_i\big ]^{2} \Bigg )\Bigg \rbrace ,\\&\sim \chi ^{-2}\Bigg (n+\nu ,\;\dfrac{\nu \lambda }{n+\nu }+\sum _{i=1}^{n} \dfrac{\big [y_{ik} - \mathbf {m}_k\varvec{\uptheta }_i\big ]^{2}}{n+\nu }\Bigg ), \end{aligned} \end{aligned}$$

(14)

where $\mathbf {m}_k$ is the $k\mathrm{th}$ row of the M matrix. In the case of priors as in (9), the posterior for $\varvec{\Sigma }$ is

$$\begin{aligned} \begin{aligned} f\big (\varvec{\Sigma }|\varvec{\uptheta }_i,\nu ,\varvec{\Xi },\mathbf {y}_i\big ) \propto \,&|\varvec{\Sigma }|^{-\frac{n+\nu +r+1}{2}}\mathrm {exp}\Bigg \lbrace -\dfrac{1}{2}\mathrm {tr} \Bigg [\bigg (\sum _{i=1}^{n}\big (\mathbf {y}_i-\mathbf {M}\varvec{\uptheta }_i\big ) \big (\mathbf {y}_i-\mathbf {M}\varvec{\uptheta }_i\big )^{\prime } {+}\, \varvec{\Xi } \bigg )\varvec{\Sigma }^{-1}\Bigg ]\Bigg \rbrace ,\\ \quad \sim&\; \mathrm {Inv.Wishart} \bigg (n+\nu ,\; \sum _{i=1}^{n}\big (\mathbf {y}_i-\mathbf {M}\varvec{\uptheta }_i\big ) \big (\mathbf {y}_i-\mathbf {M}\varvec{\uptheta }_i\big )^{\prime } + \varvec{\Xi } \bigg ). \end{aligned} \end{aligned}$$

(15)

1.3 Hyperparameters

The posterior distribution for the scale parameter $\lambda$ is

$$\begin{aligned} \begin{aligned} f\big (\lambda |\sigma _{k}^{2}, \mathbf {y}_i\big )&\propto \lambda ^{\frac{r\nu }{2}-1}\mathrm {exp}\bigg (-\dfrac{\lambda \nu }{2}\sum _{k=1}^{r} \dfrac{1}{\sigma _{k}^{2}}\bigg ),\\&\sim \mathrm {Gamma}\bigg (\dfrac{r\nu }{2},\;\dfrac{\nu }{2}\sum _{k=1}^{r} \dfrac{1}{\sigma _{k}^{2}} \bigg ). \end{aligned} \end{aligned}$$

(16)

Similarly, under the parameterization $\varvec{\Omega } = (\nu -r-1)\varvec{\Xi }$ and (11), the posterior distribution of the hyperparameter $\omega _{k}$ (Bouriga and Fern 2013) is

$$\begin{aligned} \begin{aligned} f\big (\omega _{k}|\nu ,\varvec{\Sigma },\mathbf {y}_i\big )&\propto \omega ^{\frac{r\nu }{2}}\mathrm {exp}\bigg [-\dfrac{r\nu }{2}\mathrm {tr}\big ( \varvec{\Sigma }_{k,k}^{-1}\big )\bigg ],\\&\sim \mathrm {Gamma}\bigg (\dfrac{r\nu }{2},\dfrac{\varvec{\Sigma }_{k,k}^{-1}}{2}\bigg ). \end{aligned} \end{aligned}$$

(17)