Embedding Positive Process Models into Lognormal Bayesian State Space Frameworks using Moment Matching

State Space Models, sometimes referred to as Hidden Markov Models (HMMs), are a broad class of models used to track the states of a model (some unobserved process $X_{1:T}$) through a set of observations, $Y_{1:T}$ (Durbin and Koopman,, 2012; Petris et al.,, 2009). A state space model has three components – the state component, the observation component, and additional parameters. The state component consists of an (unobserved) Markov process, $X_{t}$, being moved through time by an evolution function (or process function) $f(X_{t}\lvert X_{t-1},\Theta)$. The implication of $X_{1:T}$ following a Markov process is that the past and the future are independent conditional on the state at the current time, $X_{t}$ (Shumway and Stoffer,, 2011). The observation component consists of noisy observations of the latent process, $Y_{1:T}$, that are governed by an observation density function, $g(Y_{t}\lvert X_{t},\Theta)$. These observations are assumed to be independent of one another conditional on the latent states $X_{1:T}$. The additional parameter component, $\Theta$, contains the parameters that govern the evolution function and observation function.

The Gompertz (Gompertz,, 1825) and Moran-Ricker (Ricker,, 1954) SSMs are simple discrete density-dependent state space models with lognormally distributed process and measurement error (Dennis et al.,, 2006). The Gompertz and Moran-Ricker SSMs are popular choices for lognormal state space models because they can be easily transformed into Normal Dynamic Models (NDMs), and can then take advantage of a suite of well studied fitting methods including Kalman filtering (Kalman,, 1960), extended Kalman filtering (Julier and Uhlmann,, 1997), and Gibbs sampling (Geman and Geman,, 1984; Carter and Kohn,, 1994), easing computational difficulties associated with fitting SSMs. The latent process models for the Gompertz model can be written as:

$\displaystyle X_{t}$

$\displaystyle=X_{t-1}\exp(a+b\log(X_{t-1})+\epsilon_{t}),\text{ }\epsilon_{t}\sim\mathcal{N}(0,\phi)$

(5)

The latent process model for the Moran-Ricker model can be written as:

$\displaystyle X_{t}$

$\displaystyle=X_{t-1}\exp(a+bX_{t-1}+\epsilon_{t}),\text{ }\epsilon_{t}\sim\mathcal{N}(0,\phi)$

(6)

Letting $A=\exp(a)$, we can rewrite Equations 5 and 6, the Gompertz and Moran-Ricker process equations (respectively), in the form:

	$\displaystyle X_{t}$	$\displaystyle=A(X_{t-1})^{b+1}\exp(\epsilon_{t}),\text{ }\epsilon_{t}\sim\mathcal{N}(0,\phi)$		(7)
	$\displaystyle X_{t}$	$\displaystyle=AX_{t-1}\exp(bX_{t-1})\exp(\epsilon_{t}),\text{ }\epsilon_{t}\sim\mathcal{N}(0,\phi)$		(8)

More generally, we can think of both the Gompertz and Moran-Ricker process models as belonging to a class of lognormally distributed process models that have the form:

$\displaystyle X_{t}=f^{*}(X_{t-1}\lvert\Theta)\exp(\epsilon_{t}),\text{ }\epsilon_{t}\sim\mathcal{N}(0,\phi),\text{ }f^{*}(X_{t-1}\lvert\Theta)>0,$

(9)

where $f^{*}(X_{t-1}\lvert\theta)$ is the model of choice for the non-negative physical process, and the process error is governed by a multiplicative zero mean lognormal distribution. Key statistical properties (conditional mean, conditional variance, and conditional median) of this class of process models are:

	$\displaystyle\mathbb{E}[X_{t}\lvert X_{t-1}]$	$\displaystyle=f^{*}(X_{t-1}\lvert\Theta)\exp((2\phi)^{-1}),$		(10)
	$\displaystyle\mathbb{V}[X_{t}\lvert X_{t-1}]$	$\displaystyle=f^{*}(X_{t-1}\lvert\Theta)^{2}\exp(\phi^{-1}),\Big{(}\exp(\phi^{-1})-1\Big{)}$		(11)
	$\displaystyle\mathscr{M}[X_{t}\lvert X_{t-1}]$	$\displaystyle=f^{*}(X_{t-1}\lvert\Theta).$		(12)

When looking at the conditional expected value, we see that this class of models is biased in terms of mean process evolution, but unbiased in terms of median process evolution. Further, the variance of the latent state is controlled by the value of the process model $f^{*}(\cdot)$, with larger values of the process assumed to have larger variation than smaller ones. Therefore, in the class of process models that describe the Gompertz and Moran-Ricker models, the density dependent relationship of the latent process variance is assumed to scale quadratically with the predicted value of the process model.

