Latent trajectory models for spatio-temporal dynamics in Alaskan ecosystems

The Alaskan landscape has undergone substantial changes in recent decades, most notably the expansion of shrubs and trees across the Arctic. We developed a dynamic statistical model to quantify the impact of climate change on the structural transformation of ecosystems using remotely sensed imagery. We used latent trajectory processes in a hierarchical framework to model dynamic state probabilities that evolve annually, from which we derived transition probabilities between ecotypes. Our latent trajectory model accommodates temporal irregularity in survey intervals and uses spatio-temporally heterogeneous climate drivers to infer rates of land cover transitions. We characterized multi-scale spatial correlation induced by plot and subplot arrangement in our study system. We also developed a Polya-Gamma sampling strategy to improve computation. Our model facilitates inference on the response of ecosystems to shifts in the climate and can be used to predict future land cover transitions under various climate scenarios.


Introduction
Climate change can impact ecosystems by altering vegetation composition.A prominent example is the expansion of shrubs, coined "shrubification," due to warming climates in northern Alaska (Swanson, 2013;Brodie et al., 2019).The encroachment of woody plants reduces erosion (Tape et al., 2010), increases fire frequency (Higuera et al., 2008), decreases albedo, and furthers warming (Chapin III et al., 2005).Remotely sensed imagery provides readily available high-resolution information about land cover and are commonly used to understand landscape changes (Svenningsen et al., 2015).In this study, we analysed pairs of historic and contemporary aerial images acquired across central Alaska.Our multi-scale multivariate spatio-temporal model provides inference about the rates of climate-driven land cover transitions among 5 major ecotypes.
Discrete-time Markov Chains (DTMCs) characterize a sequence of state changes by their transition probabilities and are commonly used to model landscape change (Baker, 1989), community dynamics (Hill et al., 2004), and plant succession (Logofet and Lesnaya, 2000).
In our study, the state space is finite and specified by ecotypes, and state changes occur in discrete time because of distinct growth seasons in Alaska.However, transition probabilities of a DTMC depend directly on the observed states.In the presence of sampling irregularity, modeling transition probabilities by a DTMC would require imputation of "missing" states at a temporal resolution typically defined by the smallest common sampling interval.Our objective is to learn about transition probabilities at the climate scale of 30 years, whereas the imagery pairs were collected 25 to 32 years apart.A naive DTMC would impute intermediate states annually to account for sampling discrepancy (suggesting that over 90% of the states need to be imputed).Without regularization on the transition mechanism or techniques such as multiple imputation (Scharf et al., 2017(Scharf et al., , 2019) ) there is no guarantee that the sequence of states inferred by the DTMC is representative of the progressive changes in plant communities.
The evolutionary mechanism of plant succession often sustain dependence beyond the temporal resolution of a DTMC.Higher-order Markov chains and semi-Markov models provide additional flexibility in temporal dynamics but dependence may remain at the state level (Moore, 1990;Lazrak et al., 2010).We characterize state changes using latent spatiotemporal processes in a logit-transformed probability space.We perceive such latent processes as an ecosystem analogue to animal movement (McClintock et al., 2014;Hooten et al., 2017) and refer to the spatio-temporally evolving state probabilities as the "latent trajectories" of the ecosystems.Ecosystem trajectories are commonly used in ecology and environmental science to describe the ecosystem dynamics over time (Locatelli et al., 2017;Lamothe et al., 2019), and in our method such trajectories are represented by their ordination in the latent probability space.Our latent trajectory representation redirects dependence between states to that between latent locations in the transformed probability space, thereby circumventing imputation of intermediate states.Our approach is related to spatial process models for non-Gaussian data via generalized linear modeling (Diggle et al., 1998;Finley et al., 2009).Jin et al. (2013) and Berrett and Calder (2016) respectively demonstrated the utility of spatial generalized linear (mixed) models in land cover classification.When temporal dynamics are involved, Bradley et al. (2019) showed that multinomial spatio-temporal mixed effects models can be used to analyze high-dimensional longitudinal data.Modeling latent processes in a continuous space provides additional smoothness to modeling observed processes in a discrete space and leads to consistent inference on ecotypes from year to year.
Different trajectory models reflect various evolution mechanisms of ecosystems we are able to accommodate.For example, the development of a primary forest may be represented by a random-walk-with-drift model directing at high probability regions of the dominant community.Successional impact of invasive species may be represented by interactions between the different dimensions of a latent trajectory, and extreme weather events that trigger secondary succession may be incorporated as Dirac delta functions in the latent trajectory.We illustrate these concepts in detail in Section 3.2.Further, our model achieves sampling efficiency through Pólya-Gamma data augmentation (Polson et al., 2013).Similar to the Albert-Chib data augmentation for probit regression (Albert and Chib, 1993), sampling auxiliary Pólya-Gamma random variables in a logistic regression model promotes conjugacy of linear predictors.We extend the Pólya-Gamma approach to multinomial logistic regression and incorporate multi-scale spatial correlation for the case study.

