Clement Aldebert
email: clement.aldebert@mio.osupytheas.fr

Daniel B Stouffer

November

Community dynamics and sensitivity to model structure: toward a probabilistic view of process-based model predictions

Keywords: model predictions, Bayesian statistics, structural sensitivity, community dynamics, predation

Introduction

With the need for more accurate predictions in biology and environmental sciences [START_REF] Tr Anderson | Progress in marine ecosystem modelling and the "unreasonable effectiveness of mathematics[END_REF][START_REF] Ol Petchey | The ecological forecast horizon, and examples of its uses and determinants[END_REF][START_REF] Pennekamp | The practice of prediction: What can ecologists learn from applied, ecology-related fields?[END_REF], two philosophically different streams of research have been growing, statistical inference and mechanistic modelling. While the former aims to make predictions based on uncovering statistical relationships in large data sets, mechanistic modelling aims to make predictions based on causal mechanisms that explain observed patterns. In practice, the pros of one approach are the cons of the other, so a promising way forward would be to combine them in a "symbiotic relationship" [START_REF] Baker | Mechanistic models versus machine learning, a fight worth fighting for the biological community?[END_REF]. Here, we provide an example of such cross-fertilisation. Specifically, we use Bayesian statistics to present probabilistic predictions of a deterministic mechanistic model-built around empirical data-in way that takes into account uncertainty both in model construction and model parameterization. We therefore improve the model's predictive capability by including prediction uncertainty while maintaining the explanatory power of a mechanistic model. As our focal model, we use the Rosenzweig and MacArthur [START_REF] Ml Rosenzweig | Graphical representation and stability conditions of predator-prey interaction[END_REF] predator-prey model, which is known to be "structurally sensitive" [as defined by 6]; that is, apparently minor changes in model formulation can lead to a dramatic change in both quantitative (predicted biomasses over time) and qualitative predictions, such as prey-predator oscillations or the coexistence of alternative stable states [START_REF] Mr Myerscough | Stability, persistence and structural stability in a classical predator-prey model[END_REF][START_REF] Gf Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF][START_REF] Seo | Sensitivity of the dynamics of the general rosenzweigmacarthur model to the mathematical form of the functional response: a bifurcation theory approach[END_REF].

Structural sensitivity is a common phenomenon that emerges in mechanistic biological models, which usually aim to summarise multi-level processes into equations after adopting simple assumptions regarding the complexity of the biological system of interest. As the entire complexity can rarely, if ever, be taken into account, available empirical data may be insufficient to statistically discriminate between alternative models [START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF][START_REF] Jassby | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF].

Alternative models are known to make different predictions when the uncertain process is the infection in a host-pathogen system [START_REF] Wood | Super-sensitivity to structure in biological models[END_REF], the colimited uptake of nutrient [START_REF] Jc Poggiale | How far details are important in ecosystem modelling: the case of multi-limiting nutrients in phytoplankton-zooplankton interactions[END_REF], or predation in predator-prey and food-web models [START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF][START_REF] Mr Myerscough | Stability, persistence and structural stability in a classical predator-prey model[END_REF][START_REF] Gf Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF][START_REF] Gross | Enrichment and foodchain stability: the impact of different forms of predator-prey interaction[END_REF][START_REF] Mw Adamson | When can we trust our model predictions? unearthing structural sensitivity in biological systems[END_REF][START_REF] Mw Adamson | Bifurcation analysis of models with uncertain function specification: how should we proceed[END_REF][START_REF] Mw Adamson | Defining and detecting structural sensitivity in biological models: developing a new framework[END_REF][START_REF] Aldebert | Does structural sensitivity alter complexity-stability relationships?[END_REF][START_REF] Aldebert | Structural sensitivity and resilience in a predatorprey model with density-dependent mortality[END_REF][START_REF] Aldebert | Is structural sensitivity a problem of oversimplified biological models? insights from nested dynamic energy budget models[END_REF][START_REF] Aldebert | Three-dimensional bifurcation analysis of a predator-prey model with uncertain formulation[END_REF]. Some studies also indicate that major ocean-scale predictions, such as the dominance of phytoplankton groups, primary production and export to the deep ocean, can be deeply affected by this form of sensitivity [START_REF] Tr Anderson | Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model[END_REF][START_REF] Tr Anderson | Sensitivity of secondary production and export flux to choice of trophic transfer formulation in marine ecosystem models[END_REF][START_REF] Vallina | Maximal feeding with active prey-switching: a kill-the-winner functional response and its effect on global diversity and biogeography[END_REF]. In a way, structural sensitivity can be considered an extension of classical parameter sensitivity; moreover, Cordoleani et al. [START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF] and Adamson and Morozov [START_REF] Mw Adamson | When can we trust our model predictions? unearthing structural sensitivity in biological systems[END_REF] have provided examples where a change in model formulation has a higher effect on model predictions than an equivalent change in parameter values. However, as far as we know, such an analysis has never been statistically performed simultaneously across alternative models and all plausible parameter values.

Here, we ask whether the predictions made by a biological model are more sensitive to the uncertainty in its parameterization or to its mathematical formulation. To answer this question, we present a probabilistic view of predictions made by a deterministic predator-prey model, the Rosenzweig and MacArthur [START_REF] Ml Rosenzweig | Graphical representation and stability conditions of predator-prey interaction[END_REF] model (section 2). In particular, we explore the consequences of alternative predator functional responses in changing the behaviours predicted. We achieve the shift from determinism (section 3) to probabilism (section 4) by merging bifurcation theory for model analysis [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF][START_REF] Perko | Differential equations and dynamical systems[END_REF][START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF] and Bayesian statistics [START_REF] Claeskens | Model selection and model averaging[END_REF][START_REF] Parent | Introduction to Hierarchical Bayesian Modelling for ecological data[END_REF][START_REF] Mcelreath | Statistical rethinking: a Bayesian course with examples in R and Stan[END_REF] by harnessing the latter's ability to include uncertainty in a model and to propagate it forward into predictions. Intriguingly, our probabilistic approach -applied on three example data sets-indicates that parametric uncertainty is regularly overcome by model uncertainty (section 5), an observation that has broad implications across various other challenges faced in biology.

A simple predator-prey model

We study the Rosenzweig-MacArthur model [START_REF] Ml Rosenzweig | Graphical representation and stability conditions of predator-prey interaction[END_REF]:

               dN dτ = rN   1 - N K    -f (N )Y dY dτ = (ef (N ) -µ) Y. (1)
Here, the prey population N follows logistic growth with per capita growth rate r and environmental carrying capacity K in the absence of predators. The predator population Y has a linear per capita mortality rate µ, conversion efficiency e, and functional response f (N ). All parameters and variables are strictly positive to ensure their biological meaning.

In line with the biology, a well-behaved functional response is expected to be strictly increasing, concave, saturating, and vanishes only at 0. Here, we consider three possible alternatives for

f ∈ F = {f (H) , f (I) , f (t) },
all of which only depend upon two parameters:

f (H) (N ) = a (H) N 1 + a (H) h (H) N , f (I) (N ) = 1 h (I) 1 -e -a (I) h (I) N , f (t) (N ) = 1 h (t) tanh a (t) h (t) N , (2)
which are the Holling Type II [30, 31, later denoted Holling], Ivlev [START_REF] Vs Ivlev | Experimental ecology of the feeding of fishes[END_REF], and the Hyperbolic Tangent [START_REF] Jassby | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF],

respectively. The parameters in each have the same mathematical meaning: 1/h (•) corresponds to the maximum uptake rate and a (•) corresponds to the function's slope in the limit of no prey, but they originate from different perspectives. For Holling, a (H) and h (H) are predator attack rate and handling time. For Ivlev, 1/h (I) and a (I) h (I) are predator maximum digestion rate and satiation coefficient. The Hyperbolic Tangent is actually a phenomenological model without underlying biological assumptions, which nevertheless can at times provide a more accurate description of data than others [START_REF] Jassby | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF].

To drop one parameter and ease the analysis without loss of information, we re-scale model ( 1): y = Y /r, t = rτ , m = µ/e, ε = e/r, and write x = N . The re-scaled model reads:

               dx dt = x   1 - x K    -f (x)y dy dt = ε (f (x) -m) y. (3)
Upon re-scaling, the remaining parameters are (i) the prey carrying capacity K, (ii) the scaled predator mortality m, (iii) the time-scale factor ε, and (iv) the two parameters a (•) and h (•) from the functional response that must be estimated from data. These first two are of biological interest as they might be affected by external factors (e.g. environmental degradation, additional predator mortality due to harvesting).

3 Probabilistic predictions

Fitting functional responses on data