The generalization of the Gompertz and Moran-Ricker process models also holds true for the observation model. For example, for continuous responses $Y$, the assumed relationship between an arbitrary observation $Y_{t}$ and the corresponding latent state $X_{t}$ for the Gompertz model may be given by:

$\displaystyle Y_{t}=X_{t}\exp(\epsilon_{obs,t}),\text{ }\epsilon_{obs,t}\sim\mathcal{N}(0,\tau)$

(13)

Broadly, the Gompertz and Moran-Ricker observation models belong to the same class of lognormally distributed observation models in Equation 13 that have the form:

$\displaystyle Y_{t}=g^{*}(X_{t}\lvert\Theta)\exp(\epsilon_{obs,t}),\text{ }\epsilon_{obs,t}\sim\mathcal{N}(0,\tau),\text{ }g^{*}(X_{t}\lvert\Theta)>0,$

(14)

where $g^{*}(X_{t}\lvert\Theta)$ is now interpreted as the model used to link our observations to their respective latent states, and the measurement error is dictated by a multiplicative zero mean lognormal distribution. The conditional mean, conditional variance, and conditional median of this class of observation models are given by:

	$\displaystyle\mathbb{E}[Y_{t}\lvert X_{t}]$	$\displaystyle=g^{*}(X_{t}\lvert\Theta)\exp((2\tau)^{-1}),$
	$\displaystyle\mathbb{V}[Y_{t}\lvert X_{t}]$	$\displaystyle=g^{*}(X_{t}\lvert\Theta)^{2}\exp(\tau^{-1})\Big{(}\exp(\tau^{-1})-1\Big{)},$
	$\displaystyle\mathscr{M}[Y_{t}\lvert X_{t}]$	$\displaystyle=g^{*}(X_{t}\lvert\Theta).$

Importantly, while this class of observation models assumes that observations are unbiased in log space, they assume there is systematic observation bias in their measurement process, which is given by $B(Y_{t}\lvert X_{t})=g^{*}(X_{t}\lvert\Theta)(1-\exp((2\tau)^{-1}))$.

Though we have only examined the Gompertz and Moran-Ricker models when discussing the class of lognormal process and observation models presented here, there are many examples of lognormal modeling frameworks in the animal population modeling literature (see Buckland et al.,, 2004; Knape et al.,, 2011, for examples) and fisheries modeling literature (see Maunder et al.,, 2015; Mäntyniemi et al.,, 2015, for examples) that fall into the classes we discuss here. Including a term to correct the bias is a common method used in applications to rectify assumptions of biased process evolution and biased observations (Knape et al.,, 2011; Maunder et al.,, 2015; Mäntyniemi et al.,, 2015) but this bias correction term changes the variance structure, as the lognormal variance is a function of both $\mu$ and $\sigma^{2}$.

Instead of assuming an unbiased median process evolution and a density dependent variance structure, a modeler developing a lognormal SSM for their application may want to specify the mean evolution and the variance structure of the stochastic lognormal process model. That is, we desire a lognormally distributed stochastic evolution function such that $\mathbb{E}[X_{t}\lvert X_{t-1}]=f^{*}(X_{t-1}\lvert\Theta)$ and $\mathbb{V}[X_{t}\lvert X_{t-1}]=\phi_{t}^{-1}$. We can create a process model with these characteristics by using a moment matching transformation on the mean and precision, specifically:

	$\displaystyle\mu^{*}_{t}$	$\displaystyle=\log\Big{(}\frac{f^{*}(X_{t-1}\lvert\Theta)^{2}}{\sqrt{f^{*}(X_{t-1}\lvert\Theta)^{2}+\phi_{t}^{-1}}}\Big{)},$		(15)
	$\displaystyle\phi^{*}_{t}$	$\displaystyle=\log\Big{(}1+(f^{*}(X_{t-1}\lvert\Theta)^{2}\phi_{t})^{-1}\Big{)}^{-1}.$		(16)

This approach provides a stochastic lognormal process model with the desired properties (see Appendix 6.1 for derivation):

	$\displaystyle X_{t}\lvert X_{t-1}$	$\displaystyle\sim\text{Lognormal}(\mu^{*}_{t},\phi^{*}_{t}),$
	$\displaystyle\mathbb{E}[X_{t}\lvert X_{t-1}]$	$\displaystyle=f^{*}(X_{t-1}\lvert\Theta),$
	$\displaystyle\mathbb{V}[X_{t}\lvert X_{t-1}]$	$\displaystyle=\phi_{t}^{-1}.$

