Developing hierarchical density-structured models to study the national-scale dynamics of an arable weed

. Population dynamics can be highly variable in the face of environmental hetero- geneity, and understanding this variation is central in the study of ecology. Robust management decisions require that we understand how populations respond to management at a range of scales, and under a broad suite of conditions. Population models are potentially valuable tools in addressing this challenge. However, without adequate data, models can fail to pro- duce useful results. Populations of arable weeds are particularly problematic in this respect, as they are widespread and their dynamics are extremely variable. Owing to the inherent cost of collecting data, most studies of plant population dynamics are derived from localized experi-ments under a small range of environmental conditions, limiting the extent to which variance in population dynamics can be measured. Density-structured models provide a route to rapid, large-scale analysis of population dynamics, and can expand the scale of ecological models that are directly tied to data. Here we extend previous density-structured models to include environmental heterogeneity, variation in management, and to account for inter-population variation. We develop, parameterize, and test hierarchical density-structured models for a common agricultural weed, black-grass ( Alopecurus myosuroides ). We model the dynamics of this species in response to crop management, using survey data gathered over 4 yr from 364 fields across a network of 45 UK farms. We show that hierarchical density-structured models provide a sub- stantial improvement over their nonhierarchical counterparts. Using these models, we demonstrate that several alternative crop rotations are effective in reducing weed densities. Rotations with high wheat prevalence exhibit the most severe infestations, and diverse rotations generally have lower weed densities. However, a key outcome is that in many cases the effect of crop rota- tion is small compared to the high variability arising from spatiotemporal heterogeneity. This result highlights the need to monitor and model population dynamics across large spatial and temporal scales in order to account for variation in the drivers of plant dynamics. Our framework for data collection and modeling provides a means to achieve this.