Imagery Data
The data motivating our analysis comprise 200 aerial imagery pairs acquired across the NPS Arctic Inventory and Monitoring Network (ARCN) from 1977 to 2010.A pair of aerial images consists of a georeferenced high-resolution color digital aerial image taken between 2008-2010 and a scanned and georeferenced color-infrared aerial image taken between 1977-1985 of the same plot on a systematic grid over ARCN (Figure 1).Each image was further divided into 37 contiguous regular hexagonal subplots during processing (Figure 2) and their dominant ecotypes were visually determined according to the scheme developed for the ARCN Ecological Land Survey and Land Cover Map (Jorgenson et al., 2009).We do not account for classification error in this study because a single highly-experienced observer processed all imagery pairs (Swanson, 2013).We categorized the 44 ecotypes developed for the Land Cover Map by their vegetation biomass and composition into the following 5 major ecotypes: The above more generalized ecotype categorization accounts for 90% of the total transitions between the original 44 mapped ecotypes.The 10% of transitions unaccounted for are intracategory.Table 1 summarizes the observed transition frequencies by subplots.Most of the subplots experienced no change over the study period.Among the 5 major ecotypes used for our analysis, the most frequent transitions were from Low Shrub to Forest, Low Shrub to Tall Shrub, and Other to Barren.Empirical studies suggested that most Shrub to Forest transitions occurred by post-fire succession, and most Low Shrub to Tall Shrub transitions occurred by tundra shrub increase; most transitions between Barren and Other occurred by fluvial processes in riverine environments, and a few other transitions occurred by thermokarst (Swanson, 2013).

Data Model
We let y i,s,t be a K-dimensional standard unit vector with one denoting the observed state (ecotype) of plot i, subplot s, at time t.For i = 1, . . ., n I (number of plots), s = 1, . . ., n S (number of subplots per plot), and t ∈ N (years), we model y i,s,t as where The mapping is known as a stick-breaking transformation similar to the stick-breaking process used in the construction of Dirichlet processes (Ishwaran and James, 2001).The kth element of pi,s,t represents the conditional probability that y i,s,t is in state k given it is not in any of the states 1, . . ., k − 1.This allows us to express the probability mass function of y i,s,t as a product of conditional binomials, where N i,s,t,1 = 1 and N i,s,t,k = 1 − r<k y i,s,t,r for k = 2, . . ., K − 1.The stick-breaking transformation allows us to exploit conjugacy in each of the binomial models through Pólya-Gamma data augmentation (Linderman et al., 2015), and we describe this strategy in Section 3.3.Although the stick-breaking process implies a prior stochastic ordering, it is not of practical concern when the model is data driven and the state space finite.In addition, the stick-breaking representation of multinomial logistic regression can be more economical than the alternative representation in Holmes and Held (2006) because it does not require evaluating a proportionality constant for every state-specific parameter.