Thus, there is a framework for embedding the process model $f^{*}(X_{t-1}\lvert\Theta)$ such that the temporal evolution is unbiased, and that allows for a flexible way to model the variance of the process through time. For example, we can choose $\phi_{t}=\phi$ to model constant variance through time, $\phi_{t}=\phi/f^{*}(X_{t-1}\lvert\Theta)^{2}$ to model density dependent variance, or $\phi_{t}=\phi\exp(-t^{-1})$ to model variance that dissipates over time.

The lognormal moment matching approach can be similarly applied to the observation model to produce a lognormally distributed observation density with the properties $\mathbb{E}[Y_{t}\lvert X_{t}]=g^{*}(X_{t}\lvert\Theta)$, $\mathbb{V}[Y_{t}\lvert X_{t}]=\tau_{t}$ by using the transformation:

	$\displaystyle\mu^{*}_{t,obs}$	$\displaystyle=\log\Big{(}\frac{g^{*}(X_{t}\lvert\Theta)^{2}}{\sqrt{g^{*}(X_{t}\lvert\Theta)^{2}+\tau_{t}^{-1}}}\Big{)},$		(17)
	$\displaystyle\tau^{*}_{t}$	$\displaystyle=\log\Big{(}1+(g^{*}(X_{t}\lvert\Theta)^{2}\tau_{t})^{-1}\Big{)}^{-1}.$		(18)

With the forms of both the process model and observation model fully specified in this fashion, it is possible to write a likelihood equation for a model that includes both components, that we call the Lognormal Moment Matching model (LNM3). Suppose that we have observations $\mathbf{Y}$ at a subset of time points $I\subset\{1,\dots,T\}$. Then the likelihood for the LNM3 is given by

	$\displaystyle\mathscr{L}(X_{1:T},\Theta\lvert\mathbf{Y})=$
	$\displaystyle\prod_{t=1}^{T}\text{log}\mathcal{N}\left(\log\left(\frac{f^{*}(X_{t-1}\lvert\Theta)^{2}}{\sqrt{f^{*}(X_{t-1}\lvert\Theta)^{2}+\phi_{t}^{-1}}}\right),\log\left(1+(f^{*}(X_{t-1}\lvert\Theta)^{2}\phi_{t})^{-1}\right)^{-1}\right)$
	$\displaystyle\times\prod_{i\in I}\text{log}\mathcal{N}\left(\log\left(\frac{g^{*}(X_{i}\lvert\Theta)^{2}}{\sqrt{g^{*}(X_{i}\lvert\Theta)^{2}+\tau_{i}^{-1}}}\right),\log\left(1+(g^{*}(X_{i}\lvert\Theta)^{2}\tau_{i})^{-1}\right)^{-1}\right)$		(19)

The moment matching method can be used to embed any positive process or observation model and any variance structure into a stochastic lognormal model, and thus helps to maintain biophysical realism when modeling positive processes. Fitting these models as state space models further allows them to easily handle missing data and partition between process and measurement error. Thus the LNM3 approach provides a framework that is flexible, biophysically realistic, and statistically coherent.

Our analysis focused on comparing and contrasting the LNM3 approach to the Gompertz and Moran-Ricker approaches, two common lognormal SSMs. We estimate the latent states, precisions, and model parameters for six different models presented here using Bayesian state space models (SSMs) (Hamilton,, 1994; Durbin and Koopman,, 2012; Petris et al.,, 2009). The six different models include the Gompertz SSM, the Moran-Ricker SSM, and four different SSM formulations using the Lognormal Moment Matching (LNM3) method presented in Equations 15 – 18. Model parameters, latent states, and precisions were estimated using Markov Chain Monte Carlo (MCMC) (Robert and Casella,, 2005). MCMC is an estimation method that generates samples of parameters from their posterior distributions using Markov chains. Parameter uncertainty can be quantified using posterior samples from these Markov chains. All models were fit using JAGS (Plummer,, 2003) using the rjags package (Plummer,, 2019) in R version 4.1.0 (R Core Team,, 2016).

We designed a simulation study to assess the forecasting performance of the six models presented here, both under the cases where they are the true generating models and where the models are mis-specified. Our three primary objectives for the simulation study were: 1) assess forecasting performance of each model when it is the true generating model and when it is misspecified, when observation precision is being estimated 2) assess forecasting performance of each model when it is the true generating model and when it is misspecified, with fixed observation precision; 3) analyze the estimability of precisions for the models considered in this manuscript.