Empirical data were not considered in previous structural sensitivity analyses in Rozensweig-MacArthur model [START_REF] Gf Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF][START_REF] Seo | Sensitivity of the dynamics of the general rosenzweigmacarthur model to the mathematical form of the functional response: a bifurcation theory approach[END_REF]. As our baseline, we therefore use three examples of functional-response data to provide an overview of our approach and its practical use. To introduce the approach, we first focus on one of these data sets (until section 5). Specifically, we use data from an experiment where two individuals of Gazella thomsoni were eating hand-assembled grass swards, and their grazing rate as a function of grass biomass (i.e. a functional response) was estimated [33, figure 1]. This consumer-resource system can be modeled with generic assumptions on population dynamics: (i) grass biomass follows a logistic growth limited for instance by space and nutrients availability;

(ii) the growth rate of Gazella thomsoni population is proportional to its grass intake which is a strictly increasing, concave, saturating function of grass biomass that vanishes in absence of grass; (iii) Gazella thomsoni has a linear per capita mortality rate (e.g. ageing, harvesting). All these assumptions define Rozensweig-MacArthur model, which is often refered as a predator-prey model but is generic enough to also describe a variety of consumer-resource systems.

Previous studies on structural sensitivity compared a deterministic analysis of model predictions based on best-fitted functions [8, 9, 11, 18, 20, among others]. Fitting a function to data implies uncertainty in the parameters being inferred, which means that there is additional insight to be gained by considering both parameter and model uncertainty. To perform a probabilistic analysis of model predictions, we first estimate the parametric uncertainty while fitting each functional response to data. This is achieved by a Hamiltonian Markov Chain Monte Carlo (HMCMC) algorithm computing the probability P (θ f ) that the set of functionalresponse parameters θ f = (a (f ) , h (f ) ) of function f predict observed data (figure 1a-c), assuming a Gaussian noise around f [29, details in Appendix]. The probability P (θ f ) is proportional to the likelihood of θ f , and the set of probabilities corresponding to the inferred parameter values is called the posterior distribution in Bayesian statistics. This posterior distribution of θ f is estimated by the HMCMC algorithm by learning from available data (figure 1a-c). The maximum of this distribution occurs at parameters having the maximum likelihood and thus giving the best fit to data based on this criterion. Therefore, the HMCMC approach gives the same best-fitted functions as a maximum likelihood estimation, but it adds information about parametric uncertainty for each function by considering the whole posterior probability distribution (figure 1d-f).

Based on the posterior distribution for each function, we use the Widely-Applicable-Information-Criterion (WAIC) to derive Akaike weights which correspond to the relative predictive accuracy of each function -i.e. the probability P (f ) that function f ∈ F gives the best fit to new data, conditional to the alternative functions that we consider [START_REF] Mcelreath | Statistical rethinking: a Bayesian course with examples in R and Stan[END_REF]. This relative accuracy can be estimated using any other suitable method than WAIC, without altering our generic framework. For the first data set, the Hyperbolic Tangent has the highest probability of (d-f) Functional responses fit to experimental data (points '+'). The parametric uncertainty from the HMCMC estimation gives a confidence interval (95 %, shaded area) around the best-fitted functions (curves). Model uncertainty is derived through the relative likelihood that one function fits new data better than the others, knowing their respective parametric uncertainty. 1). Based on that, one could decide to focus only on the Hyperbolic Tangent as it seems to be the best function. However, doing this neglects the fact that the alternative functions can be a better description of new data with a significant probability (0.428). Our approach hence aims to keep all the alternative functions, and to merge their predictions based on their respective likelihood through the use of model averaging.

Overview of model predictions

Let use first present an overview of model predictions for a given functional response f and its plausible parameter values θ f . We performed a bifurcation analysis based on analytical and numerical results (mathematical details in Appendix) that reveals all the qualitative asymptotic dynamics predicted by the model for any environmental condition (K, m). We limited the range of these conditions to K ∈]0, x max ], where

x max ≈ 197 is the maximal prey abundance in functional response data, and m ∈]0, m max ], where m max :=

f ∈F P (f ) R 2 + P (θ f , D)f (x max )dθ f ≈ 5.
44 is the average mortality rate (among all functional responses) that allows predator survival for the considered range of prey abundance. The choice of m max restrain our analysis to the range of parameter values that will give interesting results. Any m > m max will lead to predator extinction in most cases, and this is a trivial result that presents no interest for our study since all models necessarily give identical predictions. 2) , y (2) ) exists. It starts to exist for higher carrying capacity values when predator mortality (i.e. the loss that must be overcome to ensure predator survival) is higher. This equilibrium is stable only in the blue and red areas. In the green area found at high carrying capacity, the paradox of enrichment [START_REF] Ml Rosenzweig | Paradox of enrichment: destablilization of exploitation ecosystems in ecological time[END_REF] has destabilised the coexistence equilibrium and there are stable prey-predator oscillations. Stable oscillations and the stable coexistence equilibrium are both possible in the red area, which is only predicted with the Hyperbolic Tangent at low predator mortality and high carrying capacity. In this bistability area, the model predicts that the system will converge to one of the two alternatives depending on its initial state. The presence of alternative stable states is of particular interest to study ecosystem resilience when facing external disturbances [START_REF] Aldebert | Structural sensitivity and resilience in a predatorprey model with density-dependent mortality[END_REF][START_REF] Scheffer | Anticipating critical transitions[END_REF]. This other form of the paradox of enrichment has been found in predator-prey models incorporating a more detailed description of organisms biology [START_REF] Aldebert | Is structural sensitivity a problem of oversimplified biological models? insights from nested dynamic energy budget models[END_REF].