Latent Trajectory Model
We specify a logit function that maps pi,s,t to η i,s,t for plot i, subplot s, at time t in the (K − 1)-dimensional real space as logit ( pi,s,t ) = η i,s,t .We model the dynamics in pi,s,t through a trajectory in the logit-transformed probability space (η-space).For illustration, Figure 3 shows a two-dimensional η-space where each location is translated to a threedimensional vector representing its probabilities in Shrub (blue), Forest (green), and Barren (yellow) states.The locations are colored by their most probable states.The simulated ecosystem trajectory starts at η 0 = (−0.4,0.8) (black point), which corresponds to state probabilities p 0 = (0.40, 0.41, 0.19) , indicating that the ecosystem is likely to be Shrub or Forest and unlikely to be Barren.At the end of the trajectory (red point), the ecosystem is at , which corresponds to state probabilities p T = (0.45, 0.07, 0.48) , indicating that the ecosystem is likely to be Shrub or Barren and unlikely to be Forest after time T .
An ecosystem in the η-space will always have non-zero probability in all states; however, state probabilities can become highly concentrated as the ecosystem departs from the origin (e.g., at η t = (7, 0) , the associated state probabilities, p t , are almost (1, 0, 0) ).
Among an array of discrete-time continuous-space trajectory models, we specify a randomwalk-with-drift model because it is a simple model that accommodates a temporal trend.The drift in the latent trajectory indicates the direction of ecological succession.For example, post-fire colonization of forests is represented by a drift vector directing from regions with high barren probabilities to regions with high forest probabilities, and the magnitude of the drift implies the rate of succession.We have, for k = 1, . . ., K − 1, where δ i,s,t,k is the spatio-temporally varying drift and e i,s,t,k represents uncertainty in the trajectory.The latent location of the ecosystem at time T is therefore We model the initial conditions, η i,s,0,k , with landscape covariates, h i,s , as follows where the Gaussian random effects, ζ i,s,k ∼ N 0, σ 2 ζ , provide additional flexibility.We model the drifts, δ i,s,t,k , using climate covariates, x i,s,t , as follows We further decompose movement uncertainty, e i,s,t,k , into two-levels of spatial random effects.
At the plot level, we account for correlation among locations on a systematic grid (Figure 1), ξ t,k , using a geostatistical model.At the subplot level, we account for correlation among the contiguous regular hexagons within the same plot (Figure 2), i,t,k , using an intrinsic conditional autoregressive (ICAR) model, so that for t = 1, . . ., T and k = 1, . . ., K − 1, where D represents geodesic distances (in kilometers) between plots, W is the adjacency matrix for subplots within a plot, and R is the diagonal matrix with row sums of W as diagonal elements (Ver Hoef et al., 2018).The ecosystem trajectory at plot i, subplot s, in the kth dimension can be summarized as The image pairs in our study have different beginning and ending years, and the dependence structure of the latent trajectories need to account for overlapping time intervals at different plots.The spatio-temporal covariance at the plot level is therefore The spatio-temporal covariance at the subplot level (within the same plot i) is where Q = (R − W ) −1 .To maintain propriety in the ICAR model, we constrain the subplot level random effects to sum to zero by sampling i,s,t,k via conditioning by kriging (Rue and Held, 2005).
We specify exchangeable Gaussian priors for the covariate coefficients, α k and β k , for k = 1, . . ., K − 1, and we specify Inverse-Gamma priors for the variance parameters, σ 2 ζ , σ 2 , and The range parameter φ is given a uniform prior bounded above by 1/3 of the maximum distance between plots.We provide a full description of the priors used for the case study in Web Appendix B.