To assess the forecast performance of our models, we used proper scoring rules (Gneiting and Raftery,, 2007). Broadly, proper scoring rules use information about the predictive distribution coupled with observations to assign a measure of agreement of the forecast and the observations (Krüger et al.,, 2021). Specifically, Gneiting and Raftery, (2007) define a scoring rule to be proper if the expected value of the score is maximized by a draw from the true forecast distribution, and show that both the Continuous Ranked Probability Score (CRPS; Matheson and Winkler,, 1976) and the Ignorance (IGN; Good,, 1952; Roulston and Smith,, 2002) score are proper scoring rules. Given a probability density function $f(\cdot)$ and corresponding cumulative distribution function $F(\cdot)$ for our forecast, and an observation $y$, the IGN and CRPS may be written:

	$\displaystyle\text{IGN}(y)$	$\displaystyle=-\log(f(y));$		(41)
	$\displaystyle\text{CRPS}(y)$	$\displaystyle=\int_{\mathbb{R}}(F(z)-\mathbbm{1}_{z\geq y})^{2}dz.$		(42)

For our simulation studies, we consider both the CRPS and IGN scores, so that we have one scoring rule defined in terms of the probability density function (IGN) and one scoring rule defined in terms of the cumulative distribution function (CRPS). We use the IGN and CRPS scores to quantitatively compare forecasts, with lower scores within a scoring rule indicating a better performance. CRPS has support over the positive real line, $[0,\infty)$, while IGN can take values between $[-\log(f(y^{*})),\infty)$, where $y^{*}=\text{argmax}_{y\in\mathbb{R}}f(y)$.

To analyze the estimability of precision parameters in each of our six models, we used the empirical coverage rate of the 95% highest posterior density credible intervals for the precision parameters. Auger-Méthé et al., (2016) show that SSMs can have difficulty recovering the process and observation precisions even when the models are linear and Gaussian. To perform a thorough analysis on the estimability of the precision parameters, we fit each model under two different scenarios. In the first scenario, we fixed the observation precision and estimated the process precision for each model. In the second scenario, we estimated both the observation precision and the process precision for each model. This approach allowed us to assess the estimability by looking at the increase in empirical coverage rate for the process precision when the observation precision is fixed compared to when the observation precision is estimated. By choosing to use coverage and proper scoring rules to quantify our model performance we follow a common approach used in the literature – using frequentist concepts to assess Bayesian models (Box,, 1980; Rubin,, 1984; Little,, 2006, 2012).

To quantify our three objectives, we performed the simulation study as follows: for each of the six models, thirty different synthetic datasets of length 575 were generated by simulating from the underlying process model and observation model, with parameter values taken from Table 2. Parameter values for each model were chosen so that the systems had similar mean dynamics for each generating model. Each synthetic dataset was fit to each of the six models discussed here. Models were initially fit with the first 365 days of data, and then forecasts were computed for days 366 – 372, a seven day forecast horizon. The average IGN and CRPS for the seven day forecasts horizons were computed using the logs_sample and crps_sample function from the scoringRules R package (Jordan et al.,, 2017). We also saved the highest posterior density credible intervals for the estimated precisions. We then re-fit the models with the first 372 days of data, and forecasts were computed for days 373–379. This process was repeated until the synthetic data was exhausted. Overall, each individual synthetic data set was fit a total of 12 times – once by each model with observation precision known and once by each model with observation precision estimated, for a total of 2,160 simulations.

This simulation study design helps us to assess our three objectives. Our first two objectives were to assess the forecasting performance of each model under model mis-specification; both in the case of estimated observation precision and in the case of fixed observation precision. Our simulation study had each model fit each dataset twice: once with observation precision fixed and once with observation precision estimated. The forecasts for each precision scenario were evaluated using the CRPS and IGN scores. By using the same synthetic datasets, we were able to isolate the differences in forecast performance for each precision scenario and quantify the difference using paired t-tests with a Holm-Bonferroni adjustment (Holm,, 1979). With thirty different synthetic datasets generated by each model and each synthetic dataset evaluated over thirty forecast horizons by each of the six models, we were able to get a comprehensive picture of how each model performed under model mis-specification for each of the precision scenarios. Our third objective was to analyze the estimability of precisions for the models considered in this manuscript. We used differences between the scores of the forecasts and coverage rates of the precision credible intervals to determine the severity of underlying precision estimation problems. In the absence of estimation problems, we expect that the empirical coverage rates of the precision parameters should achieve close to the nominal coverage rate of 95%. CRPS and IGN are proper scoring rules, and their expected scores are maximized by draws from the true forecast distribution (Gneiting and Raftery,, 2007). Thus we expect that across a large number of simulations, the true generating model should score best in terms of CRPS and IGN in the absence of estimation problems.