INTRODUCTION
Populations of many species are distributed across large spatial scales, and are subject to highly variable environments. As a consequence, they can exhibit dynamics that are variable in both time and space, and depend on local conditions and context (Levin 1992, Dunning et al. 1995, Lundberg et al. 2000, Coutts et al. 2016. Managing populations in variable environments requires detailed knowledge of environment-driven spatiotemporal dynamics, whether the focus is on balancing natural resources and conservation (Tscharntke et al. 2005, Flesch and Steidl 2010, Damschen et al. 2019), or eradicating problematic species (Freckleton et al. 2000, Bianchi et al. 2006, Ziska et al. 2011. Understanding how populations respond to environmental drivers is essential for effective management, especially considering rapid rates of global change (Sutherland 2006, IPCC 2013, Sutherland et al. 2013). However, gathering data of sufficient quality to encapsulate both the full range of population responses, and the associated environmental drivers, remains extremely challenging. Gathering data over large scales is expensive and time consuming, leading to a trade-off between data extent and quality. As a result, ecological studies are typically reliant on intensive small-scale studies, and often only capture demographic variation at a few locations under relatively static conditions (Clutton-Brock et al. 1996, Coulson et al. 2001, Bremer and Jongejans 2010, Dalgleish et al. 2010, Garnier et al. 2018. These studies may produce high-quality data, but often fail to encompass the full range of conditions that populations experience (Forman 1995, Miller et al. 2004, and there is typically a lack of large-scale and detailed demographic data representing more than a few locations (Salguero-Gomez et al. 2015, Gurevitch et al. 2016. Often, data from one or a small number of well-studied populations are used to generalize demography (Silvertown et al. 1996 with the assumption that this accurately depicts dynamics of populations over relevant scales (Burns et al. 2010, Crone et al. 2011. However, generalization of demography can be extremely problematic, as reliable extrapolation of demographic metrics is limited to local scales (Coutts et al. 2016), and demographic parameters estimated for one population may be inappropriate for others (Che-Castaldo et al. 2018). Generalizing demographic traits is unreliable, partly due uneven sampling across space and phylogeny. Sampling bias exists across taxa, which confounds analysis, as closely related species will often occupy ranges with similar environmental conditions (Coutts et al. 2016). Similarly, local population models are also limited by the availability of suitable data, as temporal heterogeneity can make it difficult to obtain sufficient data to accurately estimate the environmental variance in key parameters (Cousens 1995, Freckleton et al. 2006. Without adequate parameterization, demographic models can fail in the face of uncertainty (Freckleton et al. 2008) and result in poor forecast accuracy (Crone et al. 2013).
To effectively model population dynamics over large spatial extents, sampling must extend over multiple populations in order to accurately capture the variance and covariance in demographic parameters due to the environment (Crone et al. 2011, Coutts et al. 2016, Che-Castaldo et al. 2018, Quintana-Ascencio et al. 2018, Damschen et al. 2019. Density-structured modeling is a method that addresses this challenge by facilitating rapid data collection across multiple populations, while also permitting accurate characterization of local population dynamics (Taylor and Hastings 2004, 2017, Mieszkowska et al. 2013, Tredennick et al. 2017). Instead of a continuous measure of abundance, density-structured methods discretize local population numbers into ordinal density 'states' and model the probabilities of transition between these categories. This method addresses the problems faced by conventional methods in two ways. First, it enables fast data collection over large spatial and temporal scales, because time-consuming counts of individual plants are replaced by rapidly estimated density states. Second, model parameterization is simplified because dynamics are summarized by transition matrices. Without the need for in-depth quantification of key demographic parameters, transition matrix models are more robust to numerical instability (which may be caused by parameters that are sensitive to environmental heterogeneity) and parameterization error, both of which are potentially pathological for demographic models (Freckleton et al. 2008, 2011, Freckleton and Stephens 2009. Moreover, as these models are inherently empirical they facilitate collection of large amounts of data, which itself reduces the risk of estimation error. These features mean that density-structured models can be used to quantify a much wider range of responses to different stimuli, which are directly underpinned by empirical data. Densitystructured models enable studies to both encompass numerous locations over environmental gradients and accurately capture a large range of population responses to their environment . Although density-structured population models are a promising technique, their application has been limited to the short term (Freckleton et al. 2017) and often without accounting for location-specific effects Hastings 2004, Mieszkowska et al. 2013). In order for density-structured models to inform the management of widespread populations they need to account for variability across multiple scales. One of the limitations of density-structured models is that if there are S sites and K density states, and transition probabilities differ between sites, then S(K 2 − K) independent parameters need to be estimated, in order to measure the variation in population dynamics across all sites. This is because there are K 2 possible transitions between K density states in a single site, with sum-to-one constraints meaning that K 2 -K parameters are required to characterize dynamics in one population. Furthermore, it may be difficult to estimate transition probabilities if not all density classes are observed at every site.
A powerful solution is to model the hierarchical organization within such data, accounting for effects that reflect spatial or temporal dependence between parameters. We define hierarchical models (also known as multi-level models) as models with multiple variance components to allow partitioning of variance components belonging to different groups (Gelman and Hill 2007). Incorporating hierarchical effects into densitystructured models could improve the estimation of parameters of density-structured models in a number of ways. First, hierarchical models can capture the variation in density-structured dynamics across multiple populations, and simultaneously model the dynamics of individual populations and landscapes (e.g., Freckleton et al. 2017). Second, hierarchical models can deal with the problem of estimating large numbers of parameters, for which there is little information, through partialpooling (Gelman and Hill 2007), allowing all the data to inform the estimates for each population. In the context of density-structured models, this benefit is potentially an important advance given the likely need to estimate large numbers of parameters.
In this paper, we develop hierarchical density-structured models for populations of an arable weed subject to a range of different management regimes distributed across a landscape. The threat to arable farming posed by weeds is escalating due to evolved resistance to multiple herbicides, that makes them increasingly difficult to manage (Moss and Clarke 1994, Powles and Yu 2010, Moss et al. 2011, Jasieniuk et al. 2018, as well as exacerbated risks of invasion due to climate change (Dukes andMooney 2008, Peters et al. 2014).
It is particularly difficult to model arable weeds for several reasons. First they often occur in complex and highly fragmented landscapes, existing as small populations in fields that are subject to variable management practices (Wiese et al. 1997, Mack et al. 2000, Tilman et al. 2001, Walker et al. 2005. Second, weeds are subject to very high levels of control, with populations often close to an extinction threshold, making population dynamics potentially numerically unstable (Freckleton et al. 2008). Finally, many demographic parameters are extremely difficult to estimate from field data, due to population-specific variation in environmental variables such as soil, climate, and crop varieties, which can lead to extreme variability in dynamics (Cardina et al. 1997, Wallinga et al. 1999, Freckleton and Watkinson 2002a, b, Freckleton and Stephens 2009, Lima et al. 2012, Lutman et al. 2013. The consequence is that weed population dynamics can be very challenging to predict because data are typically limited to single populations (Gonzalez-Andujar and Fernandez-Quintanilla 1991, Buhler 1999, Freckleton et al. 2000, Buckley et al. 2003, Colbach et al. 2005, 2006, Metcalfe et al. 2018). However, density-structured methods have proved successful in monitoring and modeling these populations using large-scale survey data (Taylor and Hastings 2004, 2017.
We develop hierarchical density-structured models to estimate the effect of crop rotation, an integral part of arable farming and weed management, on the population dynamics of a common arable weed (Zacharias and Grube 1984, Liebman and Dyck 1993, Melander et al. 2005. Specifically, we have two objectives. First, we aim to find the best candidate hierarchical density-structured model for accurate characterization of population dynamics. Second, we apply these hierarchical densitystructured models to a large-scale data set to examine the impact of crop rotation on weed dynamics. We then assess the impact of alternative control strategies on one of Europe's most economically damaging weeds, blackgrass (Alopecurus myosuroides) (Moss and Clarke 1994, 2011, Lutman et al. 2013. To our knowledge, this study represents one of the largest studies of weed population dynamics to date.
Our models demonstrate extreme levels of between-field variability in weed density, relative to the effect of rotation, highlighting that quantification of spatiotemporal dynamics at such scales is vital in assessing the effectiveness of management.

Study system and survey
Densities of black-grass (Alopecurus myosuroides Huds., Poaceae) were recorded in a series of repeated surveys from 2007 to 2010. This data set includes 682 repeated field-scale surveys from 364 individual fields nested within 45 arable farms throughout the counties of Norfolk, Lincolnshire, and Bedfordshire in the UK. The density-structured survey method, detailed in Queenborough et al. (2011), involved repeated surveys of individual fields to map weed densities in a given survey year. Each field was divided up into a set (median = 438, IQR = 245) of 20 × 20 m survey quadrats, predefined using a GPS system. Researchers walked the fields recording the densities at each quadrat as one of five categories: absent (A), low (L), medium (M), high (H), or very high (VH). These categories are roughly delineated by the quartiles of a lognormal density distribution, and were chosen based on previous studies that have critically evaluated the method to demonstrate high within-and between-observer repeatability .

Modeling density-structured data
A density-structured model has the form where n is a vector of length K (where K is the number of density states) whose elements are the probabilities of each density state at time t, and T is a K × K columnstochastic matrix of transition probabilities The transition matrix, T, defines the population dynamics. Diagonal entries of T represent probabilities (p) that a quadrat in a given state will remain in that state for the next survey, and off-diagonals represent the transition between states between years. For example, p 11 is the probability that a quadrat in state 1, will remain in state 1, and p 12 is the probability that a quadrat in state 2 will transition to state 1. Eq. (1) defines a first-order Markov chain model that can be used to predict future density state distributions. Detailed explanations and evaluations of density-structured models can be found in Freckleton et al (2011).

Analyzing black-grass dynamics in response to crop rotation
We investigated the effect of management on blackgrass density by formulating density-structured models that simulate the impact of different rotations on weed density. In the context of weed populations, crop-rotation involves cyclic environmental disturbance, under such circumstances, density-structured dynamics can be investigated using periodic models (Skellam 1967, Lesnoff 1999, Cushing and Henson 2018, Mertens et al. 2002. Population dynamics under different rotations can be modeled by changing the transition matrix in successive time steps in the Markov model (Appendix S1: Eqs. S9, S10), and complex rotations can be modeled by changing the order of rotation-specific transition matrices. The only condition is that the destination crop of the previous matrix is the initial crop of the next, i.e., if one matrix models the transitions from wheat to barley, the next matrix in the sequence must model the transitions from barley to the next crop.

Model fitting
To analyze the impact of environmental variability on population dynamics through rotations, we first constructed a set of models that accounted for field-level spatiotemporal effects on transition probabilities. We parameterized transition matrices for each field-year (i.e., a field observed in a given year) observed in subsets of data representing three rotations, wheat-to-wheat, wheatto-oil seed rape (OSR), and OSR-to-wheat.
To estimate transition probabilities, we fitted latent variable ordered category logit models to our density-state data. These empirical models allowed easy conceptualization of drivers of weed densities and flexible parameterization. As such they can account for the variable dynamics present in weed populations (Freckleton et al. 2017, Gonzalez-Andujar and Hughes 2000, Freckleton and Watkinson 2002a, b) as they do not restrict transitions between non-adjacent categories (Agresti 2012:303).
In an ordered categorical model, the probability of observing a given category, k, at quadrat i is expressed in terms of a real-valued latent variable that reflects the true (unobserved) value. A linear predictor, η i , is defined for each quadrat, which is constructed from the row-vector of J quadrat-specific explanatory variables, x i , and the unknown column-vector parameter β i . β ij is therefore the effect of explanatory variable x ij on n i , such that The constraint β i1 = 0 is enforced to allow identifiability, a common practice in such models (Agresti 2012:297). The ordering of categories in this model is then enforced through a set of K − 1 (where K is the total number of categories), "cut-point" parameters, c i , where c 1 < c 2 < . . .c K−2 < c K−1 (Appendix S1: Eq. S6).
Although η is unobserved, we categorize outcomes according to the following rules, where Θ ik gives the probability of observing state k at quadrat i As there is potentially uncertainty when it comes to estimating parameters in these models with collinear observational data, we employed a Bayesian framework using the probabilistic programming language Stan (Stan Core Development Team 2017) to allow flexibility in parameterization and to account for this uncertainty. Full hierarchical and prior specifications are detailed in Appendix S1.

Alternative models
We constructed a series of five ordered category logistic regressions to compare alternative formulations for estimating transition probabilities in a hierarchical density-structured framework.

Model I: Global, nonhierarchical model
Model I is the formulation presented in Eqs. 3 and 4. This model formed the baseline for all subsequent models and incorporated the effect of source state (i.e., density state of quadrat i at time t) as covariates x i1 . . .x i5 , which take the form of indicator variables. Eqs. 3 and 4 describe the model for the probability of observing state k, conditional on source state j Therefore, p K1 in Eq. 2 is equivalent to Θ k where j = 1 in Eq. 3. Model I was included in the analyses as a baseline reference for Fieldlevel comparison and had no hierarchical variance components. Models II-IV (summarized in Table 1) took different approaches to modeling the hierarchy present in our data set.

Model II: Field-level intercept
The simplest model that accounts for between-field variability included a "field" effect in the construction of the linear predictor. The scalar intercept term γ f represented the field-level effect on the linear predictor within field f. Here the cut-point parameters c k . . . c K−1 remained as in Eq. 4

Model III: Field-level source state effects
The logical extension of Model II was to permit more flexibility in the construction of the linear predictor by allowing the source-state effect to vary between fields. The global effect of source state β ij and global cutpoints, c 1 . . .c K−1 , remained as in Model II. The source state effects governing the probability of transition between density states, β ij , are likely to be determined by the same process across a landscape. However, environmental heterogeneity mean they may vary at the field level. To account for variance in field-level source state effects we estimated field-specific slopes for the effect of source density state. This can be expressed via an additional the column-vector parameter, γ if , which represents the effect of source state j in field f on η. The addition of γ if allowed the effect of source state to vary between fields and aimed to account for the various drivers that affect changes in black-grass densities between surveys

Model IV: Field-level cut-points
Many applications of hierarchical modeling use the approach we have outlined above, accounting for grouplevel variation by including a term for group-level effects in the construction of the linear predictor. In an ordered category logistic regression, an alternative approach is to allow cut-point parameters to vary between each group (in our case field-years). The advantage of this method is greater flexibility in the estimation of transition probabilities, because the cut-points themselves, which control the conditional probability of an observation being in state 1, K, are able to vary between field. We implemented this method by estimating a set of cut-points for each field: Model V: Field-level cut-points and field-level intercept The final model is a combination of Models III and IV, with both random cut-points and a random intercept included in the linear predictor. As such the linear predictor was the same as in Model II, and the cut-point parameters were the same as in Model IV. For the purposes of simplicity in model selection, we use data from the three most common rotations, wheat-to-wheat, wheat-to-OSR, and OSR-to-wheat, to compare our models using a variety of cropping systems with different dynamics.

Assessing predictive performance and posterior checks
To assess model performance across the three rotational subsets selected for model fitting, we used leave-one-out cross validation (LOO-CV) implemented in R package loo version 1.1.0 (Vehtari et al. 2017), using LOO-IC (LOO information criterion) as a measure of relative predictive error. We also visualized model performance via graphical posterior predictive checks to assess any systematic prediction errors due to differences in parameterization between models. We simulated field scale density distributions from posterior probabilities and compared them to the observed distributions in each corresponding field. We compared the full distribution of density states as well as mean density states calculated for each individual field. Although the density-states are categorical variables they delineate the true underlying continuous distribution of black-grass, the mean density-state is therefore approximately proportional to the geometric mean of the true abundance and a useful measure for comparison.

Ecological analyses
To explore the impact of spatiotemporal variability on weed dynamics across different crop rotations, we used our best performing model (Model IV) to fit ordered category logistic regressions to each of observed rotations in Fig. 1. In each of these models, we estimate field-year matrices for the transitions observed in each field in a given year, so that each matrix incorporates both spatial and temporal variability in transition probabilities. For example, we parameterized seven matrices from observed transitions in fields rotated from barley into wheat, and 38 matrices of wheat into barley. As the data are fragmented and often lack uninterrupted observations within a single field, we use permutations of field-level matrices within rotations to allow us to simulate and compare more complicated rotations than we observe in the data. To examine black-grass dynamics under various crop rotations and environmental contexts we conducted three analyses using density-structured models.

Asymptotic and transient dynamics
Given enough time, density-structured transition models (with primitive matrices) will always converge on a stable density structure for any given point in the rotation cycle, because the sum-to-one constraints ensure that the dominant eigenvalue of any transition matrix is always one (Caswell 2001). A general approach to studying the dynamics of these models is to compare the stable density structures and the rates of convergence of different rotations. To investigate the dynamics of particular cropping alternatives we compared net transition matrices for a series of two-step rotations. The asymptotic dynamics of a two-step rotation, analogous to running the Markov chain for a large number of iterations, can be examined via net transition matrices, defined as the product of two field-level matrices in a given rotation where subscripts A and B on T represent the sequence of crops modeled in each transition matrix, and superscript denotes the field. T ij ABA is therefore the matrixmodeling transition probabilities of weeds between density states during the rotation of crops A to B to A, with the superscript denoting the specific permutation of field-years i and j. Permutation of field-year matrices as in Eq. 8 is required as many crop sequences were unobserved within a single field. This allows the simulation and assessment of the impact of site-specific effects on weed densities across a larger subset of rotations (Mertens et al. 2002(Mertens et al. , 2006. We generated every possible model of the form T ij WBW , where subscript W denotes a wheat crop, and B, denotes an intermediate crop (Table 2). We also parameterize models for continuous wheat (T ij WWW ) and continuous barley (T ij BBB ). For each of these rotations, we generate net matrices from each possible permutation of field-level matrices; thus, for the wheat to barley to wheat example, there are 38 field-level wheat to barley matrices (T i WB ) and seven barley to wheat matrices (T j BW ), making a total of 266 (38 × 7) net matrices (T ij WBW ). For each net matrix, we calculated two summary statistics that allowed us to examine the dynamics of different cropping systems.
FIG. 1. A transition matrix illustrating the rotation space covered by models fit to observed data. Models were fit to each nonzero entry of this transition matrix. The top number in each cell is the number of fields, and the bottom number is the overall density-state observations in that rotation. The color represents the number of fields in each observed rotation. OSR, oil-seed rape.
The stable density structure s, gives the proportion of sites in a field in each density-state at the asymptote of each rotation and is a measure of the distribution of overall population density. For each net matrix, s can be calculated from the leading eigenvector of a net transition matrix where v 1 is the right eigenvector of the matrix, and i indicates the element. s is a K length vector representing the field-level density state distribution in terms of proportion of the field in each category. To summarize asymptotic dynamics, we calculate the mean density state from the stable density structure for each net matrix, i.e., the proportion of a population occupied by each state multiplied by the integer value (1-5) of each state category. The rate of convergence to this structure is a measure of how quickly the population would reach equilibrium after a disturbance, and is governed by the relationship between the dominant and subdominant eigenvalues (noting that the dominant eigenvalue is 1) where λ 1 is the dominant eigenvalue and λ 2 is the second largest eigenvalue. The parameter ρ is the damping ratio and gives a measure of sensitivity in the face of perturbation, or the rate at which a population will approach its stable density structure. The higher this ratio, the faster the convergence (Caswell 2001:95).

Short-term projections using two-step rotations
From an agronomic perspective, prediction of shortterm dynamics is important because weed management objectives (e.g., leading to outcomes in terms of yield and profit) are typically measured on very short time scales. We therefore analysed population dynamics on timescales of 2 yr (e.g., following Freckleton et al. 2017). As with the asymptotic and transient dynamic analyses we constructed models for all possible rotations starting and ending with wheat, as well as continuous rotations of barley and wheat. Simulating wheat to wheat transitions ensures that transitions are meaningful, first as the severity of an infestation can be judged against a common crop, but also because wheat is the most economically valuable crop in the UK (Department for Environment, Food and Rural Affairs [DEFRA] 2018). We made two-step projections for each possible combination of field-level matrices, for three initial starting densities (i.e., the density-state distribution of the field at the beginning of the simulation), and represented typi- Thus, for the wheat-barleywheat example, there are 266 matrices for a single initial density and 798 (266 × 3) outcomes in total. We compared the average outcome and distribution of each rotation from each initial starting density, as well as the relative change compared to wheat. Notes: Includes the length of the rotation in years, the number of years a field spent in a break crop and the wheat dominance (proportion spent in wheat) of a rotation. Rotation denotes the sequence of wheat crops (W) and break crops (B), where B1 and B2 represent the first and second break crops in systems with two different break crops. To deconstruct the effects of environment and management, we used a life-table response-experiment (LTRE; Caswell 1989Caswell , 2001 to analyze the variation in the change in weed density between years due to local spatiotemporal effects (or field identity), initial densities, and rotation. This method uses the change in blackgrass density from the two-step projections described above as the response variable in a linear mixed-effects model (e.g., Freckleton et al. 2017), to account for variation in population structure and intrinsic dynamics In Eq. 11, Y i is the change in mean density state between years for an individual simulation outcome i.
Parameters μ ð1Þ f , μ ð2Þ j ,μ ð3Þ k , represent the effects of matrixpair f (i.e., the permutation of field-level matrices used in the projection), initial density j, and rotation k, respectively. Variances associated with these parameters are represented by θ 1,2,3 , and ε and σ are the residual error and its associated variance. The parameter α 0 represents the global intercept and was set to 1 in this analysis. To evaluate the variance components of field-level matrix pair, initial density, and rotation we estimate parameters μ ð1Þ f ,μ ð2Þ j , μ ð3Þ k as random effects. Although there were only three levels in initial density, it was estimated as a random effect for the sake of efficiency and still provides the estimate of the variance component as would be obtained from treating it as a fixed-effect.

Long-term dynamics: stochastic projections
To investigate the impact of rotational diversity and type of break crop on weed density in the context of spatiotemporal variability we implement a series of stochastic models where T i AB denotes the field-level matrix that models transitions between crop A and B in field i. Stochasticity is implemented through an iid lottery, so that the "fieldlevel" matrix for a particular step in a rotation is sampled randomly as an independent and identically distributed variable. For example, for rotation AB, T i AB is selected randomly from all fields in rotation AB, and for rotation BA, T j AB , is selected randomly from all fields in rotation BA. Modeling dynamics in this way allows us to construct complicated rotational strategies that incorporate the large scale spatiotemporal variation captured by our surveys. We project the density state forward 10,000 time steps (to ensure convergence) for every model described in Table 2, which includes rotations with multiple break crops. All models were started at a "low" density distribution, and mean density states were calculated at each time step. Across time series, we then compared means, variances, and coefficients of autocorrelation, for a series of rotations with different break crops, and rotational diversity (measured as the proportion of the dominant crop, wheat, in a rotation). From this, we examined overall rotational effects on blackgrass density, variability, and whether observed densities are likely to persist year to year.

Model fitting
Hierarchical models (Models II-V) that incorporate field-level effects into the transition probability have better predictive accuracy than our nonhierarchical model (Fig. 2a). All hierarchical models provide similar levels of predictive performance. LOO-CV provided the most support for the models that incorporate hierarchical effects through cut-point parameters. Models that only incorporated hierarchical structure in the linear predictor fared slightly worse than field-level cutpoint models, while there was no difference between in the field-level cut-point model (IV) and the combined linear predictor/cut-point model (V). Rotational subsets of our data all displayed the same order of model preference, suggesting that each model performs similarly under different cropping systems (Appendix S1: Fig. S1).
Through comparing observed and predicted fieldscale mean density states it is again apparent that hierarchical models provided a better fit to the data than our baseline model (Fig. 2b). Although slight improvements in terms of root-mean-square error (RMSE) were seen from Models IV and V compared to II and III, all hierarchical models have similar performance as in all cases the prediction error and 80% equal-tailed credible intervals were close to zero.
There were, however, notable differences in the predictive accuracy of the nonhierarchical models between rotational subsets. Variance in field-scale predictive error was much higher in fields rotating from wheat or OSR into wheat (Fig. 2b right and left panels), than from wheat into OSR (Fig. 2b middle panel). The lower RMSE resulting from the predictions of Model I was also accompanied by a slight tendency to overestimate field-scale density. This was not the case with the hierarchical implementations, where all models displayed smaller error and higher correlations between observed and predicted densities.
Hierarchical models had lower error across densitystate distributions for the wheat to wheat rotation (Fig. 2 c), and similar trends were evident in the other rotational subsets (Appendix S1: Fig. S1). Models II and III tended to underpredict the frequencies of absent states within fields, while having higher error around the frequencies errors for mean density state by rotational subset. Hollow dots are median errors and vertical bars are 80% equal-tailed credible intervals. Numbers next to each median are the root-meansquare-error (RMSE) of the difference between predicted and observed densities. Individual colored points represent difference between observed and predicted values for an individual field. Errors are relative to each rotational subset. A total of 15 points were excluded from the extremes of Model I in this plot to allow clearer visualization but were included in calculations of RMSE. (c) Field-level error between observed and predicted frequencies for each density state category in the wheat-to-wheat subset. Black dots are median field-scale differences between predicted and observed frequencies of each density state across all fields, vertical bars are 80% equal-tailed credible intervals. Colors represent each of our five models, and individual colored points represent the difference in observed vs. density-state frequencies for individual fields. A total of 43 points were removed from the extremes of Model I in this plot to allow clearer visualization. Density states are absent (A), low (L), medium (M), high (H), or very high (VH). of low and medium states. Models IV and V, on the other hand, had error distributions in all states close to zero and exhibited little evidence of systematic prediction error in the estimation of density-state distributions ( Fig. 2b and c).
The composite rotational matrices demonstrate that fields of continuous wheat exhibited higher probability of transition into higher density states than continuous barley (Fig. 3, top row). An absent quadrat in continuous wheat had a 0.28 probability of being occupied by low densities of black-grass the next year, but when rotating from barley to barley the same transition was 0.11. A similar tendency was apparent in matrices that modeled rotations into wheat (Fig. 3, bottom row). For example, there was a 0.18 probability that a patch of very high density black grass would remain in that state when rotating from peas to wheat. Conversely, rotations into an intermediate crop (Fig. 3, middle row) showed the bulk of the transition probability into lower density states.
Analysis of transient and asymptotic dynamics revealed that rotated systems stabilized faster but in general had slightly lower equilibrium densities (Fig. 4a, b). Equilibrium density states for continuous wheat were generally higher and more variable than in rotations that included break crops or continuous barley (Fig. 4a). Using barley or peas as a break crop resulted in the highest equilibrium density for rotated systems, while potatoes have the lowest (Fig. 4a). Continuous barley produced lower equilibrium densities than all rotations using break crops.
Populations of black-grass that had been subjected to a rotation demonstrated higher damping ratios than cropping of continuous barley or wheat, (Fig. 4b) and therefore faster convergence on a stable state distribution. Faster convergence was generally accompanied by lower equilibrium densities, with the exception of barley and pea rotations where the pairing suggested faster convergence on higher equilibrium densities. Although generally high damping ratios were paired with low equilibrium densities, there was considerable variability in these relationships (Appendix S1: Figs. S2, S3).

Short-term dynamics: two-step rotation
There were clear but weak patterns in the effect of rotation on mean density state in short-term simulations. After two-step rotations, field-scale mean density states differed between, but were highly variable within, each rotation (Fig. 5a). Continuous wheat and, to a lesser extent, continuous barley, showed sensitivity to the initial density, indicating that the effectiveness of rotation as a management tool depended on the starting conditions. However, some rotations seemed to be invariant to FIG. 3. Average transition probability matrices for given crop rotations. Matrices display the probability of transitioning from a state in year 1 (x-axis) to any other state in year 2 (y-axis). The darker the color the higher the probability of transition. Numbers in each cell are the estimated probabilities. The first row is the matrices for continuous barley and wheat, the second row are matrices for rotations out of wheat, and the final row are rotations back into wheat. initial state, for example rotation into OSR produced the same outcome regardless of initial conditions. For most rotations low initial densities resulted in higher final densities (Fig. 5b), populations with medium initial densities tended not to change, and populations at high initial densities saw drastic reductions. Relative to continuous wheat, populations at low initial densities in rotations of beans and sugar beet offer low reductions, while OSR and barley offered virtually no reduction. Conversely rotating to peas increased weed density. At higher starting densities all rotations, except into peas, offer considerable reduction in relative density state.
Within two-step rotations it was clear that spatiotemporal effects and initial density (i.e., the field-specific conditions) explained the most variation in final density state, while the actual management intervention (rotation) was marginal in its contribution (Fig. 6). Large- FIG. 4. Results from analysis of transient dynamics using damping ratios and mean density states from stable state distributions. Distribution of (a) field-level log damping ratios and (b) mean density states for each rotation, and continuous crops of barley and wheat (C. barley, C. wheat). Hollow points are the median value across all permutations of net matrices, red and black bars are 90% and 50% equal-tailed credible intervals, respectively.

August 2021
MODELING NATIONAL-SCALE PLANT DYNAMICS Article e01449; page 11 scale spatial and temporal heterogeneity (i.e., inter-field scale variation) had a large effect on the outcome of management. We estimated the variance in change in density due to the identity of the pair of matrices used in the model (spatiotemporal effects) as 0.45, variance of the initial density conditions as 0.44, and the variance of each of the rotations as 0.03, with a residual variance of 0.08.

Long-term dynamics: stochastic projections
Reducing the proportion of wheat in a rotation decreased overall black-grass density, but variability between break crops is evident (Fig. 7), and there is considerable variability in the trajectories of weed density of rotations when increasing proportion of wheat. For example, potatoes increase steeply in mean density from wheat proportions of 0.6 to 0.75, while the same increases in barley rotations were shallower. There was noticeable variation in the effect of break-crop on outcome density (Fig. 8a). Using potatoes as a break crop produced significantly lower overall black-grass density, while OSR and peas produced comparably high levels. Potatoes also produced significantly lower between-year variability (Appendix S1: Fig. S4a), while OSR, peas and barley produced the highest levels of between-year variability in black-grass densities. Similar trends were evident in the relationship between break crop and autocorrelation (Fig. S4b).

DISCUSSION
Density-structured models are a promising method to assess the impact of environmental drivers on population dynamics over large spatial scales. Hierarchical implementations of these models allow simultaneous modeling of local and landscape-scale dynamics. Hierarchical density-structured models provided considerable improvement compared to nonhierarchical model performance in terms of point-wise predictive capacity. The key finding is that the most flexible models (Models IV and V) performed the best; they had the lowest error in terms of out-of-sample prediction, and next to no systematic error in predicting field-scale density-state distributions. This result reflects the fact there is a large amount of among field variation in the data. These models incorporated field-level effects via hierarchical cutpoints, which allows more flexibility and better predictive performance than models that account for field specific effects in the linear predictor. These more flexible models may also account for more variation in factors that act at the field scale such as observation error (e.g., the effect of different developmental stages of the crop or observation conditions; Queenborough et al. 2011), or the ecology of the system (e.g., the effect of crop condition on transition probabilities). However, they will not account for quadrat-level observation, such as misclassification of density states. How variation in the observation process and more subtle aspects of population ecology affect outcome in hierarchical cut-point formulations should be considered when fitting these models in future applications.

Density-structured modeling framework
This study provides evidence for the capability of hierarchical density-structured models for interrogating landscape-scale data sets on population responses to management. Analysis of asymptotic, transient and stochastic dynamics all consistently illustrated that populations responded to management. Populations that are spread over large areas are naturally subject to a wide range of environmental effects and, in this case, management practices, both of which are both likely to contribute to the variability in control that we observed here. The variability in these populations demonstrates the necessity of large-scale monitoring, as effective management relies on our understanding of the responses of populations to controls across the range of environmental contexts in which they exist. Crucially, our densitystructured approach provides the data and modeling tools that permit this.
An important recent conceptual advance has been the genesis of "landscape demography" (Gurevitch et al. 2016), which has demonstrated the importance of environmental heterogeneity for population dynamics, and the prevalence of scale dependencies in ecology (Hui et al. 2017, Quintana-Ascencio et al. 2018, Damschen et al. 2019. Although these approaches are still heavily reliant on intensive data collection, they provide detailed information on the causes and consequences of demographic variation across a landscape. Moreover, with the increase in accessibility of remote-sensing technology, it is becoming easier to collect expansive and high-resolution population-level data, which could begin to facilitate large-scale and in-depth studies of population FIG. 6. Contribution to the variance in the difference between initial and final densities by field ID/ matrix origin (first panel), initial density (second panel), and rotational control (third panel). Each point represents the size of the effect that variable had on the change in density. Points in the initial density panel are labelled accordingly (low, medium, high). Points in the rotation panel correspond to density-structured models of continuous barley, potatoes, beans, sugar beet, OSR, barley, continuous wheat, and peas, respectively. dynamics. Some studies already demonstrate the promise of remote-sensing data in large-scale spatiotemporal models (Conn et al. 2015, Tredennick et al. 2016. However, these studies rely on high quality population assessment, which is often difficult to achieve with remotely sensed data across large spatial extents (Lambert et al. 2018(Lambert et al. , 2019, and spatiotemporal models still remain computationally intensive compared to densitystructured approaches. Considering the importance of environmental heterogeneity that we have highlighted, the computational and data intensive nature of most large-scale modeling techniques will limit their utility to expedite management at the relevant scales. Densitystructured models are a useful component of the rapidly expanding set of tools in large-scale ecology. They offer a good alternative to many more complex approaches to fulfil the need for empirical studies of population dynamics across landscapes, managements, and environmental gradients. Although useful alternatives to traditional approaches, density-structured models may sometimes be unable to fully capture the underlying continuous dynamics of the population ). As such we need to better understand the limitations of density-structured models in order to bridge currents gaps in our knowledge of their ability to model a range of complex dynamics.

Management implications
Predicting the densities of weeds in the context of realworld environmental variability is vital to understanding the effects of management. Diversifying management options is necessary to maintain control of arable weeds, with increasing emphasis needed on non-chemical options (Chauvel et al. 2001. Interrogation of large-scale data sets on the ecology of arable weeds provides a route to improve our understanding of the responses of weeds to interventions at local and larger scales , Baucom and Busi 2019, Comont et al. 2019. Weed populations are inherently variable and there is a significant body of literature dedicated to understanding the causes and consequences of this variability (Freckleton and Watkinson 1998, 2002, Gonzalez-Andujar and Hughes 2000, Freckleton and Stephens 2009. Field-scale variation in weed density can have various origins. For example, persistent seedbanks, soil type, and local climatic variables will all substantially impact black-grass populations throughout the season (e.g., Colbach et al. 2006, Metcalfe et al. 2018). Monitoring and predicting highly variable dynamics is challenging, but the first step to tackling the issue is capturing the range of responses exhibited by a population. Our models capture the high inter-population variability across a national scale and within all the cropping systems we considered.
The considerable variation in black-grass density between fields, regardless of rotation used, has implications for both the large-scale modeling and management of black-grass. High levels of variability suggest that in some cases, any benefit from an applied control may be overwhelmed by location-specific effects positively influencing black-grass growth (Freckleton et al. 2017). Without the necessary data to model the drivers of the field-level effects of black-grass dynamics, it becomes extremely difficult to predict outcomes for individual populations. Additional covariates, such as climate (Lima et al. 2012, García De León et al. 2014, soil type (Colbach et al. 2006, Metcalfe et al. 2017, 2018, herbicide resistance status , 2011) and more specific management options (Holst et al. 2007, Harker and O'Donovan 2013, Metcalfe et al. 2017 are obvious areas where the predictive performance of these models could be improved, and should form the basis for future applications.
Despite high levels of variability, we have provided evidence for the effectiveness of rotation for reducing weed infestations. Weed densities were related to different aspects of crop rotation, with control varying by break crops and the proportion of wheat in a rotation. This result agrees with the current research on how cropping systems can provide control of problematic weeds (Zacharias and Grube 1984, Liebman and Dyck 1993, Chauvel et al. 2001, Melander et al. 2005, and demonstrates that this approach provides correct assessments of dynamics. In our models, cropping continuous wheat resulted in high density, highly variable infestations of black-grass that are likely to persist. Wheat is well known to be particularly susceptible to black-grass ): due to overlapping germination profiles, control is limited by the risk of damage to the crop itself (Thurston 1964). Our simulations demonstrated that many populations converged on medium densities, this suggests that alternative managements for moderate weed infestations may only produce small long-term changes, which emphasizes the importance of transient metrics in future modeling efforts. Although differences in weed observability in different crops may also contribute to some of the periodic differences we see across rotations, we show that rotation decreased not only the average density of black-grass, but also its variability and autocorrelation, FIG. 8. Mean density states, variance, and lag-1 autocorrelation in stochastic models with different break crops aggregated across levels of wheat dominance. Only models with a single break year and multiple levels of wheat dominance were included. Vertical bars represent two standard errors.
suggesting that weed populations will be more predictable and less likely to persist in rotated systems. Colonization of wheat fields via margins, farm machinery or crop seed are likely drivers of establishment of weed populations, but it is unlikely that they contribute to the differences in density we observe here, as black-grass typically sheds the majority of its seed before harvest and exists in low densities in field margins .
The benefits of rotational controls are widely appreciated in the literature and have various modes of action (Zacharias andGrube 1984, Liebman andDyck 1993). Primarily, rotation allows opportunities to apply controls without risking damage to crops. As black-grass emergence usually occurs during autumn (Thurston 1964, Moss 1990), rotating into a spring crop (known as spring cropping) is often cited as an effective control measure. Spring cropping can reduce black-grass abundance by facilitating targeted herbicide application, seed bed preparation, and cultivation during a period where the field is empty of crop, but during the germination period of the weed (Moss and Clarke 1994, Chauvel et al. 2001, Lutman et al. 2013). Many of the cropping systems with the lowest black-grass densities from our analyses include spring crops. For example, broad-leaf crops such as sugar beet, beans, and potatoes are generally planted in spring and are also resistant to grass-specific herbicides.
Control can also be achieved through direct and indirect competition for resources. Competitive cultivars such as barley or OSR can suppress weed populations through rapid accumulation of biomass and exclusion from nutrients and sunlight (Nicholas 1991). We see reductions in black-grass density from crops often cited by farmers as competitive, namely OSR and barley (Lutman et al. 2013). Compared to continuous wheat these crops showed noticeable reductions in density, but densities were generally higher than most alternatives. Some of the benefit of competitive cultivars, however, comes from resistance to yield penalties rather than reduction of seed return (Andrew et al. 2015), which may have accounted for the previous popularity of OSR despite the continued abundance of black-grass. An important consideration for future modeling is balancing the costs of infestation and controls. As wheat is the most valuable crop in the UK (DEFRA 2018), rotations that reduce the prevalence of wheat will reduce income. However, this economic loss will have to be balanced against the potential loss due continued infestation, the viability of alternative crops, and costs of additional controls.
An important limitation of density-structured models in the context of weed management is that they can only accurately provide summary descriptions of field-level densities. Spatial structure is an important driver of population dynamics and this is especially true for weeds. Numerous studies have investigated how the role of spatial structure in weed populations influences the persistence, spread, and control of weed populations (Rew and Cousens 2001, Holst et al. 2007, Metcalfe et al. 2017, Somerville et al. 2017, Gonzalez-Andujar et al. 2018. It is possible that the greater influence of fieldlevel effects compared to rotational management in our simulations is partly due to specific spatial configurations of weeds within fields. Within field dynamics, demography, and spatial structure are also important in dictating larger scale patterns of weed abundance (Cardina et al. 1997, Freckleton andWatkinson 2002a, b). Incorporating information on spatial interactions and plant life cycles will, therefore, be an important step for extending the utility of density-structured models and has two distinct advantages. First, it will improve the overall predictive capacity of density-structured models, allowing more accurate descriptions of landscape-scale dynamics. Second, it will allow predictions of individual populations and planning of targeted and more efficient management options.
Monitoring, understanding, and predicting large-scale population dynamics is an essential step for effective management. We demonstrate that the use of a rapid and accessible survey methodology in conjunction with density-structured modeling can directly tie empirical observations to predictions of landscape scale population responses to management and environmental heterogeneity. With these tools we show that the densities of an economically important weed were driven by rotational management, but all systems we studied were subject to high degrees of between-field variability in the degree of control. From this work, we have demonstrated that development of integrated controls that are effective over a wide range of environmental conditions will be necessary to limit the detrimental effects of arable weeds. Perhaps more importantly, this study demonstrates that we must extend the scale of sampling in ecology in order account for environmentally driven variability in dynamics. Density-structured models are a useful tool in this endeavor and offer a route to achieve robust and inexpensive analysis of spatially extensive populations.