Pólya-Gamma Data Augmentation
Every binomial component of the multinomial likelihood in (3) can be expressed as y i,s,t,k ∼ Binom(N i,s,t,k , logit −1 (η i,s,t,k )) for i = 1, . . ., n.By Theorem 1 in Polson et al. (2013), the following integral identity holds for the binomial likelihood [y i,s,t,k |η i,s,t,k ], where κ i,s,t,k = y i,s,t,k − N i,s,t,k /2 and ω ∼ PG(N i,s,t,k , 0).Because the right hand side of (11) contains a Gaussian kernel when conditioned on the Pólya-Gamma random variable ω, normal conjugacy holds for η i,s,t under a normal prior, η i,s,t ∼ N (µ i,s,t , Σ i,s,t ).The posterior where Conjugacy also holds for the Pólya-Gamma random variables.The posterior distribution of for k = 1, . . ., K − 1.The Pólya-Gamma approach aligns well with our stick-breaking representation of the multinomial likelihood and facilitates conjugacy in the linear predictors and the spatial random effects of the latent trajectory model.Although Johndrow et al. (2018) suggested that Metropolis-Hastings algorithms may be more efficient than data augmentation schemes including Pólya-Gamma with imbalanced categorical data, empirical examination of the traceplots in our study did not suggest significant autocorrelation.The data augmentation approach circumvents tuning, and is therefore particularly helpful to implement complex models such as ours.We provide a detailed description of the MCMC algorithm developed for this study in Web Appendix A.

Transition Matrix
Transition matrices can be used to compare our model inference to the observed transition frequencies and illustrate land cover changes under future climate scenarios.Our model allows us to infer instantaneous state probabilities from latent locations in the η-space, and we obtain elements of the transition matrix as a derived quantity.Because the transition mechanism of latent trajectory models differs entirely from that of DTMCs, the derived transition matrices do not possess Markovian properties such as phase-type distributions and stationary distributions.Nonetheless, our model answers relevant ecological questions such as species persistence and climax community using posterior predictive realizations of states.A posterior predictive realization of the transition matrix over time T given the qth posterior sample, η i,s,0 and ∆ (q) i,s,T (see (8) for definition), for q = 1, . . ., Q, is derived as follows, 1. Sample a posterior predictive realization of the initial state as for i = 1, . . ., n I and s = 1, . . ., n S ; 2. Sample a posterior predictive realization of the state at time T as i,s,0 + ∆ (q) i,s,T for i = 1, . . ., n I and s = 1, . . ., n S ; 3. Calculate a posterior predictive realization of the transition probability from state k 0 to k 1 over time T , for k 0 = 1, . . ., K and k 1 = 1, . . ., K, as , where M (q) is the qth posterior predictive transition matrix.
The posterior mean predictive transition matrix is evaluated as q) , and element-wise credible intervals can be constructed by applying quantile functions to the Q realizations of M .
We use simulation to illustrate that our model is able to recover the covariate coefficients and spatial parameters in Web Appendix C.