We examined the differences in predictive performance between two different LN-SSM formulations by modeling Leaf Area Index (LAI) at University of Notre Dame Environmental Research Center (UNDE), a National Ecological Observatory Network (NEON) site. LAI is the total surface of leaves per area of ground and is a metric of photosynthetic capacity of an ecosystem. We obtained estimates LAI at UNDE from Oak Ridge National Laboratory Distributed Active Archive Center using their fixed subsets feature (DAAC,, 2018). The dataset contains estimates of LAI at four day intervals, computed using reflectance estimates the from Moderate Resolution Imaging Spectroradiometers (MODIS) satellite as inputs to a phenological model that predicts LAI (Yang et al.,, 2006). The two statistical models that we considered were a biased LN-SSM, based on Equations 9 and 14, and a moment matching LN-SSM, based on Equations 15, 16, 17, and 18. Our analysis was designed with two primary questions in mind: 1) can we construct informed out-of-sample predictions for LAI for long horizons (e.g., multiple seasons) by embedding a process-based ecosystem model into a lognormal state space framework?; 2) do the moment matching and biased formulations show differences in out-of-sample predictive performance for long horizons, when measured by CRPS and IGN?

To fit the LAI data, we used a reduced version of the DALEC2 process-based ecosystem model (Bloom and Williams,, 2015). In DALEC2, the LAI is modeled as a function of the carbon stored in foliage. Rather than fitting the full model (a six dimensional dynamical system with 23 process parameters) we only considered the interaction between foliage carbon ($C_{f}$) and the labile carbon ($C_{lab}$). This reduced the model to a two dimensional dynamical system with eight process parameters. The form of the reduced process model is given as:

	$\displaystyle\mathbf{C}^{(t)}=M_{t}\mathbf{C}^{(t-1)}+p_{t},\text{ where}$		(43)
	$\displaystyle M_{t}=\begin{bmatrix}1-\Phi_{f}^{(t)}&\Phi_{o}^{(t)}\\ 0&1-\Phi_{o}^{(t)}\end{bmatrix},\mathbf{C}^{(t-1)}=\begin{bmatrix}C_{f}^{(t-1)}\\ C_{lab}^{(t-1)}\end{bmatrix},p_{t}=\begin{bmatrix}G(\mathbf{D^{(t)}},c_{lma},c_{eff})f_{f}\\ G(\mathbf{D^{(t)}},c_{lma},c_{eff})f_{lab}\end{bmatrix}$
	$\displaystyle\text{LAI}^{(t)}=\frac{C_{f}^{(t)}}{c_{lma}}$		(44)

	$\displaystyle\Phi_{f}(t,d_{f},c_{r},c_{lf})=\sqrt{\frac{2}{\pi}}\cdot\frac{-\log(1-c_{lf})}{c_{rf}}\exp\left(-\left(\sin\left(\frac{t-d_{f}+\psi_{f}}{s}\right)\frac{\sqrt{2}s}{c_{rf}}\right)^{2}\right),$		(45)
	$\displaystyle\Phi_{o}(t,d_{o},c_{ro})=\sqrt{\frac{2}{\pi}}\left(\frac{6.9088}{c_{ro}}\right)\exp\left(-\left(\sin\left(\frac{t-d_{o}+.6245c_{ro}}{s}\right)\frac{\sqrt{2}s}{c_{ro}}\right)^{2}\right),$		(46)

where $s=365.25/\pi$ and $\psi_{f}=-\sqrt{W_{0}\left(\left(2\pi\log(1-c_{lf})^{2}\right)^{-1}\right)}/\sqrt{2}$, where $W_{0}$ is the principal branch of the Lambert W function (Lambert,, 1758), and $G(\mathbf{D^{(t)}},c_{lma},c_{eff})$ is the output of the Aggregated Canopy Model (ACM) for gross photosynthetic production (Williams et al.,, 1997). $\mathbf{D}^{(t)}$ represents meteorological driver variables for day $t$ that include: daily minimum and maximum temperatures (^∘C), daily incoming shortwave radiation (gCm^-1), and atmospheric carbon (CO₂ ppm). Minimum temperatures, maximum temperatures, and daily shortwave radiation were obtained from the National Ecological Observatory Network (NEON; National Ecological Observatory Network,, 2020; National Ecological Observatory Network , 2022b, NEON; National Ecological Observatory Network , 2022a, NEON). We imputed any missing NEON observations using a piece-wise linear interpolation. We took monthly measurements of atmospheric carbon from the Scripps Project (Keeling et al.,, 2005), and interpolated daily measurements by assigning the monthly values to each day within the month.

