Estimation and classiﬁcation of temporal trends to support integrated ecosystem assessment

We propose a trend estimation and classiﬁcation (TREC) approach to estimating dominant common trends among multivariate time series observations. Our methods are based on two statistical procedures that includes trend modelling and discriminant analysis for classifying similar trend (common trend) classes. We use simulations to evaluate the proposed approach and compare it with a relevant dynamic factor analysis in the time domain, which was recently proposed to estimate common trends in ﬁsheries time series. We apply the TREC approach to the multivariate short time series datasets investigated by the ICES integrated assessment working groups for the Norwegian Sea and the Barents Sea. The proposed approach is robust for application to short time series, and it directly identiﬁes and classiﬁes the dominant trends underlying observations. Based on the classiﬁed trend classes, we suggest that communication among stakeholders like marine managers, in-dustry representatives, non-governmental organizations, and governmental agencies can be enhanced by ﬁnding the common tendency between a biological community in a marine ecosystem and the environmental factors, as well as by the icons produced by generalizing common trend patterns.


Introduction
The integrated ecosystem assessment (IEA) is one approach to organizing scientific information at multiple scales and across sectors to support ecosystem-based fisheries management (EBFM) (Levin et al., 2009). IEA aims to analyse and synthesize information on a wide range of ecosystem components and pressures and to identify status, changes, relationships, and processes at the ecosystem level (WKINTRA report, ICES 2018). The outcome of an IEA can take multiple forms, which generally include descriptions of the main interacting ecosystem, human components, and past changes in these components. Eventually, it also gives an assessment of the risks associated with possible future trajectories of the ecosystem.
This information can then be fed back into the design of dedicated observational efforts (monitoring plans), the definition of multi-sectorial management objectives, or the establishment of new indicators and associated reference points (Levin et al., 2009(Levin et al., , 2014DePiper et al., 2017).
There is a long history of making quantitative assessments of individual fish stocks, and with it, common vocabularies and practices that support efficient communication between natural scientists, managers, the fishing industry, and other interested stakeholders (Hilborn and Walters, 1992;ICES, 2013). Even when the numerical methods used to reconstruct individual fish stock histories are complex, the outputs of fish stock assessments are communicated in a standardized manner that can generally V C International Council for the Exploration of the Sea 2020. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. be understood and used by different stakeholders. The same situation does not apply to IEAs, for which no methodological standards yet exist. Instead, many IEAs are still predominantly in the developmental phase, and a variety of methods have been explored. IEA results can reflect various aspects of an ecosystem beyond dynamics, status, and future risks, and they tend to be reported in an ad hoc fashion. For example, WKINTRA ICES (2018) reports using integrated trend analyses (ITA), which are used for IEA as a way to summarize changes that have occurred in recent decades in ecosystems in the north Atlantic and to highlight the possible connections among the physical, biological, and human ecosystem components. These methods cover graphical analyses as well as univariate and multivariate statistical analyses. Some methods have focused on individual components, while others have provided simplified representations of major past trends in the system or conveyed additional information on possible connections and causalities in the system. Unlike stock assessment models for which there are established benchmarking practices, no such practice exists for IEA. This can lead to ambiguous results or spurious interpretations, as was shown in the case of multivariate analysis applied to time series inappropriately (Kawasaki, 2004;Vanhatalo and Kulahci, 2016;Planque and Arneberg, 2018;Hallin et al., 2018) or in the case of a statistical method applied to a short time series dataset (Hardison et al., 2019;Solvang and Subbey, 2019). Therefore, to describe past changes in individual ecosystem components, more robust and theoretically correct trend analysis is required in IEA, and then the approach should be validated by benchmarking practices towards becoming a methodological standard.
Temporal changes in ecosystems can take the form of long-run movements-sometimes referred to as long-term trends or drifts-as well as short-/middle-run cyclic terms and noise components. The long-run changes entail a drift in time, a characteristic of temporal processes with a non-stationary mean (Kitagawa and Gersch, 1996). On the other hand, the cyclic and noise terms are considered stationary processes in the practical sense of time series analysis. When constructing statistical time series models, it is often useful to identify non-stationary trends and cyclic terms separately. Several structural time series models of this kind have been proposed in psychology and economics since the mid-1980s (Harvey, 1989;Kitagawa and Gersch, 1996;West and Harrison, 1997). When building multivariate time series models, it is also possible to decompose a non-stationary mean and cyclic components or seasonal components using a time series model (Kato et al., 1994(Kato et al., , 1996. Multivariate time series analysis has recently been used to model marine ecosystem dynamics (Solvang et al., 2017) and to investigate population's dynamics through causal inference (Solvang and Subby, 2018;Solvang and Subbey, 2019).
In this article, we focus on modelling the non-stationary mean trend (hereafter called "trend"), which represents the long-run movements in observations. In particular, we focus on finding similar trends, called Common trends, which refer to similar longterm tendencies across ecosystem components. Identifying common trends can be useful as a diagnostic tool to reveal past changes and to explore the relationships among biological communities and between these communities and environmental conditions. Individual and common trends are not observed directly; instead, these trends are assumed to exist and are represented as latent components in time series models. Such unobserved common components are termed "factors" (Helmut, 1993) in classical factor analysis. For the analysis of time series, dynamic factor analysis (DFA) has been applied to summarizing the information in macroeconomic analysis and to forecasting in a datarich environment (Darné et al., 2013;Ward et al., 2019). The factor models are described via a spectral approach in the frequency domain (Hallin and Lippi, 2013). This approach works well when the observed time series consists of a large number of observations that can be transferred from the time domain into the frequency domain using a Fourier transformation. However, data from marine ecosystem monitoring programmes often consist of relatively short time series, for which it is thus difficult to apply frequency-based methods. As an approach to facing this challenge, Zuur et al. (2014) suggested the use of DFA to estimate common trends in the time domain using structural time series modelling. The approach is implemented as a general modelling framework for state space representation in the MARSS library of the R language (Holmes et al., 2018). The state space form is a useful way to express several hidden components in the time domain, and the state is estimated by a Kalman filter or the expanded filtering theory (Shumway and Stoffer, 1982;Harvey and Pierse, 1984;Kitagawa and Gersch, 1996). Using DFA, Zuur et al. (2003) investigated the relationships between the estimated common trends for a biological community and environmental variables. In that work, an environmental variable was considered as an explanatory variable, and it was used for only the predicted biological responses that are fitted to observations. Zuur et al. (2003) used canonical correlations between the sea surface temperature series and common biological trends to explore possible relationships. In a canonical correlation, variations in a multivariate system are assumed to follow a stationary Gaussian process, although some formulations have proposed correcting for certain departures from the Gaussian assumption (Min and Tsay, 2005). In the case of biological and environmental variables in large marine ecosystems, it is reasonable to assume that multidecadal time series are a non-stationary process, not a stationary Gaussian process, in response to ongoing climate change and the associated biological responses.
For such non-stationary mean time series data consisting of environmental variables and biological communities, we propose a more direct approach to estimating individual trends, classifying these trends, and extracting a set of underlying common trends. We term this method TREC, for trend estimation and classification. The non-stationary mean time series are modelled using polynomial trend models (Solvang et al., 2008) or a stochastic trend model by state space representation, commonly used as detrending models in time series analysis (Kitagawa and Gersch, 1996). Then, classification for common configurations of the trends is performed using discriminant analysis (Solvang et al., 2008). The common configurations are further subdivided into sets of underlying common trends based on the target trends selected from the estimated trends. The common trends identified for each class are representative of dominant non-stationary patterns for these classes and are interpreted as common variations of biological and environmental data. TREC can be applied to short time series data, and it is not necessary to apply a stationary Gaussian assumption to the estimated trends to investigate the relationships among them. We expect that the results can be communicated in an easily accessible form to serve the needs of multiple stakeholders. We evaluate the proposed approach using simulated datasets, including short and long time series, and apply the method to multidimensional time series observations from the ICES IEA working group for the Barents Sea (WGIBAR, ICES 2016) and the Norwegian Sea (WGINOR, ICES, 2015).

Trend models
The observation model of a time series is given by (1) where tðnÞ is the trend component and uðnÞ is the residual component at time step n. In this article, we apply two different kinds of parametric trend models: a polynomial trend model using a polynomial regression model (Solvang et al., 2008) and a stochastic trend model using a dth order difference equation model (Kato et al., 1994(Kato et al., , 1996. The polynomial trend model was combined with the discriminant trend analysis examined in Solvang et al. (2008), given by The apostrophe symbol in vector b denotes transposition. The least squares estimator for b is given byb LSE ðZ0ZÞ À1 ZY; (2) where Z ¼ 1; 1; 1 2 ; Á Á Á; 1 pÀ1 1; 2; 2 2 ; Á Á Á; 2 pÀ1 . . .
The trend component is estimated by which is defined on all time steps f1; 2; Á Á Á; N g. Assuming that the residual uðnÞ obeys a normal distribution with variance r 2 , the estimated variance is given bŷ yðnÞ Àb 0 Àb 1 n À Á Á Á Àb pÀ1 n pÀ1 2 yðnÞ Àt ðnÞ 2 : The log-likelihood of the polynomial trend model for this is given by The stochastic trend model is defined by the dth order difference equation, which was posed as the smoothing problem by Whittaker and Robinson (1924). This model allows for more flexible trends than does the polynomial regression model. The stochastic trend model is expressed in the following way: where r is a difference operator rtðnÞ ¼ tðnÞ À tðn À 1Þ and v t ðnÞ is assumed to be a white noise sequence. If d ¼ 1, tðnÞ % tðn À 1Þ and the trend is known as a random walk model. If d ¼ 2, tðnÞ À 2tðn À 1Þ þ tðn À 2Þ % 0 (Kitagawa and Gersch, 1996). Provided that the variance of v t ðnÞ is sufficiently small, tðnÞ yields a smooth trend. The model can be represented in state space form as yðnÞ ¼ HzðnÞ þ wðnÞ; where zðnÞ is the state vector corresponding to tðnÞ, vðnÞ is the system noise vector with mean 0 and unknown variance r 2 v , F, G, and H indicate integer or matrices, and wðnÞ is observation error with mean 0 and unknown variance r 2 w . The state zðnÞ is taken as a latent factor, since we cannot directly observe it from the data. Corresponding to the above dth order difference equation, when d ¼ 1, When d ¼ 2, the state vector and matrices are as follows: The observation error corresponds to uðnÞ in (1) in this case. The trend component is estimated using a Kalman filter, which is a powerful numerical algorithm that recursively operates the state estimation, prediction and filtering (Kato et al., 1994(Kato et al., , 1996: Prediction: Filtering: Here, zðnjn À 1Þ and V ðnjn À 1Þ are the conditional mean and conditional variance, R is the observation error, and K is Kalman gain. This trend model includes the parameter vector h ¼ ðd; r 2 v ; r 2 w Þ. The log-likelihood function lðhÞ of this model is given by where Y ðn À 1Þ ¼ ðyð1Þ; yð2Þ; Á Á Á; yðn À 1ÞÞ, DyðnÞ ¼ yðnÞÀ Hzðnjn À 1Þ, and RðnÞ ¼ HðnÞV ðnjn À 1ÞH 0 ðnÞ þ r 2 w . The flexibility of the estimated trend depends on r 2 v . The optimum r 2 v can be determined by maximum likelihood within an arbitry variance range. The variance r 2 w can be directly set to the variance of the observation. Denoting the number of parameters by c, the Akaike Information Criterion (AIC) (Akaike, 1974) value of the model, given the estimate of optimum orders p for the polynomial trend model and d for the stochastic trend model, is The log-likelihood function can also be represented by Laplace approximation, and the optimum parameters can be estimated by a numerical optimization tool such as the AD Model Builder package (Fournier et al., 2012) or TMB (Kristensen et al., 2016).

Discrimination analysis of trends
For the estimated trend, we apply discriminant analysis, as proposed by Solvang et al. (2008). First, let us consider the following simple example of two-category discrimination: P 1 : model ð1Þ with trend function T 1 ðnÞ P 2 : model ð1Þ with trend function T 2 ðnÞ: Suppose that we observe new data and that the trend is estimated ast ðnÞ, which is assumed to belong to P 1 or P 2 . The classification is performed with the following divergence measure: If Lðt : T 2 Þ > Lðt : T 1 Þ, the retained category is P 1 , otherwise it is P 2 . The probability of misclassifying the observation from P i into P j ði 6 ¼ jÞ converges to zero as n ! 1. We define the discriminate function D Lðt : T 2 Þ À Lðt : T 1 Þ for use as a discriminant score (Solvang et al., 2008). In practical use, we fix two reference trends corresponding to T 1 and T 2 and make the respective trend estimatorsT 1 ðnÞ andT 2 ðnÞ using model (1). For common trend classification, it is considered in practice that T 1 ðnÞ andT 2 ðnÞ represent different (opposite) shapes, namely increasing and decreasing. In other words, we assume that we can obtain a rising tendency and a declining tendency for annual changes. The different shapes result in a large distance between the two categories. Then, the discriminant function for observation y ðjÞ ðnÞ, j ¼ 1 to 2, is defined bŷ Àt j ðnÞg 2 À fT 1 ðnÞ Àt j ðnÞg 2 : IfD j > 0, category P 1 is chosen, otherwise, category P 2 is chosen. Applying the nearest neighbour method, we can classify groups according to similarD j , e.g. upward, downward, and flat, including convex or concave. In this article, we first apply a twocategory discriminant analysis to the estimated trends to roughly divide them into three groups representing configurations related to upward, flat, and downward. If it is necessary to classify them into groups of more concrete patterns from the three roughconfiguration groups, this two-category classification can be easily extended to multiple-category discrimination. The problem is then specified in the following way: P j : model ð1Þ with trend function T j ðnÞ; j ¼ 1; 2; Á Á Á; k: The divergence measure in this case is given by We provide a divergence measurement vector given by The classification rule is defined as the requirement that the estimated trendt belong to P l to satisfy Lðt : T l Þ ¼ minðfÞ: In practice, the reference trend T j should be predefined in the three configurations groups obtained by the two-category discriminates. Finally, the reference trend T j is assigned as a general reference, called an icon, which is an easily accessible form that can be used to serve the needs of stakeholders.
The entire numerical procedure and the icons we predefine in this study are summarized in Figure 1. TREC is implemented using the MATLAB code (MATLAB ver. R2018b), which is available on request.

Common trend approach by dynamic factor model
As mentioned in the Introduction, DFA is an appropriate method for analysing time series data and investigating the common trends across ecosystem components. Therefore, we refer to the approach provided by Zuur et al. (2003) and compare it with TREC. This comparison is considered in detail in the Supplementary Text.
We then derive six variables as follows: x 1 ðnÞ ¼ 0:3 Â t 1 ðnÞ þ r 1 ðnÞ; x 2 ðnÞ ¼ t 1 ðnÞ þ x e ðnÞ þ r 2 ðnÞ; x 3 ðnÞ ¼ À0:5 Â t 1 ðnÞ þ r 3 ðnÞ; x 4 ðnÞ ¼ 0:0004 Â t 2 ðnÞ þ 0:5 Â x e ðnÞ þ r 4 ðnÞ; x 5 ðnÞ ¼ 0:00045 Â t 2 ðnÞ þ r 5 ðnÞ; x 6 ðnÞ ¼ 1:5 Â x e ðnÞ þ r 6 ðnÞ; where r 1 ðnÞ; Á Á Á; r 6 ðnÞ are random variables normally distributed with zero mean and unit variance. We simulated series with lengths of 50-and 200-time steps (Supplementary Figure S2). We first applied the polynomial and stochastic trend models to x 1 ðnÞ; Á Á Á; x 6 ðnÞ and x e ðnÞ. Supplementary Table S1a summarizes the calculated log-likelihood and AIC obtained by applying the stochastic trend model for d ¼ 1 and d ¼ 2. The results suggest that the model fits best for d ¼ 1. The optimum order p (based on AIC) for the polynomial trend model and the optimum variance Q for the stochastic trend model are summarized in Supplementary Table S1b. For optimizing Q, we set 0:01 Q 0:1 as a search range. This procedure was conducted for each time series, and it resulted in a set of estimated polynomial trends (Supplementary Figure S3a and b, red lines) and stochastic trends (Supplementary Figure S3c and d, red lines). The filtering of state vector zðnÞ in the stochastic model (Supplementary Figure S3c and d, blue lines) was generally more variable than the associated polynomial trend. This is because the state vector can follow data fluctuations according to each time point by using a Kalman filter. The prediction values usually indicated more fluctuation than the filtering values.
Next, we performed a two-category discrimination analysis. The estimated trends for x 2 and x 3 were defined as the reference trendsT 1 ðnÞ andT 2 ðnÞ to calculate the discriminant function. The value obtained by the function was then used as a distance metric to perform hierarchical clustering, applying the unweighted centre of mass distance as the linkage between clusters. In Supplementary Figure S4, the dendrograms for polynomial (a and b) and stochastic trends (c and d) indicate a discrimination between upward, flat, and downward trends, which are coloured by red, blue, and black thick lines, respectively. Based on the reference trends, two categories for downward and upward were classified by D > 0 or D < 0. Some D values around 0 belonged to the flat group. In the case of data using 50 time points (Supplementary Figure S4a and c), it might be appropriate to classify the polynomial trend x 5 ðnÞ into the upward configuration group rather than to the flat configuration group, while the stochastic trend of Verification for discriminant function and divergence measure for j > 2 Solvang et al. (2008) mentioned the possibility that several common configuration groups were classified by the discriminant function for only the two-category setting. However, if the estimated trends present an exactly symmetrical configuration, e.g. convex and concave, both trends may be appropriate to deal with classifying in the flat condition. We demonstrated such a case using the artificial trend series generated by quadratic and exponential functions shown in Supplementary Figure S5. By selecting number 33 and 30 trends as reference trendsT 1 ðnÞ andT 2 ðnÞ, the two-category discriminant was considered. The bar plots for D j are illustrated in Supplementary Figure S6. The bars in groups for upward and downward correspond to positive and negativê D j values, that is, all data can be divided into two main categories.D j corresponding to patterns such as convex or concave is clearly indicated around 0 (Supplementary Figure S6, trends 1 -20). This can be resolved using multi-category discrimination, with j > 2, as defined in previous works (18)-(20). In this case, eight-category discrimination was performed, with trends 33, 30, 1, 2, 29, 21, 38, and 32 set as the reference trends. Applying this eight-category discrimination, all data were classified into certain categories that include similar configurations as shown in Supplementary Table S2. The predefined icons set in Section 2 were also assigned according to the reference trend's patterns.

Field observations
We applied the proposed method to the time series compiled by the ICES integrated assessment working groups for the Barents Sea (WGIBAR, ICES, 2016) and the Norwegian Sea (WGINOR, ICES, 2015). Time series lengths vary among variables, so we selected the periods with continuous records (no missing data) until the final year of observation for all biotic and abiotic data. The WGIBAR dataset consists of 35 annual time series for the period 1980-2015, including 7 abiotic, 18 biotic, and 8 human impact variables. The WGINOR dataset consists of 24 annual time series     Table 1, including the data that are missing less than three time points.

Results and discussion
The stochastic trend models with d ¼ 1 outperformed the stochastic trend model for d ¼ 2, for both datasets, as indicated by the AIC values (Supplementary Figure S7a for WGIBAR and Supplementary Figure S7b for WGINOR). The estimated polynomial and stochastic trends visible in Figure 2 (WGIBAR) and Figure 3 (WGINOR) and the optimum numbers of parameters are summarized in Supplementary Table S4. While the order of polynomial trend models can vary among individual time series, the values of optimum Q are mainly equal to 0.1, reflecting the fact that flexible stochastic trends are predominantly selected. The clustering of trends is visualized in Supplementary Figures S8 (WGIBAR) and S9 (WGINOR). The classification procedure via discriminant analysis requires the pre-selection of reference trends. As mentioned, the reference trends for the two-category discriminant should be clearly opposite shapes, such as upward and downward. For WGIBAR, trends 16 (number of lumpfish age 1þ) and 30 (landings of harp seals) were used as references for both polynomial and stochastic trend models. For WGINOR, trends 5 (Maximum chlorophyll a level in Norwegian basin) and 9 (recruitment of mackerel) were used as references for polynomial trend models, and trends 2 (sub-polar gyre index from satellite ssh data) and 18 (puffin stock size) were used as references for stochastic trend models.
The trends were first divided into two categories and further classified into three categories: upward, flat, and downward. Tables 2 and 3 list the abiotic and biotic data for the common configurations. For most individual trends, the classification into upward, flat, and downward groups was consistent between the polynomial and the stochastic trend models. A few trends were exceptions to this general pattern: while in the WGIBAR analysis polynomial trend 15 (ages 1 and 2 herring stock biomass) was convex and polynomial trend 31 (Landings of minke whales) was concave, the corresponding stochastic trends 15 and 31 display more fluctuations and are not discriminated into flat configurations. Similar tendencies were observed in trends 8 (recruitment of blue whiting per year class at age 1), 15 (weight at age 6 in herring stock), and 16 (weight at age 6 in blue whiting catch) for WGINOR. The simulation study in the case of 50 time points also indicated inconsistent classification in Supplementary Figure  S4a and b. This suggests that a polynomial trend could be handled in a simpler way to assign trend estimates to several icon configurations using multi-category discrimination.
Using polynomial trend estimates, we set more reference trends in the three groups and classified all trends precisely into groups by multi-category discrimination. The subdivided outputs in three common trends are summarized in Figure 4 for WGIBAR and Figure 5 for WGINOR. Tables 4 and 5 list the details for classified data in each category. To each classified trend, we associated an icon representation that could be used to simply and efficiently communicate information on long-term changes in variables. The estimated trends by the polynomial model present rough long-term changes and may be more useful for assigning icons. On the other hand, the estimated trend by the stochastic trend model can follow more precise changes according to data fluctuations. Therefore, analysers could select the polynomial or stochastic trend model to investigate long-term changes in their data depending on their aim or interest.
The icons presented in Table 4 highlight the observation that in the Barents Sea there has been a continuous increase in temperature. It is visible in the increasing temperature trends in the Bear Island trough (4) and the Kola section (6) and in the increasing extent of the area occupied by Atlantic water (2) as well as the decreasing trend in the extent of area occupied by Arctic water (1). These changes have been paralleled by increases in the biomass of jellyfishes (8), the abundance of juvenile lumpsuckers (17), and haddock recruitment (24) and by decreases in fishing mortality of shrimps (29) and landings of shrimps (33) and harp seals (30). The increasing trends in krill biomass (9) and Table 2. Results of discrimination analysis for WGIBAR data: three common configuration groups of trends by two-category discriminates and the assigned data.  Table 3. Results of discrimination analysis for WGINOR data: three common configuration groups of trends by two-category discriminates and the assigned data.  (27) has been accelerating in recent years. A similar understanding of trends in the main abiotic, biotic, and human factors in the Norwegian Sea can be achieved based on the icons presented in Table 5. As seen by these interpretations of the outputs, these icons can be used in ecological "dashboards" displaying a list providing summaries of ecological states and trends.  The combination of trend identification, classification, and icon representations provides an easily accessible representation of the trends in abiotic, biotic, and human factors, which can support discursive interpretations while still being rooted in an objective statistical approach.
In reality, the majority of marine ecosystems do not have complete datasets for a full time series to be analysed. For such missing data, the provided trend models support he interpolation of them. In the case of a polynomial trend model, Solvang et al. (2008) showed the interpolated trend estimates and the classification. Furthermore, the stochastic trend model can fill the missing data in the Kalman filter recursions (Kitagawa and Gersch, 1996).
Considering the meaning of common trend, the dominant trend's configuration presents long-run movements towards comprehensive data and it may be useful for a biologist or ecologist if the fluctuations in environmental factors and fish communities could be easily interpreted as upward, flat, or downward. In the case of analysis by the DF model (Zuur et al., 2014) summarized in Supplementary Text, the obtained common factors include several patterns mixing long-term trend and cyclic fluctuations, and thus, it may be difficult to interpret the physical meaning of such factors. Since most relevant observations have a small sample size of data points in the time series data, comparing configurations may provide a simple and useful way to gain a preliminary understanding of an ecosystem. TREC helps to provide a simple interpretation of which trend pattern is dominant over all relevant data for abiotic and biotic cases. When we predefined the reference trends used for discrimination, no prior knowledge was considered, that is, it was simply done in an arbitrary manner in this article. However, it would also be possible to predefine reference trends based on the prior knowledge of marine biologists or ecologists.

Conclusion
TREC is proposed as an approach to analysing common trends in a marine ecosystem. It consists of two procedures: (i) estimating trends with statistical trend models and (ii) classifying the estimated trends into categories of common configurations. The classification step enables us to find specific configurations by representing the estimated trend patterns. A simulation clarified the performances of two different trend models and their flexibility for use with several representative patterns, which can be predefined as icons. We applied TREC to two real time series datasets provided by the ICES integrated assessment working groups. Abiotic, biotic, and human impact data were classified into common trend groups. The proposed TREC focuses on long-term trends of data, and it works for any length of time steps. TREC could become a methodology established for ITA in IEA by validation through benchmarking practices. Studies conducted to investigate precise ecosystem functions using TREC are expected as further extensions of this work.

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.

Funding
This work has been funded by the Research Council of Norway, titled "Barents-RISK; Assessing risks of cumulative impacts on the Barents Sea ecosystem and its services" (Project number # 288192), and the Project "causality and food web modelling in the Barents Sea" (# 14565) from the Institute of Marine Research, Norway. Table 4. Results of discrimination analysis for WGIBAR data: assigned icons by multi-category discriminates. Table 5. Results of discrimination analysis for WGINOR data: assigned icons by multi-category discriminates.