Case Study
The landscape and the climate variables only vary at the plot level in our application due to their large-scale spatial resolutions.Summer temperature (July in Northern Hemisphere) and soil moisture are the two major drivers of tundra vegetation biomass and composition (Murray and Miller, 1982;Elmendorf et al., 2012).For the initial condition model in (4), we used the covariate vector, h i,s , that includes an intercept, mean July temperature in the first year (Celcius), aspect degree (azimuth clockwise from north), slope degree, interaction between slope and aspect, and elevation (feet).The interaction term quantifies potential for Low Shrub (α 03 ) indicate a high initial probability in the state.We estimated a positive temperature coefficient (α 11 ) and a negative elevation coefficient (α 51 ) for Forest, suggesting that forest ecotypes are likely to occupy plots at low elevations with warm growing seasons.
The model fit also suggested that barren ecotypes are likely to occupy plots at high elevations (α 54 ) with steep slopes (α 34 ).
The β coefficients explain the temporal dynamics of, and the effect of climate change on, land cover state probabilities.The time coefficients for Forest (β 01 ) and Tall Shrub (β 02 ) have positive estimated 95% credible intervals, suggesting that probability mass will accumulate in these states over time.The coefficients for Low Shrub (β 03 ) and Barren (β 04 ) both have zero-overlapping estimated 95% credible intervals; however, as the probabilities in Forest and Tall Shrub grow, probabilities in other states will likely decline.We estimated a positive temperature coefficient for Forest (β 11 ) and a negative temperature coefficient for Barren (β 14 ), suggesting that warmer climates will lead to more frequent forest ecotypes and cooler climates will lead to more frequent barren ecotypes.We estimated negative precipitation coefficients for Tall Shrub (β 22 ) and Low Shrub (β 23 ), suggesting that drier climates are conducive to more frequent shrub ecotypes.We obtained posterior predictive realizations of transition probabilities following Section 2.5 for a variety of climate scenarios.Although inference on the covariate coefficients depends on the order of ecotypes, the derived transition matrices demonstrate consistent patterns regardless of ecotype ordering.Figure 5a shows the posterior mean predictive transition matrix for the case study, which is validated by the empirical transition frequencies in Table 1.
Temperature and precipitation are highly variable across Alaska.Nonetheless, the state has experienced an overall warming more than twice as fast as the contiguous U.S. in recent decades, with the most dramatic changes in spring and winter (Stewart et al., 2013).Studies project that the annual mean temperature will increase from 4 to 10 degrees Celsius by the end of this century under higher emission scenarios or from 2 to 6 degrees Celsius under lower emission scenarios (Stewart et al., 2017).The annual precipitation in Alaska is also projected to increase by 10% or more by mid-century (Stewart et al., 2017).We projected land cover transitions under the assumption of uniform climate change over the study area for the purpose of demonstration.Our model inference can be used in conjunction with spatially detailed climate forecasts to obtain more realistic predictions on landscape transformation.
Figure 6a shows the posterior mean predictive transition matrix under a high emission scenario where we assume that July temperature increases by 8 degrees Celsius and daily precipitation increases by 2mm (equivalent to a 730mm annual increment) uniformly in the study region from c. 1980 to 2100.As suggested by the β 0 estimates, probabilities will accumulate in Forest from all other states in c. 120 years.Most predicted transitions into Forest are from Other, Tall Shrub, and Low Shrub, possibly due to succession facilitated by warming climates.In comparison, there are fewer transitions from Barren to Forest, possibly due to landscape factors (e.g., temperature and elevation) that limit forest expansion.There are also few transition from Low Shrub to Tall Shrub, possibly because low shrubs undergoing succession have passed the state of Tall Shrub and reached the state of Forest after over a century.Figure 6b shows the posterior mean predictive transition matrix under a low emission scenario where we assume that July temperature increases by 4 degrees Celsius and daily precipitation increases by 2mm uniformly in the study region from c. 1980 to 2100.Table 4 (Web Appendix D) summarizes the posterior predictive mean transition probabilities and their associated uncertainty under the high and the low emission scenarios, respectively.The posterior mean transition probabilities are generally lower under the low emission scenario than those under the high emission scenario, although the differences between the two scenarios are not great.Our predicted transition patterns agree with the observed patterns in Figure 5a where the most frequent transitions are those from low biomass categories to high biomass categories.The transition probabilities are magnified under hypothetical warming scenarios and serve to support the association between Arctic greening and climate change on the global scale (Epstein et al., 2004;Harsch et al., 2009).Lastly, our predictions of landscape transitions do not account for uncertainty in climate forecasts and the variability in our predictions increases with the length of projection into the future.