To fit the reduced DALEC2 as a LN-SSM, we embedded the process-based model (Equation 43) as the state component and the LAI equation (Equation 44) as the observation component. To model the variance structure of the process evolution function and observation function, we used a density dependent variance. We chose a density dependent variance structure because we expect more process variation and measurement error for foliage carbon when it is large during the early spring and summer months, and less when it is low in the winter months. We considered two different LN-SSM formulations with density dependent variance structures. The first model, the biased model, embeds the process and observation components in a biased manner using Equation 9 and Equation 14. The second model that we consider embeds the process and observation components in an unbiased manner using a moment matching approach. The process evolution function and the observation function for the moment matching model are given by:

	$\displaystyle\mathbf{C}^{(t)}\lvert\mathbf{C}^{(t-1)}\sim\text{MV}\log\mathcal{N}\left(\log\left(M_{t}\mathbf{C}^{(t-1)}+p^{(t)}\right)-.5\log(\mathbf{11}^{T}+\Omega_{m})\mathbf{1},\log(\mathbf{11}^{T}+\Omega_{m})\right),$
	$\displaystyle\text{LAI}^{(i)}\lvert\mathbf{C}^{(i)}\sim\log\mathcal{N}\left(\log\left(\frac{C_{f}^{(i)}}{c_{lma}}\right)-.5\log(1+\tau_{m}^{-1}),\log(1+\tau_{m}^{-1})\right),$

where $M_{t},p_{t}$ and $\mathbf{C^{(t)}}$ take on the same values from Equation 43, $\mathbf{1}$ represents a $2\times 1$ vector of ones, $\Omega_{m}=\text{Diag}(\omega^{2}_{f,m},\omega^{2}_{lab,m})$, and the $\log$ in the process evolution variance represents taking the component-wise logarithm of each element of the matrix/vector.

The process evolution function and observation function for the biased model are given by:

	$\displaystyle\mathbf{C}^{(t)}\lvert\mathbf{C}^{(t-1)}\sim\text{MV}\log\mathcal{N}\left(\log\left(M_{t}\mathbf{C}^{(t-1)}+p^{(t)}\right),\Omega_{b}\right),$
	$\displaystyle\text{LAI}^{(i)}\lvert\mathbf{C}^{(i)}\sim\log\mathcal{N}\left(\log\left(\frac{C_{f}^{(i)}}{c_{lma}}\right),\tau_{b}^{-1}\right),$

where $\Omega_{b}=\text{Diag}(\omega^{2}_{f,b},\omega^{2}_{lab,b})$, and the $\log$ once again represents a component-wise logarithm. The subscripts $b$ and $m$ are used to distinguish between parameters of the biased and moment matching models, respectively. Additional information on expected values and variances for the two different DALEC2 LN-SSM formulations can be found in Table 6 in the Appendix.

The equations above that we use to define our LN-SSM also highlight an interesting difficulty of our analysis: we are fitting a two state dynamical system, but we only have observations for one of the states, ($C_{f}$). Although we are primarily interested in predicting LAI, which is a function of $C_{f}$, we still need to estimate labile carbon so that we can model it’s contribution to foliage during the period of leaf regrowth. This also highlights an advantage of our state space modeling approach: the uncertainty from our lack of labile carbon measurements is propagated through time, giving us a more complete picture of the uncertainty in the foliage carbon.

We treated two model parameters as fixed: the Leaf Mass per Area ($c_{lma}$) and the density dependent observation precision parameter ($\tau$). We fixed $c_{lma}$ to avoid identifiability issues with the canopy efficiency parameter ($c_{eff}$), as the two parameters appear exclusively together in the ACM (Williams et al.,, 1997). For both models, we fixed $c_{lma}$ to a value of 75, based on empirical results from Serbin et al., (2019) that were calibrated specifically at UNDE. We estimated the density dependent observation precision parameter using historical data from MODIS, which reports both the standard deviation and the mean of the LAI estimates. For the moment matching model, we took the median value of the means divided by the standard deviations. This gave us a value of $\hat{\tau}_{m}\approx 4$. For the standard lognormal SSM formulation, we used the form of the variance for the lognormal distribution (Equation 10) to obtain $\hat{\tau}_{m}\approx 4.18$. For the DALEC2 process parameters, we used uniform prior distributions over the range of acceptable values taken from Bloom and Williams, (2015). For the density dependent process variance components, we used a $\text{Uniform}(0,1)$ distribution as the prior. We chose this because the $\omega$ parameters for the biased and moment matching parameters have the interpretation that they roughly control the average proportion of process error at each time point. All prior distributions for parameters were identical for the moment matching model and the biased model. A full table detailing prior distributions and fixed parameters can be found in Appendix 6.3.