E (2) = (x (
All of the above qualitative predictions are independent of the time scale ε, except for the bistability area that we computed numerically for ε = 1; this area starts at lower (higher) K for higher (lower) ε, which does not affect our conclusions.

As a general conclusion, using Holling and Ivlev functional responses lead to similar patterns of predictions and always one stable state whereas using the Hyperbolic Tangent allows us to predict bistability. Notably, the transitions between different qualitative predictions occur in different regions of parameter space for different models and correspond to the following biological phenomena: predator invasion (transcritical bifurcation: extinction -coexistence), onset of oscillations through enrichment paradox (supercritical Hopf bifurcation: coexistence -oscillations), and potential catastrophic shifts (i.e. tipping points) with a change in the number of alternative stable states (subcritical Hopf bifurcation: bistability -oscillations ; and limit point of cycles bifurcation: bistability -coexistence). Note that if the number of alternative stable states is affected by structural sensitivity, it becomes hard to estimate ecosystem resilience under disturbances [START_REF] Aldebert | Structural sensitivity and resilience in a predatorprey model with density-dependent mortality[END_REF].

Introducing parametric uncertainty into predictions

To introduce parametric uncertainty in our analysis, we now look at the probability that model (3) with function f predicts the different qualitative dynamics X, which can be any of "extinction", "equilibrium", "oscillations", or "bistability". For fixed values of the model parameters α = (K, m, ) that do not relate to f , the probability that functional-response f leads to the prediction X is:

P (X|f ) = R 2 + M (X, α, f, θ f )P (θ f )dθ f , (4)
where M (X, α, f, θ f ) = 1 if the model based on function f and with parameter values (θ f , α) predicts dynamics X and 0 otherwise, and the integral is over the range of the posterior distributions of the two parameters in θ f . In practice, this probability ( 4) is estimated by performing model analysis for a finite sample of θ f values (here 1000 parameter sets) drawn randomly from the posterior distribution P (θ f ) [START_REF] Mcelreath | Statistical rethinking: a Bayesian course with examples in R and Stan[END_REF]. Computing the probability (4) for different model parameters α therefore provides a probabilistic bifurcation analysis of model predictions based on one functional response. Translating the probabilities of each predictions into colour gradients gives the "fuzzy" transitions between qualitative dynamics (figure 2a-c

Introducing model uncertainty into predictions

To go further, model uncertainty can be introduced through a weighted mean over the alternative functions.

This allows us to compute the probability P (X) that the model generally predicts the different qualitative dynamics X:

P (X) = f ∈F P (f )P (X|f ). (5)
The resulting probabilistic predictions are shown in figure 2d. Not surprisingly, the average prediction looks most similar to that made with the Hyperbolic Tangent since that model has a probability of 0.572. For low mortality values, the figure looks greener, with a red-green and a blue-green areas. In these areas, the Hyperbolic Tangent leads to bistability or a stable equilibrium, respectively, whereas the two other functions lead to oscillations and have a cumulative probability of 0.428. A representative example to illustrate the utility of these calculations occurs at m * = 1.08 and K * = 42.6 (point "+" in figure 2), where all functional responses predict completely different qualitative dynamics: (i) the Holling model predicts oscillations, (ii) the Hyperbolic Tangent predicts equilibrium coexistence, and (iii) the Ivlev predicts both of these dynamics with roughly equal probability (Table 2, upper half).

Identifying the sources of uncertainty

The probabilistic predictions presented in figure 2 present a form of predictive uncertainty, in the sense that predictions can change due to the uncertainty in the choice of a function and parametric uncertainty. As is commonly done in other studies, let us assume that we have chosen one function f among the possible ones. We will now explore the consequences of this choice, both in terms of parametric uncertainty and model uncertainty.

For the parametric uncertainty, consider two parameter sets θ f for function f drawn independently from their posterior distribution. We know from (4) the probability P (X|f ) that one parameter set predicts the dynamics X. Moreover, predicting dynamics X or not is a binary outcome. Therefore, the joint probability that one set of parameters predicts the dynamics X and the other does not follows a Binomial distribution and is

P (X|f ) [1 -P (X|f )
]. This probability is maximal (0.25) if there is an equal probability that f predicts dynamics X or not (P (X|f ) = 1 -P (X|f ) = 0.5), which corresponds to the highest uncertainty on predicting dynamics X. Conversely, the probability is null if there is no uncertainty in predicting dynamics X, i.e. P (X|f ) = 0 (f never predicts X) or P (X|f ) = 1 (f always predicts X). Thus, the total probability that these two parameter sets ever predict different dynamics is:

P param (f ) = X∈Q P (X|f ) [1 -P (X|f )] , (6)
where the sum is across possible qualitative dynamics. Similarly, we can define the probability that different dynamics are predicted by one parameter set θ f of function f and one parameter set of any of the other alternative functions:

P model (f ) = g∈F,g =f P (g) 1 -P (f ) X∈Q P (X|f ) [1 -P (X|g)] . (7)
The sum over alternative functions g is weighted in order to take into account the respective weight of each alternative to function f . The probabilities (6-7) quantify the probabilities of making different predictions knowing that we are looking at the parametric or the model uncertainty. This implies that we can define:

U param (f ) = P (f )P param (f ), (8)
which equals the parametric uncertainty that is proportional to the probability that f is the best function to describe new data, and:

U model (f ) = 1 -P (f ) |F | -1 P model (f ), (9)
which equals the model uncertainty that is proportional to the probability that another function than f is the best function to describe new data. Note that this second probability is averaged over the |F |-1 functions other than f , to give an equal weight to parametric and model uncertainty if all the alternative functions are equally plausible, independent of their number. In other words, we are comparing f against one of its alternatives, not f against all its alternatives, because in the latter case the approach would always indicate a high model uncertainty if the number of alternative functions is sufficiently high. If one chooses the Holling function in our example, the parameter uncertainty as defined in ( 8) is small in comparison to model uncertainty defined by ( 9), as the Holling function is very unlikely to be the best function (probability of 0.079) in comparison to the two alternative functions (figure 3a-b). According to this analysis, the total uncertainty associated to the choice of function f is:

U tot (f ) = U param (f ) + U model (f ), (10)
with the example for the Holling function shown in figure 3c.

To have a global overview, one can look at the average uncertainty across all the alternative functions that are considered (figure 3d):

Ūtot = f ∈F P (f )U tot (f ). (11)
To go deeper, we propose to quantify the respective importance of parametric and model uncertainty in the predictions made with a function f with:

S(f ) = 2    U model (f ) U tot (f ) -0.5    . (12)
Note that this index can be computed only if there is uncertainty in predictions (U tot (f ) > 0). It is positive if model uncertainty is greater than parametric uncertainty and negative otherwise. In addition, S(f ) = 1 if the function is not a suitable candidate (P (f ) = 0) or its parametric uncertainty does not affect the predicted dynamics, which are affected (in any amount) by model uncertainty. Conversely, S(f ) = -1 if predictions are affected by parametric uncertainty, and there is no model uncertainty i.e. P (f ) = 1. Finally, S(f ) = 0 can occur if (i) all alternative functions have an equal likelihood and they all lead to the same uncertain predictions (i.e. P (X|f ) = P (X|g) for all dynamics X and alternative functions g); or more generally (ii) the function f is more likely than others but others lead to predictions that are different enough to get

U param (f ) = U model (f ).
Based on the index ( 12) computed for each function, one can look at its average value across the alternative functions to get a global overview of the source of uncertainty in the model predictions:

S = f ∈F P (f )S(f ). ( 13
)
Figure 4a-c shows the index for the three alternative functional responses and our first data set. If Holling is chosen, model uncertainty is always higher than the parametric uncertainty, as alternative equations are more plausible and lead to different predictions. Conversely, if the Hyperbolic Tangent is chosen, parametric uncertainty is higher on average than model uncertainty (S(f (t) ) < 0 over 66.6 % of the parameter space).

This equation has the highest likelihood, which explains why its own parametric uncertainty appears to be more important. The Ivlev represents an intermediate case where model uncertainty is higher on average than parametric uncertainty (S(f (I) ) > 0 over 69.0 % of the parameter space). On average, choosing one equation creates a higher predictive uncertainty due to model uncertainty than parametric uncertainty in 66.3 % of the parameter space. It is worth noting that parametric uncertainty is higher for high mortality rates, which are close to the predator maximum growth rate that itself is proportional to 1/h • . Thus, this parameter value has a strong impact on the predicted dynamics, explaining why parametric uncertainty is higher. As a representative example, we show detailed calculations of these probabilities for (m * = 1.08, K * = 42.6) in table 2.

(a) Source of uncertainty with Holling

predator mortality (m) prey carrying capacity (K) 5 Overview of the method's possible outcomes and use

The underlying idea behind our approach is that the cost of choosing one model to make a prediction increases if (i) there is a high probability that another model also fits available data well and (ii) the two models make different predictions. Thus, this approach extends earlier studies on structural sensitivity [8, 9, among others] by considering the fact that, even if alternative models predict different dynamics, one model may outperform others in fitting available data. Again with the Rosenzweig-MacArthur model, we illustrate this possibility in figure 5 by using data on the ingestion rate (i.e. a functional response) of copepods (Calanus pacificus) feeding on diatoms (centric sp.) to parameterize the alternative functional responses [36, figure 4]. One functional response -the Hyperbolic Tangent-fits the data better than others, with a high probability (0.853) to be the best description of new data (table 3, upper rows). Thus, the predictions made with the Hyperbolic Tangent functional response strongly drive the average predictions across alternative functional responses (figure 5b).

Also, the areas of highest predictive uncertainty in figure 5c are due to the parametric uncertainty of the dominant function, as indicated in blue in figure 5d. Conversely, predictive uncertainty is low in areas where it is mostly due to model uncertainty (in red in figure 5d). This example illustrates the idea that, if a change in equations leads to different predictions, it mostly matters whether or not alternative equations are actually likely to be chosen.

As a third and last example, we fit the functional responses to the data on the ingestion rate of starved copepods Calanus pacificus feeding on diatoms (Thalassiosira fluviatilis) [36, figure 2]. Here, in contrast to the earlier examples, all three alternative functional responses have roughly the same likelihood (table 3, lower rows). As a result, the predictive uncertainty of the population model is mostly due to the model uncertainty (figure 6). This last example shares similarities to the one we used to introduce our approach. There, two alternative equations had a high likelihood in comparison to the third. Thus, one may imagine removing the most unlikely function from the analysis, and to keep only the equations giving an almost equally good fit. Doing so will end up either in one of two cases. First, something like figure 6 where some functions remain equally likely; or second, something like figure 5 where one of the alternative functions has a high likelihood in comparison to all others, and we can imagine only considering this best-fit function. In this last case, one may remove all functions except the one giving the best fit. By doing so, our approach simplifies in a parameter sensitivity analysis, as the scientist decides that the uncertainty in model construction -defined for the set of alternative equations that we choose to consider-can be neglected. Of course, deciding to neglect the uncertainty in model construction and choosing the best function (according to available data) can only be done after considering alternative models to fit to data, and having determined that one seems better than others.

Discussion

Here we present probabilistic predictions of system dynamics in the classical Rosenzweig-MacArthur model that take into account both model and parametric uncertainty. Overall, we show that uncertainty in model formulation regularly leads to a larger predictive uncertainty than does usual parametric uncertainty. Moving forward, it is worth noting that uncertainty likely also involves other parameter values and processes. For given predator mortality and prey carrying capacity, the resulting prediction uncertainty can be estimated directly from figure 3. Conversely, a known uncertainty based on data on the time-scale parameter ε or the intrinsic dynamics of the prey (here specified as logistic growth) can be taken into account by also sampling their posterior distributions simultaneously with those of the functional response and its parameters. Knowing the respective contribution of different biological processes to the resulting predictive uncertainty would be a major step forward. To go this way, sampling a higher-dimensional parameter space is feasible with existing algorithms [START_REF] Mcelreath | Statistical rethinking: a Bayesian course with examples in R and Stan[END_REF],

but performing each model analysis might be computationally prohibitive. Nevertheless, the cost of conducting such an analysis may ultimately be lower than the cost of making the wrong biological/environmental decision due to an unreasonable faith in model predictions.

Decreasing predictive uncertainty can be achieved by decreasing the uncertainty in either model construction, model parameterization, or both. To do this, the naive idea of improving the experimental data collection might help to greatly decrease the model uncertainty. In our three examples, we are studying a process that is assumed to be a saturating function of prey abundance. However, if this saturation is not present in the collected data (figure 6a), it necessarily becomes harder to identify the best model. Indeed, the model parameter defining the plateau is less constrained by data (in comparison to figures 1d-f,5a), increasing the parametric uncertainty of each alternative equations. Thus, alternative equations tend to have highly overlapping confidence intervals, subsequently increasing the uncertainty in model selection. Therefore, our study highlights the importance of designing experiments in a way that maximises the constraints on parameter values of alternative models to decrease the uncertainty in selecting the "best" model.

Finding the "best" model in a given situation may imply arguments beyond the simple fit of models to data, such as the fact that one of the alternative equations is a well-established model in the literature that is based on valuable theoretical arguments (e.g. underlying mechanisms). Though we did not do so here, such non-quantitative arguments can be taken into account in our quantitative analysis thanks to the concept of prior probabilities. Prior probabilities are used in Bayesian statistics to give an a priori weight to some parameter values, or here alternative models, before estimating their likelihood from data [START_REF] Mcelreath | Statistical rethinking: a Bayesian course with examples in R and Stan[END_REF]. With prior probabilities, one could thus give more weight to the well-established Holling functional response instead of the phenomenological Hyperbolic Tangent. As a consequence, this function might become the best model in the example in the table 3 (upper rows), where all functions are equally likely based solely on the data. Conversely, Holling functional response might remain the worst model in the example in table 3 (lower rows), where the Hyperbolic Tangent is about 400-fold more likely according to data.

Our framework can also be easily extend to consider quantitative as well as qualitative differences in predictions (e.g. equilibrium vs. oscillations). Quantitative differences are important for topics such as resource management, disease and pest control, or species invasions. Here, we combined Bayesian statistics with a qualitative analysis of model predictions (bifurcation analysis). Depending on model complexity and the question at stake, one can combine Bayesian statistics with quantitative predictions (e.g. numerical simulations) or include quantitative aspects in the bifurcation analysis. Quantitative aspects of structural sensitivity have already been considered in earlier studies [START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF][START_REF] Aldebert | Uncertainty in predictive ecology: consequence of choices in model construction[END_REF]. However, including quantitative aspects in the bifurcation analysis might become computationally prohibitive. Thus, getting the quantitative predictions of the model of interest is the challenge in extending our framework, which is well-suited to consider uncertainty in both quantitative and qualitative predictions.

The predictions made with the predator-prey model we used as examples could not be compared to data on the temporal dynamics of the studied system. Indeed, we extracted functional-response data from studies on grazing/ingestion rates of a few predator individuals -functional response sensu stricto-but these studies did not follow the temporal dynamics of a prey-predator system over many generations (i.e. the scale of the population model). Conversely, some experiments on population dynamics are not combined with functional response experiment, and functional-response parameters are optimised so that the predicted system dynamics fit the temporal data [e.g. [START_REF] Bw Kooi | The transient behaviour of food chains in chemostat[END_REF][START_REF] Curtsdotter | Ecosystem function in predator-prey food webs-confronting dynamic models with empirical data[END_REF]. However, if one has access to both types of data, the additional information on temporal dynamics can be used to constrain the probabilistic analysis. One way to do so is to remove candidate models that never predict the observed qualitative dynamics. Another complementary way is to perform the parameter estimation with constraints coming from fitting both the functional-response data and data on the temporal dynamics. This approach might be the best way to solve the issue of structural sensitivity. However, it is worth noting that this is a truly idealistic case. Indeed, for organisms with a long lifetime (months, years), collecting temporal data on their population dynamics would require a long-term monitoring (at least many years). Therefore, making probabilistic predictions prior to such an experiment would still be of interest, given the time needed to acquire data.

Though a similar idea has been used to improve parameter inference [START_REF] Ac Babtie | Topological sensitivity analysis for systems biology[END_REF], as far as we know this is the first time the full approach proposed here has been used. An intriguing application of our probabilistic approach would be to conduct structural-sensitivity analysis for additional data sets corresponding to other prey-predator species. This might allow the classification of organisms with population dynamics that are inherently more or less predictable, either because of parameter values or the uncertainty around their functional-response data. This approach can also be extended to different models, for example to test the recent hypothesis that mass-balanced prey-predator models with maintenance are less structurally sensitive [START_REF] Aldebert | Is structural sensitivity a problem of oversimplified biological models? insights from nested dynamic energy budget models[END_REF]. More generally, our proposed approach may benefit from further cross-fertilisation with the approach of partially-specified models [START_REF] Mw Adamson | When can we trust our model predictions? unearthing structural sensitivity in biological systems[END_REF][START_REF] Mw Adamson | Bifurcation analysis of models with uncertain function specification: how should we proceed[END_REF][START_REF] Mw Adamson | Defining and detecting structural sensitivity in biological models: developing a new framework[END_REF]. That approach provides generic results but cannot currently make probabilistic predictions of the sort we make here. Another way to include uncertainty would be to use Stochastic Differential Equations where the uncertain process is randomly drawn to include all data variability into model predictions.

To conclude, we have conducted a proof-of-concept study outlining a novel approach that considers both parametric uncertainty and uncertainty in biological model construction while presenting model prediction. We achieved this by bringing elements of Bayesian statistics into the analysis of deterministic dynamical systems.

Here, the prey-predator system considered is small enough to get analytical results on bifurcations, as this provides a comprehensive overview of qualitative model predictions. More complex models can also be studied by the proposed approach, by adapting the automatic model analysis to model complexity. For instance, a full overview of qualitative model predictions can be obtained if a numerical bifurcation analysis can be performed, or a sample of numerical simulations of different scenarios can be used for models as large as global ecosystem models. Some of these models are known to make different large scale predictions such as the global dominance of phytoplankton groups or the ocean primary production [START_REF] Tr Anderson | Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model[END_REF][START_REF] Tr Anderson | Sensitivity of secondary production and export flux to choice of trophic transfer formulation in marine ecosystem models[END_REF][START_REF] Vallina | Maximal feeding with active prey-switching: a kill-the-winner functional response and its effect on global diversity and biogeography[END_REF]. Those results were based on comparisons between alternative models, and the incorporation of our approach would allow to make probabilistic predictions based on different plausible biological models. Indeed, these small changes in model construction affect predictions at the ocean scale, but also the coexistence of alternative stable states and the predicted resilience at the scale of a few populations in interactions. Therefore, extending our approach to complex operational models, together with theoretical analyses to rank the processes and species according to the level of predictive uncertainty they produced at the ecosystem level, would be critical advances toward a better knowledge of the uncertainty and forecast horizon [as defined in 2] of model predictions in environmental sciences.

of which are managed by the Royal Society Te Apārangi (RDF-13-UOC-003 and 16-UOC-008).

Figure

  Figure2a-c presents model predictions based on each functional response and including parametric uncertainty. The effect of parametric uncertainty will be introduced in the subsection and is shown by the blurred colours in the figure. For the moment, the reader can have an overview of model predictions for given parameter values θ f by looking at the areas with sharp colours. In this model, there is a trivial extinction equilibrium E (0) = (0, 0) without any species, which can be reached only if the prey is initially absent. If both species are initially present, different asymptotic dynamics are predicted depending on the parameters and form of the functional response. An extinction equilibrium E (1) = (K, 0) without the predator always exists, but it is only stable in the white area in the figure. In other areas, prey carrying capacity is high enough to sustain the predator population and a coexistence equilibrium E (2) = (x(2) , y(2) ) exists. It starts to exist for higher

Figure 3 .

 3 Figure 3. Uncertainty in predictions made with Rosenzweig-MacArthur model. Details of the calculations if one chooses to use Holling functional response: uncertainty due to functional response's parameter values (a), to the choice of the functional response (b), and their sum (c). (d) Total uncertainty (like in (c)) averaged over the three alternative functional responses, weighted by their respective likelihood. Calculations at point "+" are detailed in table 2.

Figure 4 .

 4 Figure 4. Source of uncertainty in predictions made with Rosenzweig-MacArthur model. Relative importance of parametric (negative value, blue) and model (positive value, red) uncertainty in the resulting total predictive uncertainty (grey area: total prediction uncertainty lower than 0.01). (a,b,c) Source of uncertainty if one of the alternative functional responses is chosen. (d) Average over the three alternative functions, weighted by their respective likelihood. Calculations at point "+" are detailed in table 2.

Figure 5 .Figure 6 .

 56 Figure 5. Example where one equation is almost certainly the best one among candidates. Overview of the analysis with data on copepods (Calanus pacificus) feeding on diatoms (centric sp.). All panels are drawn similarly as in earlier figures: (a) functional responses are fitted to data by a HMCMC algorithm; (b) average qualitative predictions (probabilistic bifurcation diagram); (c) average total uncertainty in predictions; (d) source of uncertainty in predictions.

Table 1 .

 1 Comparison of the three functional responses fit to data on Gazella thomsoni feeding on grass swards, based on WAIC (Widely-Applicable-Information-Criterion). WAIC estimates are given plus or minus standard error.

	functional response	WAIC	difference with best WAIC weights P (f ) ranking
	Holling	-31.5 ± 9.2	4.0 ± 2.6	0.079	3
	Ivlev	-30.5 ± 9.1	1.0 ± 0.9	0.350	2
	Hyp. Tangent	-27.5 ± 9.1	0.0	0.572	1
	0.572 compared to 0.350 for Ivlev and 0.079 for Holling (exact values may change due to algorithm stochasticity,
	table

Table 2 .

 2 Example to detail probabilities related to uncertainty, for each functional response (knowing f ) and then averaged together (f -independent) based on WAIC weights from table1. This example corresponds to α = (K * = 42.6,m * = 1.08), the "+" point in figures 2-4.

	), showing how parametric uncertainty propagates

Calculations at point "+" are detailed in table 2.

Table 3 .

 3 Comparison of the three functional responses fit to data on copepods (Calanus pacificus) feeding on diatoms (centric sp. and Thalassiosira fluviatilis), based on WAIC (Widely-Applicable-Information-Criterion). WAIC estimates are given plus or minus standard error.

	prey	functional response	WAIC	difference with best WAIC weights ranking
		Holling	-39.4 ± 7.7	12.4 ± 2.9	0.002	3
	centric sp.	Ivlev	-48.2 ± 7.7	3.5 ± 1.2	0.145	2
		Hyp. Tangent	-51.7 ± 8.2	0.0	0.853	1
	Thalassiosira fluviatilis	Holling Ivlev Hyp. Tangent	-22.6 ± 4.8 -22.8 ± 4.5 -23.0 ± 3.5	0.4 ± 2.2 2.0 ± 1.6 0.0	0.299 0.332 0.370	3 2 1

Acknowledgements

We acknowledge David Nerini, Jean-Christophe Poggiale, Mathias Gauduchon, Claude Manté, Gerald Gregori and Andrew Letten for stimulating discussions. We also thank the anonymous referees for constructive comments that helped to improve the manuscript.

Funding Statement

CA received funding from European FEDER Fund under project 1166-39417. DBS received support from a Rutherford Discovery Fellowship and the Marsden Fund Council from New Zealand Government funding, both

Data Accessibility

Experimental data can be found in the references cited in the manuscript. The codes used to perform the analysis are available upon request to CA.

Author Contributions

Original idea by CA and DBS. CA performed the research and wrote a first version of the article. Both authors contributed to discussion, wrote the final manuscript and contributed to revisions.

Competing Interests

We have no competing interests.


Alexandre Garcia
email: <algarcia@enst.fr>.

Slim Essid

Chloé Clavel

Florence D'alché-Buc

Structured Output Learning with Abstention: Application to Accurate Opinion Prediction

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Up until recent years, opinion analysis in reviews has been commonly handled as a supervised polarity (positive vs. negative) classification problem. However, understanding the grounds on which an opinion is formed is of highest interest for decision makers. Aligned with this goal, the emerging field of aspect-based sentiment analysis (Pontiki Proceedings of the 35 th International Conference on Machine Learning, Stockholm, Sweden, PMLR 80, 2018. Copyright 2018 by the author(s). et al., 2016) has evolved towards a more involved machine learning task where opinions are considered to be structured objects-typically hierarchical structures linking polarities to aspects and relying on different units of analysis (i.e. sentence-level and review-level) as in [START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF]. While this problem has attracted a growing attention from the structured output prediction community, it has also raised an unprecedented challenge: the human interpretation of opinions expressed in the reviews is subjective and the opinion aspects and their related polarities are sometimes expressed in an ambiguous way and difficult to annotate [START_REF] Clavel | Sentiment analysis: from opinion mining to human-agent interaction[END_REF][START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF]. In this context, the prediction error should be flexible and should integrate this subjectivity so that, for example, mistakes on one aspect do not interfere with the prediction of polarity.

In order to address this issue, we propose a novel framework called Structured Output Learning with Abstention (SOLA) which allows for abstaining from predicting parts of the structure, so as to avoid providing erroneous insights about the object to be predicted, therefore increasing reliability. The new approach extends the principles of learning with abstention recently introduced for binary classification (Cortes et al., 2016) and generalizes surrogate least-square loss approaches to Structured Output Prediction recently studied in [START_REF] Brouard | Input output kernel regression: supervised and semi-supervised structured output prediction with operator-valued kernels[END_REF]Ciliberto et al., 2016;[START_REF] Osokin | On structured prediction theory with calibrated convex surrogate losses[END_REF]. The main novelty comes from the introduction of an asymmetric loss, based on embeddings of desired outputs and outputs predicted with abstention in the same space. Interestingly, similarly to the case of Output Kernel Regression [START_REF] Brouard | Input output kernel regression: supervised and semi-supervised structured output prediction with operator-valued kernels[END_REF] and appropriate inner productbased losses (Ciliberto et al., 2016), the approach relies on a simple surrogate formulation, namely a least-squares formulation followed by the resolution of a new pre-image problem. The paper is organized as follows. Section 2 introduces the problem to solve and the novel framework, SOLA. Section 3 provides statistical guarantees about the excess risk in the framework of Least Squares Surrogate Loss while section 4 is devoted to the pre-image developed for hierarchical output structures. Section 5 presents the numerical experiments and Section 6 draws a conclusion.

Structured Output Labeling with Abstention

Let X be the input sample space. We assume a target graph structure of interest,

G = (V = {⌫ 1 , . . . , ⌫ d }, E : V ⇥ V ! {0, 1}
) where V is the set of vertices and E is the edge relationship between vertices. A legal labeling or assignment of G is a d-dimensional binary vector, y 2 {0, 1} d , that also satisfies some properties induced by the graph structure, i.e. by E. We call Y the subset of {0, 1} d that contains all possible legal labelings of G. Given G, the goal of Structured Output Labeling is to learn a function f : X ! Y that predicts a legal labeling ŷ given some input x. Let us emphasize that x does not necessarily share the same structure G with the outputs objects. For instance, in Supervised Opinion Analysis, the inputs are reviews in natural language described by a sequence of feature vectors, each of them representing a sentence. Extending Supervised Classification with Abstention (Cortes et al., 2016), Structured Output Learning with Abstention aims at learning a pair of functions

(h, r) from X to Y H,R ⇢ {0, 1} d ⇥ {0, 1} d
composed of a predictor h that predicts the label of each component of the structure and an abstention function r that determines on which components of the structure G to abstain from predicting a label. If we note Y ? ⇢ {0, 1, a} d , the set of legal labelings with abstention where a denotes the abstention label, then the abstention-aware predictive model f h,r : X ! Y ? is defined from h and r as follows:

f h,r (x) T = [f h,r 1 (x), . . . , f h,r d (x)], f h,r i (x) = 1 h(x)i=1 1 r(x)i=1 + a1 r(x)i=0 . (1)
Now, assuming we have a random variable (X, Y ) taking its values in X ⇥ Y and distributed according to a probability distribution D. Learning the predictive model raises the issue of designing an appropriate abstention-aware loss function to define a learning problem as a risk minimization task. Given the relationship in Eq. ( 1), a risk on f h,r can be converted into a risk on the pair (h, r) using an abstentionaware loss a :

Y H,R ⇥ Y ! R + : R(h, r) = E x,y⇠D a (h(x), r(x), y). (2)
In this paper, we propose a family of abstention-aware losses that both generalizes the abstention-aware loss in the binary classification case (see (Cortes et al., 2016)) and extends the scope of hierarchical losses previously proposed by [START_REF] Cesa-Bianchi | Hierarchical classification: combining bayes with svm[END_REF] for Hierarchical Output Labeling tasks. An abstention-aware loss is required to deal asymmetrically with observed labels which are supposed to be complete and predicted labels which may be incomplete due to partial abstention. We thus propose the following general form for the a function:

a (h(x), r(x), y) = h wa (y), C a (h(x), r(x))i, (3)
relying on a bounded linear operator (a rectangular matrix) C : R p ! R q and two bounded feature maps:

a : Y H,R ! R p devoted to outputs with abstention and wa : Y ! R q , devoted to outputs without abstention. The three ingredients of the loss a must enable the loss to be non negative. This is the case for the following examples.

In Binary classification with abstention, we have Y = {0, 1} and the abstention-aware loss bin a is defined by : bin a (h(x), r(x), y) = 8 > < > :

1 if y 6 = h(x) and r(x) = 1 0 if y = h(x) and r(x) = 1 c if r(x) = 0
, where c 2 [0, 0.5] is the rejection cost; with r(x) = 0, in case of abstention and 1, otherwise. This can be written with the corresponding functions wa and a defined as:

wa (y) = ✓ y 1 y ◆ , C = ✓ 0 1 c 1 0 c ◆ , a (h(x), r(x)) = 0 @ h(x)r(x) (1 h(x))r(x) 1 r(x) 1 A .
H-loss (hierarchical loss): now we assume that the target structure G is a hierarchical binary tree. Then, E is now the set of directed edges, reflecting a parent relationship among nodes (each node except the root has one parent).

Regarding the labeling, we impose the following property : if an oriented pair (⌫ i , ⌫ j ) 2 E, then y i y j , meaning that a child node cannot be greater that his parent node. The Hloss [START_REF] Cesa-Bianchi | Hierarchical classification: combining bayes with svm[END_REF] which measures the length of the common path from the root to the leaves between these assignments is defined as follows:

H (h(x), y) = d X i=1 c i 1 h(x)i6 =yi 1 h(x) p(i) =y p(i) ,
where p(i) is the index of the parent of i according to the set of edges E, and c i is a set of positive constants nonincreasing on paths from the root to the leaves.

Such a loss can be rewritten under the form:

H (h(x), y) = h wa (y), C wa (h(x))i wa (z) = ✓ z Gz ◆ , C = ✓ 2diag(c) diag(c) diag(c) 0 ◆ ,
G is the adjacency matrix of the underlying binary tree structure and c the vector of weights defined above. The case of the Hamming loss can also be recovered by choosing:

wa (y) = ✓ y 1 y ◆ , a (h(x), r(x)) = ✓ 1 h(x) h(x) ◆ , C = I 2d ,
where I 2d is the 2d identity matrix.

Abstention-aware H-loss (Ha-loss): By mixing the H-loss and the abstention-aware binary classification loss, we get the novel Ha-loss which we define as follows:

Ha (h(x), r(x), y) = d X i=1 c Ai 1 {f h,r i =a,f h,r p(i) =y p(i) } | {z } abstention cost (4) + c Aci 1 {f h,r i 6 =yi,f h,r p(i) =a} | {z } abstention regret + c i 1 {f h,r i 6 =yi,f h,r p(i) =y p(i) ,a6 =f h,r i } | {z } misclassification cost
, where c Ai and c Aci can be chosen as constants or be function of the predictions. Thus, we have designed this loss so it is adapted to hierarchies where some nodes are known to be hard to predict whereas their children are easy to predict. In this case, the abstention choice can be used at a particular node to pay the cost c A for predicting its child. If this prediction is still a mistake, the price c Aci is additionally paid and acts as a regret cost penalizing the unnecessary abstention chosen at the parent. Acting on c A and c Ac provides a way to control the number of abstentions not only through the risk taken by predicting a given node but also its children. For sake of space, the dot product representation with wa and a of this loss is detailed in the supplementary material.

Empirical risk minimization for SOLA

The goal of SOLA is to learn a pair (h, r) from a i.i.d.

(training) sample drawn from a probability distribution D that minimizes the true risk:

R(h, r) = E x,y⇠D a (h(x), r(x), y), = E x,y⇠D h wa (y), C a (h(x), r(x))i.
We notice that this risk can be rewritten as an expected valued over the input variables only:

R(h, r) = E x hE y|x wa (y), C a (h(x), r(x))i.
This pleads for considering the following surrogate problem:

• Step 1: we define g ⇤ (x) = E y|x wa (y) = min g2(X !R q ) E x,y k wa (y) g(x)k 2 | {z } surrogate risk
. g ⇤ is then the minimizer of a square surrogate risk.

•

Step 2: we solve the following pre-image or decoding problem:

( ĥ(x), r(x)) = arg min

(y h ,yr)2Y H,R hg ⇤ (x), C a (y h , y r )i,
Solving directly the problem above raises some difficulties:

• In practice, as usual, we do not know the expected value of wa (y) conditioned on x: E y|x wa (y) needs to be estimated from the training sample {(x i , y i ), i = 1, . . . , n}. This simple regression problem is referred to as the learning step and will be solved in the next subsection.

• The complexity of the arg min problem will depend on some properties of a . We will refer to this problem as the pre-image and show how to solve it practically at a later stage.

These pitfalls, common to all structured output learning problems, can be overcome by substituting a surrogate loss to the target loss and proceeding in two steps:

1. Solve the surrogate penalized empirical problem (learning phase):

min g 1 n n X i=1 k wa (y i ) g(x i )k 2 + ⌦(g), (5)
where ⌦ is a penalty function and a positive parameter. Thus, get a minimizer ĝ which is an estimate of E y|x wa (y).

2. Solve the pre-image or decoding problem:

( ĥ(x), r(x)) = arg min

(h(x),r(x))2Y H,R hĝ(x), C a (h(x), r(x))i. (6)

Estimation of the conditional density E y|x wa (y) from training data

We choose to solve this problem in H ⇢ F(X , R q ), a vectorvalued Reproducing Kernel Hilbert Space associated to an operator-valued kernel K : X ⇥ X ! L(R q ). For the sake of simplicity, K is chosen as a decomposable operatorvalued kernel with identity:

K(x, x 0 ) = Ik(x, x 0 )
where k is a positive definite kernel on X and I is the q ⇥ q identity matrix. The penalty is chosen as ⌦(g) = kgk 2 H . This choice leads to the ridge regression problem:

arg min g2H n X i=1 kg(x i ) wa (y i )k 2 + kgk 2 H , (7)
that admits a unique and well known closed-form solution [START_REF] Micchelli | On learning vector-valued functions[END_REF][START_REF] Brouard | Input output kernel regression: supervised and semi-supervised structured output prediction with operator-valued kernels[END_REF].

As ĝ(x) is only needed at the prediction stage, within the preimage to solve, it is important to emphasize the dependency of ĝ(x) on the feature vectors wa (y i ):

ĝ(x) = n X i=1 ↵ i (x) wa (y i ), (8)
where ↵(x) is the following vector:

↵(x) = K x (K + I qn ) 1 , (9)
where

K x = [K(x, x 1 ), . . . , K(x, x n )]
. K is the qn ⇥ qn block matrix such that K i,j = K(x i , x j ) and I qn is the identity matrix of the same size. ↵ i (x) is the block i of ↵(x).

Learning guarantee for structured losses with abstention

In this section, we give some statistical guarantees when learning predictors in the framework previously described.

To this end, we build on recent results in the framework of Least Squares Loss Surrogate (Ciliberto et al., 2016) that are extended to abstention-aware prediction.

Theorem 1. Given the definition of a in (3), let us denote (h, r), the pair of predictor and reject functions associated to the estimate ĝ obtained by solving the learning problem stated in Eq. ( 7):

(h(x), r(x)) = arg min (y h ,yr)2Y H,R
hC a (y h , y r ), ĝ(x)i.

Its true risk with respect to a writes as:

R(h, r) = E x hC a (h(x), r(x)), E y|x wa (y)i.
The optimal predictor (h ⇤ , r ⇤ ) is defined as:

(h ⇤ (x), r ⇤ (x)) = arg min (y h ,yr)2Y H,R
hC a (y h , y r ), E y|x wa (y)i.

The excess risk of an abstention aware predictor (h, r): R(h, r) R(h ? , r ? ) is linked to the estimation error of the conditional density E y|x wa (y) by the following inequality:

R(h, r) R(h ? , r ? )  2c l q L(ĝ) L(E y|x wa (y)), (10) where L(g) = E x,y k wa (y) g(x)k 2 , and c l = kCk max y h ,yr2Y H,R k a (y h , y r )k R p .
The full proof is given in the Supplements. Close to the one in (Ciliberto et al., 2016), it is extended by taking the sup of the norm of a over Y H,R . Moreover when the problem ( 7) is solved by Kernel Ridge Regression, (Ciliberto et al., 2016) have shown the universal consistency and have obtained a generalization bound that still holds in our case since it relies on the result of Theorem 1 only. As a consequence the excess risk of predictors built in the SOLA framework is controlled by the risk suffered at the learning step for which we use off the shelf vector valued regressors with their own convergence guarantees.

In the following, we specifically study the pre-image problem in the SOLA framework for a class of output structures that we detail hereafter.

Pre-image for hierarchical structures with Abstention

In what follows we focus on a class of structured outputs that can be viewed as hierarchical objects for which we show how to solve the pre-image problems involved for a large class of losses.

Hierarchical output structures Definition

1. A HEX graph G = (V, E h , E e ) is a graph consisting of a set of nodes V = {v 1 , . . . , v n }, directed edges E h ⇢ V ⇥ V , and undirected edges E e ⇢ V ⇥ V , such that the subgraph G h = (V, E h ) is a directed acyclic graph (DAG) and the subgraph G e = (V, E e ) has no self loop. Definition 2. An assignment (state) y 2 {0, 1} d of labels V in a HEX graph G = (V, E h , E e
) is legal if for any pair of nodes labeled (y (i) , y (j) ) = (1, 1), (v i , v j ) / 2 E e and for any pair (y (i) , y (j) ) = (0, 1), (v i , v j ) / 2 E h .

Definition 3. The state space SG ✓ {0, 1} d of graph G is the set of all legal assignments of G.

Thus a HEX graph can be described by a pair of (1) a directed graph over a set of binary nodes indicating that any child can be labeled 1 only if its parent is also labeled 1 and

(2) an undirected graph of exclusions such that two nodes linked by an edge cannot be simultaneously labeled 1. Note that HEX graphs can represent any type of binary labeled graph since E h and E e can be empty sets. In previous works, they have been used to model some coarse to fine ontology through the hierarchy G h while incorporating some prior known labels exclusions encoded by G e [START_REF] Deng | Large-scale object classification using label relation graphs[END_REF][START_REF] Bentaieb | Topology Aware Fully Convolutional Networks for Histology Gland Segmentation[END_REF] While the output data we consider consists of HEX graph assignments , our predictions with abstention (h(x), r(x)) belong to another space Y H,R ✓ {0, 1} d ⇥{0, 1} d for which we do not restrict h(x) to belong to Y but rather allow for other choices detailed in the next section.

Efficient solution for the preimage problem

The complexity of the preimage problem is due to two aspects: i) the space in which we search the solution (Y H,R ) can be hard to explore; and ii) the a function can lead to high dimensional representations for which the minimization problem is harder.

The pre-image problem involves a minimization over a constrained set of binary variables. For a large class of abstention-aware predictors we propose a branch-and-bound formulation for which a nearly optimal initialization point can be obtained in a polynomial time. Following the line given by the form of our abstention aware predictor f h,r defined in Section 2, we consider losses involving binary interaction between the predict function h(x) and the reject function r(x), and suppose that there exists a rectangular matrix M such that a (h(x), r(

x)) = M 0 @ h(x) r(x) h(x) ⌦ r(x) 1 A
where ⌦ is the Kronecker product between vectors. Such a class takes as special cases the examples presented in Section 2. We state the following linearization theorem under binary interaction hypothesis: Theorem 2. Let l ha be an abstention-aware loss defined by its output mappings wa , a and the corresponding cost matrix C.

If the a mapping is a linear function of the binary interactions of h(x) and r(x) i.e. there exists a matrix

M such that 8(h(x), r(x)) 2 Y H,R a (h(x), r(x)) = M 0 @ h(x) r(x) h(x) ⌦ r(x)
1 A , then there exists a bounded linear operator A and a vector b such that 8 x 2 R p the pre-image problem:

( ĥ(x), r(x)) = arg min

(y h ,yr)2Y H,R h a (y h , y r ), x i,
has the same solutions as the linear program: ĥ(x), r(x) = arg min

(y h ,yr)2Y H,R [y T h y T r c T ]M T x s.t. A 0 @ y h y r c 1 A  b.
Where c is a d 2 dimensional vector constrained to be equal to y h ⌦ y r .

The proof is detailed in the supplementary material.

The problem above still involves a minimization over the structured binary set Y H,R . Such a set of solutions encodes some predefined constraints:

• Since the objects we intend to predict are HEX graph assignments, the vectors of the output space y 2 Y should satisfy the hierachical constraint : y i  y p(i) with p(i) the index of the parent of i according to the hierarchy. When predicting with abstention we relax this condition since we suppose that a descendant node can take the value y i = 1 if its parent was active y p(i) = 1 or if we abstained from predicting it r p(i) = 0. Such a condition is equivalent to the constraint

y i r p(i)  y p(i) r p(i) . (11)
• A second condition we used in practice is the restriction of the use of abstention for two consecutive nodes: structured abstention at a layer must be used in order to reveal a subsequent prediction which is known to be easy. Such a condition can be encoded through the inequality:

r i + r p(i)  1. (12)
In our experiments, the structured space Y H,R has been chosen as the set of binary vectors (h(x), r(x)) 2 Y H,R that respect the two above conditions. These choices are motivated by our application but note that any subset of {0, 1} d ⇥ {0, 1} d can be built in a similar way by adding some inequality constraints:

A Y H,R 0 @ h(x) r(x) h(x) ⌦ r(x) 1 A  b Y H,R
. Consequently, the Y H,R constraints can be added to the previous minimization problem to build the canonical form:

( ĥ(x), r(x)) = arg min

(y h ,yr) [y T h y T r c T ]M T x s.t. A canonical 0 @ y h y r c 1 A  b canonical , (y h , y r ) 2 {0, 1} d ⇥ {0, 1} d ,
where

A canonical = ✓ A A Y H,R ◆ and b canonical = ✓ b b Y H,R ◆ .
The complexity of the problem above is linked to some properties of the A canonical operator. [START_REF] Goh | Structured Prediction by Conditional Risk Minimization[END_REF] have shown that in the case of the minimization of the H-loss with hierarchical constraints, the linear operator A canonical satisfies the property of total unimodularity [START_REF] Schrijver | Theory of linear and integer programming[END_REF] which is a sufficient condition for the problem above to have the same solutions as its continuous relaxation leading to a polynomial time algorithm. In the more general case of the Ha-loss, solving such an integer program is NP-hard and the optimal solution can be obtained using a branch-and-bound algorithm. When implementing this type of approach, the choice of the initialization point can strongly influence the convergence time. As in practical applications, we expect the number of abstentions to remain low, such a point can be chosen as the solution of the original prediction problem without abstention [START_REF] Goh | Structured Prediction by Conditional Risk Minimization[END_REF]. Moreover since the abstention mechanism should modify only a small subset of the predictions, we expect this solution to be close to the abstention aware one.

Numerical Experiments

We study three subtasks of opinion mining, namely sentencebased aspect prediction, sentence-based joint prediction of aspects and polarities (possibly with abstention) and full review-based star rating. We show that these tasks can be linked using a hierarchical graph similar to the probabilistic model of [START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF] and exploit the abstention mechanism to build a robust pipeline: based on the opinion labels available at the sentence-level, we build a two-stage predictor that first predicts the aspects and polarities at the sentence level, before deducing the corresponding review-level values.

Parameterization of the Ha-loss

In all our experiments, we rely on the expression of the Haloss presented in 4. The linear programming formulation of the pre-image problem used in the branch-and-bound solver is derived in the supplementary material and involves a decomposition similar to the one described in Section 2 for the H-loss. Implementing the Ha-loss requires choosing the weights c i , c Ai and c Aci . We first fix the c i weights in the following way :

c 0 = 1 c i = c p(i) |siblings(i)| 8i 2 {1, . . . , d}.
Here, 0 is assumed to be the index of the root node. This weighting scheme has been commonly used in previous studies [START_REF] Rousu | Kernel-based learning of hierarchical multilabel classification models[END_REF][START_REF] Bi | Hierarchical multilabel classification with minimum bayes risk[END_REF] and is related to the minimization of the Hamming Loss on a vectorized representation of the graph assignment. As far as the abstention weights c Ai and c Aci are concerned, making an exhaustive analysis of all the possible choices is impossible due to the number of parameters involved. Therefore, our experiments focus on weighting schemes built in the following way:

c Ai = K A c i c Aci = K Ac c i
The effect of the choices of K A and K Ac will be illustrated below on the opinion prediction task. We also ran a set of experiments on a hierarchical classification task of MRI images from the IMAGECLEF2007 dataset reusing the setting of [START_REF] Dimitrovski | Hierchical annotation of medical images[END_REF] where we show the results obtained for different c i weighting schemes. The settings and the results have been placed in the supplementary material.

Learning with Abstention for aspect-based opinion mining

We test our model on the problem of aspect-based opinion mining on a subset of the TripAdvisor dataset released in [START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF]. It consists of 369 hotel reviews for a total of 4856 sentences with predefined train and test sets. In addition to the review-level star ratings, the authors gathered the opinion annotations at the sentence-level for a set of 11 predefined aspects and their corresponding polarity.

Similarly to them, we discard the "NOT RELATED" aspect and consider the remaining 10 aspects with the 3 different polarities (positive, negative or neutral) for each. We propose a graphical representation of the opinion structure at the sentence level (see Fig. 1). Objects in the output space y 2 Y consist of trees of depth 3 where the first node is the root, the second layer is made of aspect labels and the third one is the polarities corresponding to each aspect. The corresponding assignments are encoded by a binary matrix y 2 Y where y is the concatenation of the vectors indicating the presence of each aspect (depth 2) and the ones indicating the polarity.

An example of y encoding is displayed in Fig. 1. Based on the recent results of [START_REF] Conneau | Supervised learning of universal sentence representations from natural language inference data[END_REF], we focus on the InferSent representation to encode our inputs. This dense sentence embedding corresponds to the inner representation of a deep neural network trained on a natural language inference task and has been shown to give competitive results in other natural language processing tasks.

We test our model on 3 different subtasks. In Exp1, we first apply our model (H Regression InferSent) to the task of opinion aspect prediction and compare it against two baselines and the original results of [START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF].

In Exp2, we test our method and baselines on the problem of joint aspect and polarity prediction in order to assess the ability of the hierarchical predictor to take advantage of the output structure. On this task we additionally illustrate the behavior of abstention when varying the constants K A and K Ac . In Exp3, we illustrate the use abstention as a mean to build a robust pipeline on the task of star rating regression based on a sentence-level opinion predictor. results than [START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF]. Consequently, the next experiments will not consider them. Even though H Regression was trained in order to predict the whole structure, it obtains results similar to logistic regression and linear chain CRF.

Exp2. Joint polarity and aspect prediction with abstention. We take as output objects the assignments of the graph described (Fig. 1) and build an adapted abstention mechanism. Our intuition is that in some cases, the polarity might be easier to predict than the aspect to which it is linked. This can typically happen when some vocabulary linked to the current aspect has been unseen during the training or is implicit whereas the polarity vocabulary is correctly recognized. An example is the sentence " We had great views over the East River" where the aspect "Location" is implicit and where the "views" could mislead the predictor and result in a prediction of the aspect "Other". In such a case, [START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF] underline that the interannotator agreement is low. For this reason, we want that our classifier allows multiple candidates for aspect prediction while providing the polarity corresponding to them. We illustrate this behavior by running two sets of experiments in which we do not allow the predictor to abstain on the polarity.

In the first experiment, we want to analyze the influence of the parameterization of the Ha-loss. Following the parameterization of c Ai and c Aci previously proposed, we generated some predictions with varying values of K A 2 [0, 0.5] and K Ac 2 {0.25, 0.5, 0.75}. We displayed the Hamming loss between the true labels and the predictions as a function of the mean number of aspects on which the predictor abstained (Fig. 2) and handle two cases : modified : in the left figure, all nodes except the one on which we abstained were used to compute the Hamming loss. In the right one, all nodes except the aspect on which we abstained and their corresponding polarity were used to compute the Hamming loss. The H Strict results correspond to a predictor for which the original hierarchical constraint is forced: y (i)  y p(i) and the three other curves have been obtained with the generalized constraint hypothesis y (i) r p(i)  y p(i) r p(i) .

We additionally ran our model H Regression without abstention and our two baselines logistic regression for which we measured a similar Hamming loss of 0.03 (corresponding to 0 abstention on the left Figure 2). Concerning the micro-averaged F1 score, the H Regression retrieved a score of 0.54 being slightly above the logistic regression which scored 0.53 and the linear chain CRF with 0.52.

Two conclusions can be raised. Firstly, the value of K Ac and the choice of the hypothesis H Strict have little to no influence on the scores computed in the two cases previously described. Secondly, increasing the number of abstentions on aspects helps reducing the number of errors counted on the aspects nodes when the predictor abstains on less than 3 labels. After this point, the quality of the overall prediction decreases since the error rate on the remaining aspects selected for abstention is less than the one on the polarity labels Subsequently, we examine the Hamming loss on the polarity predictions situated after an aspect node to understand the influence of the c Ac coefficients and the relaxation of the H Strict hypothesis in Figure 3. The orange curve gives the best score when the mean number of abstentions is between 2 and 4 per sentence. The only difference with the H strict hypothesis is the ability to predict the polarity of an aspect candidate for abstention even if the predictor function does not select it. This behavior is made possible by the fact that Figure 3. Hamming loss computed on polarity nodes located after an aspect for which the predictor abstained our prediction does not respect the Y constraints but instead belong to the more flexible space Y H,R Finally we show how abstention can be used to build a robust pipeline for star-rating regression.

Exp3. Star rating regression at the review level based on sentence level predictions. In the last round of experiments, we show that abstention can be used as a way to build a robust intermediate representation for the task of opinion rating regression [START_REF] Wang | Latent aspect rating analysis without aspect keyword supervision[END_REF] which consists in predicting the overall average star rating given by each reviewer on a subset of six predefined aspects. The figure below illustrates the different elements involved in our problem. The procedure is split in two steps. Firstly, we learn a sentence-level opinion predictor that takes advantage of the available annotations. This step corresponds to the one studied in the previous experiment. Then a vector-valued regressor (star regressor in Figure 4) is built. It takes as input the component-wise average of the sentence level opinion representations, and intends to predict the star ratings at the review level. For each of the five overall aspects a separate Ridge Regressor is trained based on the true labels available. Once learned, the regressors take as input the prediction of the first step in a pipelined way Similarly to [START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF], we rescale the star ratings on a (-1,0,1) scale and report the macro-averaged mean average error on the test-set in

Conclusion

The novel framework, Structured Learning with Abstention, extends two families of approaches: learning with abstention and least-squares surrogate structured prediction. It is important to notice that beyond ridge regression, any vectorvalued regression model that writes as (8) is eligible. This is typically the case of Output Kernel tree-based methods [START_REF] Geurts | Kernelizing the output of tree-based methods[END_REF]. Also, SOLA has here been applied to opinion analysis but it could prove suitable for more complex structure-labeling problems. Concerning Opinion Analysis, we have shown that abstention can be used to build a robust representation for star rating in a pipeline framework. One extension of our work would consist in learning how to abstain by jointly predicting the aspects and polarity at the sentence and text level.

Figure 1 .

 1 Figure 1. Graphical representation of the opinion structure Exp1. Aspect prediction. In this first task, we aim at predicting the different aspects discussed in each sentence. This problem can be cast as a multilabel classification problem where the target is the first column of the output objects y for which we devise two baselines. The first relies on a logistic regression model (Logistic Regression InferSent) trained separately for each aspect. The second baseline (Linear chain Conditional Random Fields (CRF) (Sutton et al., 2012) InferSent) is inspired by the work of (Marcheggiani et al., 2014) who built a hierarchical CRF model based on a handcrafted sparse feature set including one-hot word encoding, POS tags and sentiment vocabulary. Since the

Figure 2 .

 2 Figure 2. Hamming loss as a function of the number of aspect labels where the predictor abstained itself.

Figure 4 .

 4 Figure 4. Star rating regression pipeline

  opinion representations are assumed to be known on the test set and fed in the text-level opinion regressors to find back the star ratings. The Hierarchical CRF line corresponds to the best results reported by[START_REF] Marcheggiani | Hierarchical multi-label conditional random fields for aspect-oriented opinion mining[END_REF] on the two tasks. H Regression is our model without abstention used as a predictor of the sentence-level representation in the pipeline shown in Fig 4.Finally for the H Regression with abstention, we used as a sentence-level representation : y a = h(x) (1 r(x)). Since the only non-zero components of (1 r(x)) correspond to aspects on which we abstained, subtracting them from the original prediction results in a reduction of the confidence of the regressor for these aspects and biasing the corresponding polarity predictions towards 0. H Regression strongly outperforms Hierarchical CRF on both tasks. We do not report the score for H Regression with abstention since it is dependent on the number of abstentions but show that it improves the results of the H Regression model on the text-level prediction task. The significance of the scores has been assessed with a Wilcoxon rank sum test (p-value 10 6).

Table 5

 5

	.2 below shows the results in terms of micro-averaged
	F1 (µ-F1) score obtained on the task of aspect prediction.
	The three methods using InferSent give significantly better
	method	µ-F1
	H Regression InferSent	0.59
	Logistic Regression InferSent	0.60
	Linear chain CRF InferSent	0.59
	Linear chain CRF sparse features Marcheggiani et al.	0.49
	Hierarchical CRF sparse features Marcheggiani et al.	0.49

Table 1 .

 1 Experimental results on the TripAdvisor dataset for the aspect prediction task.

Table 2 .

 2 Table 5.2 below under the column MAE text level. We additionally include the MAE error measured on polarity predictions at the sentence level counted when the underlying aspect predicted is a true positive. The first row is our oracle: the sentence-level Experimental result on the TripAdvisor dataset for the polarity prediction task

	method	M AE sentence level	M AE text level
	Oracle: regression with true sentence labels	0	0.38
	Hierarchical CRF	0.50	0.50
	H Regression	0.30	0.45
	H Regression with Abstention	C	0.43

Acknowledgements

This work has been funded by the french ministry of research and by the chair Machine Learning for Big Data of Télécom ParisTech.