Discussion
We presented a latent trajectory model for landscape change while accounting for spatiotemporal dependence.Our model characterizes dependence in the logit-transformed probability space and leverages computational efficiency through Pólya-Gamma data augmentation.We demonstrated dependence structures that manifest evolutionary mechanisms; therefore, our model is most suitable when the transition process takes place at a lower temporal frequency than that of the data collection.The latent trajectory models provide insight into the cumulative effect of long-term warming on Alaskan landscape and allow inference about state distributions and transition probabilities over flexible time intervals.
We decomposed the variance in the latent trajectory process into three sources: the uncertainty in the initial conditions, σ 2 ζ , the plot level uncertainty, σ 2 ξ , and the subplot level uncertainty, σ 2 .The three parameters are identifiable because σ 2 ζ is informed by both the historic and the contemporary images, whereas σ 2 ξ and σ 2 are informed by the contemporary images and modeled at different spatial scales.The estimated σ 2 ζ is larger than the other two variance parameters (Table 4) because the effects of σ 2 ξ and σ 2 are scaled by time.The large initial variance may be attributed to few observations in the temporal domain and scarcity of transitions in our case study.The uncertainty in the initial conditions could hinder prediction at unobserved locations and potentially confound with the latent movement processes.We regularize the estimation of initial conditions by incorporating spatial structure and leveraging the fact that the observations at different plots were staggered in time.Further, the uncertainty may be reduced with more temporal replicates at each plot and more diverse transition types.The α estimates are substantial in magnitude.They produce initial conditions with highly concentrated probabilities, which demonstrates our learning about the temporal dynamics of the transition process.An uneven initial condition (with a significant probability in one state and negligible probabilities in the other states) indicates that the ecosystem is unlikely to experience any systematic change in the immediate future, possibly due to limiting environmental factors explained by the covariate vector h i .On the other hand, the β estimates representing movement along the latent trajectories are smaller because fundamental changes in ecosystems usually take place over the course of centuries.As such, our model was able to account for small-scale temporal irregularity within the large-scale ecological process.
There are several ways to extend our latent trajectory model.A constant time coefficient (β 0 ) implies that, over time, probability mass will likely converge in one state as an ecosystem moves toward regions with high probabilities in that state.Such behaviors may be appropriate for some types of succession (e.g., primary/secondary successions where the dominant state is stable), but not otherwise (e.g., seasonal successions where the dominant state alternates frequently).A more flexible process model could use higher-order terms, interactions, or basis functions to represent the temporal trend.For succession mechanisms that are characterized by the coexistence or competition between species, we can explicitly model interactions between states through the covariance of the K − 1 dimensions of η.Lastly, environmental events such as fire and flooding often result in substantial changes in the ecosystem.These events can be important predictors of transformation and incorporating these events into dynamic models is a potential area of future research.In our application, however, only 2 out of the 200 plots had a record of fire prior to the collection of historic images, and a larger spatio-temporal domain is required to capture the effects of fires.
Our hierarchical framework aids in integrating various data sources, some of which may arise from other studies on Alaskan vegetation change.For example, Scharf et al. (2022) and Raiho et al. (2022) developed a model to identify landscape factors that resist climatedriven vegetation change.Although limited by the absence of temporal replicates, Raiho et al. (2022) made inference over a greater spatial domain and performed statistical model selection to distinguish key variables.Their analysis output -a measurement of ecosystem robustness to climate change given its landscape covariates -could be used in our latent trajectory model because "robust" ecosystems should move less in the η-space due to climate change.On the other hand, our model is based on temporal replicates, and we could compare our inference with the results from Raiho et al. (2022) in a future study to visualize how our estimated climate effects relate to the robustness scores.

Figure 1 :
Figure 1: A pair of images collected in 1979 and 2008, respectively, in the Noatak National Preserve, Alaska, as represented by the red dot on the inset map showing all plot locations for the case study.The historic image (AHAP color-infrared photo, left) consists mostly of herbaceous and low shrub vegetation.The contemporary image (small-format true color photo, right) consists mostly of tall shrubs.Source: Swanson, 2013.

Figure 3 :-Figure 4 :
Figure 3: A simulated trajectory in a two-dimensional logit-transformed probability space.Each location is colored by the most probable state.A simulated trajectory shows changes in state probabilities as an ecosystem travels across the latent space.

Table 1 :
Summary of transition frequencies by subplots.Rows indicate categorization in c. 1980 and columns indicate categorization in c. 2010.