We fit both of the models using particle MCMC (pMCMC; Andrieu et al.,, 2010). We ran each model for a total of 881 days, from September 5th, 2019 to February 2nd, 2022. We used four day MODIS LAI measurements from September 5th, 2019 to September 6th, 2021 to fit each model, and we used the remaining MODIS LAI measurements from September 7th, 2021 to February 2nd, 2022 to assess out of sample prediction. To assess out of sample predictions, we used the pMCMC samples of the latent states to generate samples from the posterior predictive distribution for the observations. We then used these samples to validate against the out of sample MODIS LAI measurements using CRPS and IGN using the scoringRules package (Jordan et al.,, 2017) in the R programming language (R Core Team,, 2016). By doing this, we are scoring on $Y\lvert X$ rather than directly on the observation $Y$, and acknowledge that the LAI observations that we are using for validation have measurement error (Ferro,, 2017; Bessac and Naveau,, 2021). We implemented pMCMC using the R package pomp (King et al.,, 2016). We ran each model for a total of 100,000 iterations, with a burn in of 50,000 iterations, 500 particles, and an adaptive multivariate normal proposal distribution (Andrieu and Thoms,, 2008; Rosenthal,, 2009) that began using a scaled empirical covariance matrix after 1,000 samples are accepted. To obtain initial parameter estimates that start in a region of high posterior density, we used a Gaussian process surrogate model optimization (Gramacy,, 2020) using TRBO (Eriksson et al.,, 2019), with the particle filter marginal log-likelihood as the objective function.

LAI is ubiquitous for forecasting changes in carbon stored in different components of a terrestrial ecosystem under different projections of climate. For our analysis, we chose to predict over multiple seasons, in an attempt to emulate a forecasting scenario of LAI response to predicted changes in climate, such as a particularly hot or cold winter. Our analysis also serves as an efficacy test for predictive modeling of LAI using a state space framework. While much work has been done on LAI prediction using ecosystem process-based models (see Mahowald et al.,, 2016; Ercanli et al.,, 2018, for examples), to our knowledge there has been little work done on predicting LAI by embedding a mechanistic process-based ecosystem model as the process component of a statistical state space model.

The first objective of our simulation study was to assess the forecasting performance of each model under model mis-specification when observation precision is estimated. We expected that the true generating model would perform best on average for its thirty synthetic datasets, using the average IGN and CRPS scores over each seven day forecast horizon. For the case where both the process precision, $\phi$, and the observation precision, $\tau$, were estimated, the Gompertz model (Gomp in Table 3) had the best forecasting performance among the density dependent models (MR, Gomp, LMRD, LGD) for both CRPS and IGN. The unbiased moment matching analog to the Gompertz, LGD, also had a strong performance, scoring second highest for IGN on all four density dependent variance models, and scoring second highest for CRPS on three out of the four. For the constant variance models (LMRC, LGC), the constant variance Moran-Ricker model (LMRC) had the best forecasting performance for both CRPS and IGN.

The second objective of our simulation study was to assess the forecasting performance of each model under model mis-specification when observation precision was fixed. For the simulations where the observation precision, $\tau$, was fixed and the process precision, $\phi$, was estimated, we found higher consistency between our scoring rules and the generating model. For the simulations that estimated both precision parameters (Table 3), the Gomp and LGD models consistently outperformed the other density dependent variance models for both scoring rules. In contrast, the simulations where $\tau$ was fixed and $\phi$ was estimated (Table 4) showed a more equal representation among the density dependent models. For the four density dependent generating models, the Moran-Ricker and its unbiased moment matching counterpart scored the best for CRPS, with each model producing the lowest score twice. When measuring performance with IGN, the Moran-Ricker and Gompertz models each scored lowest twice. The density dependent models all scored similarly, agnostic to the choice of true generating model, both for CRPS (max difference $\approx.0052$) and for IGN (max difference $\approx.006$).

Our third objective was to analyze the estimability of precisions for each of the six models we used. Our first investigation for the estimability of the precision parameters was to quantify the differences in IGN and CRPS scores between the simulations where the observation, $\tau$ was fixed and the simulations where $\tau$ was estimated. This investigation is also intimately related to the performance of the models under model mis-specification, and therefore provides insight for all three of our objectives. To evaluate the impacts statistically, we used a paired t-test. To do this, we took the average IGN and CRPS scores over each seven day forecast horizon for the two different scenarios, and treating them as ”before” and ”after”. We justify this by noting that each precision scenario was fit using identical synthetic datasets, with the only difference being whether $\tau$ was estimated or not. Unsurprisingly, we found that fixing the observation precision helped to improve forecasting performance for nearly all models. The LGC and LMRC models performed significantly better in terms of CRPS when $\tau$ was fixed (LGC: $p=4.7\text{e-}05$; LMRC: $p=.0032$), but did not perform significantly better in terms of IGN (LGC: $p=.68$; LMRC: $p=.36$). The Moran-Ricker, LMRD, and LGD models all had statistically significant decreases in IGN (MR: $p<2\text{e-}16$; LMRD: $p=2.1\text{e-}09$; LGD: $2.8\text{e-}05$) as well as CRPS (MR: $p<2\text{e-}16$; LMRD: $p<2\text{e-}16$; LGD: $1.4\text{e-}05$) when $\tau$ was fixed. Fixing the observation precision did not significantly impact the performance of the Gompertz model for IGN ($p=.84$) or CRPS ($p=.82$).

Our second investigation for the estimability of the precision parameters was to use the empirical coverage rate of the HPD intervals. The empirical coverage of the highest posterior density (HPD) credible intervals for $\phi$ unanimously increased for the simulations where $\tau$ was fixed, and all six models had empirical coverage that falls close to the nominal rate. For the simulations where both $\phi$ and $\tau$ are estimated, Table 5 shows that the Gompertz model and LGD models produce the best empirical coverage for the precision parameters. This is likely to be related to their excellent performance in these simulations, where other models struggled to consistently produce precision estimates that contained the ground truth in their 95% HPD intervals. The empirical undercoverage of the Moran-Ricker (MR) model for both $\phi$ and $\tau$ may also explain its poor performance in the forecast results from Table 3, where it came in last place for CRPS for all six generating models and last place for IGN for three out of six generating models, including the case where it was itself the true generating model.

We found that both the moment matching LN-SSM and the biased LN-SSM produced predictive distributions that captured the dynamics of both the in sample and the out of sample LAI observations. Both models showed similar fits for the in sample LAI observations when looking at medians (Figure 3) and 90% highest posterior density intervals (Figure 3). This was not surprising to us, as both models used an identical process model and prior distributions, and only differed slightly in the formulation of process evolution and observation functions. For the out of sample LAI predictive distributions, the models behaved differently. The biased model (Figure 3, top panel) had lower variance at the start of the predictive horizon, and then began to tail off at the end of the horizon. The moment matching model (Figure 3, bottom panel) had larger predictive variance at first, but then leveled off and accumulated slowly at the end of the horizon. The median predicted values for the moment matching model are slightly larger than the median predicted values for the biased model over the entire horizon. This is an interesting result: for identical process parameter values, the median of the moment matching models should be strictly smaller than the median of the biased model. This tells us that there is some mismatch between the parameters being estimated for the biased model and the moment matching model. This is something that we must be cautious about, especially since we are using parameters that have physical interpretations and modeling the evolution by embedding an ecosystem model that was designed to be deterministic.

We also found differences in performance between the two models when measured by IGN and CRPS. When measured by mean IGN across the 32 out of sample LAI measurements, the two models had similar performance, with the moment matching model performing slightly better (mean ${\text{IGN}}_{mm}=.3566$; mean ${\text{IGN}}_{biased}=.3850$). The moment matching model had marginally better IGN scores for the observations where neither model performed well (e.g. out of sample measurements two, three, and four, Figure 4, top panel), and the models had similar IGN score performance otherwise. When measured by mean CRPS across the 32 out of sample LAI measurements, the moment matching model showed a much better performance (mean ${\text{CRPS}}_{mm}=.2389$; mean ${\text{CRPS}}_{biased}=179.93$). Towards the end of the prediction horizon the CRPS for the biased model quickly increases, while the CRPS for the moment matching model stays comparatively small. We believe that this happens because of an accumulation of bias in the biased model combined with the heavy tails of the lognormal distribution.

Suppose that we are interested in finding a transformation for random variable $X$, $X\sim\text{Lognormal}(\mu^{*},\sigma^{*2})$ such that $\mathbb{E}[X]=\mu,\mathbb{V}[X]=\sigma^{2}$, with pdf given by: