P Caillou

J Gil-Quijano

X Zhou

Automated observation of multi-agent based simulations A statistical analysis approach

Keywords: Complex systems simulation, multi-agents systems, automated observation, automated characterization, clustering, value-test Studia Informatica Universalis

Multi-agent based simulations (MABS) have been successfully used to model complex systems in different areas. Nevertheless a pitfall of MABS is that their complexity increases with the number of agents and the number of different types of behavior considered in the model. For average and large systems, it is impossible to validate the trajectories of single agents in a simulation. The classical validation approaches, where only global indicators are evaluated, are too simplistic to give enough confidence in the simulation. It is then necessary to introduce intermediate levels of validation. In this paper we propose the use of data clustering and automated characterization of clusters in order to build, describe and follow the evolution of groups of agents in simulations. These tools provides the modeler with an intermediate point of view on the evolution of the model. Those tools are flexible enough to allow the modeler to define the groups level of abstraction (i.e. the distance between the groups level and the agents level) and the underlying hypotheses of groups formation. We give an online application on a simple NetLogo library model (Bank Reserves) and an offline log application on a more complex Economic Market Simulation.

Introduction

Multi-agent systems (MAS) are specially well suited to represent complex phenomena from the description of local agent behaviours. The simulation of complex systems using MAS is a cyclic process : the modeler introduces his/her knowledge into the model, runs simulations, discovers bugs, pitfalls or unwanted effects, corrects the model and eventually his/her knowledge, and the cycle restarts. The cycle is over when it is not possible to further improve the model because of technical or knowledge limitations.

Once the agents behaviors are defined, the cyclic modeling process usually focuses on fitting global simulation parameters and/or the initial states of agents in order to reproduce global behaviors observed in the modeled empirical phenomena. The global behaviors are usually reproduced and tested with global variables (for example, in the case of socio-spatial models, variables as populations density and growing rates, or evacuation time and number of death for panic simulations). Those global variables are evaluated both on simulation and on empirical data. The calibration of the simulation is achieved by finding the right set of values (or values intervals) for the agents and global simulation parameters which lead to minimize the difference between trajectories of global variables in simulation and in empirical data. This optimization-like approach is implemented by several existing frameworks (for example in GAMA [START_REF] Taillandier | Gama : a simulation platform that integrates geographical information data, agent-based modeling and multi-scale control[END_REF], see section 2 for an extended presentation of these tools).

Nevertheless, this traditional approach may be too simplistic in order to characterize the dynamics of complex systems. Indeed, in a complex system, different phenomena may simultaneously occur at different levels (at the agents and at the global levels, but also at intermediate levels) and influence each other [START_REF] Gil-Quijano | From biological to urban cells : lessons from three multilevel agent-based models[END_REF]. For instance, groups of agents (flocks of birds, social groups, coalitions, etc.) following similar trajectories of states may appear, evolve and disappear. To describe and evaluate the evolution of that type of groups, the observation of global variables is not enough. Moreover, because of the emergent properties of complex systems, those groups may be unexpected, and their presence may even be unnoticed because no global variable or any other adapted observation mechanism is provided in the simulator. The significance and even the existence of groups may then be hidden by the usually huge amount of information generated in a MAS simulations.

In this paper, we introduce the use of statistical-based tools to assist the modeler in the discovering, following the evolution and describing groups of agents. After an overview of the state of the art (section 2, in section 3, we present two complementary tools, data clustering used to discover and build groups and value test used to automatically describe those groups. In section 4, we present an observation model that uses those tools in order to produce automated analysis of the evolution of groups in MAS simulations illustrated on a NetLogo simple model. We present a more complex offline application in section 5 and conclude in section 6.

State of the art

Multi-agent based simulations have been used with a large number of economic, geographic or social applications. There are several available development frameworks for simulation, some of them userfriendly with specific coding language, such as NetLogo [LM], and with the possibility to interface with Java code parts (like GAMA [START_REF] Taillandier | Gama : a simulation platform that integrates geographical information data, agent-based modeling and multi-scale control[END_REF]). Others use only generic language (usually java or C#), such as MO-DULECO [Pha] or Repast [START_REF] Michael | Experiences creating three implementations of the repast agent modeling toolkit[END_REF][START_REF] Railsback | Agent-based simulation platforms : Review and development recommendations[END_REF]. However, none of these platforms integrate any module for automatic group analysis. Based on these platforms, some analyzer tools, such as LEIA [START_REF] Gaillard | A reverse engineering form for multi agent systems[END_REF], SimExplorer [LR] and [START_REF] Caillou | Automated multi-agent simulation generation and validation[END_REF] were developed to generate and analyze automatically simulations.

LEIA [START_REF] Gaillard | A reverse engineering form for multi agent systems[END_REF] is a parameter space browser for the IODA simulation framework [START_REF] Kubera | Interaction-oriented agent simulations : From theory to implementation[END_REF]. It allows the user to instantly make a visual comparison of numerous simulations by seeing all their results in parallel. It provides the user with a set of transformation and generation tools for model, and a set of tests to browse the simulations space. Scoring rules are applied to help the user in identifying interesting configurations (such as cyclic or regular behaviors).

SimExplorer/OpenMole [LR] is a software, which aims at providing a generic environment for programming and executing experimental designs on complex models. The goals are multiple : (1) to externalize the development of the model exploration, in order to make available some generic methods and tools which can be applied in most of the cases for any model to explore ; (2) to favor the reuse of available components, and therefore lower the investment for good quality model exploration applications ; (3) to facilitate a quality insurance approach for model exploration.

[Cai10] is a tool to automatically generate and run new simulations until the results obtained are statistically valid using a chi-square test. It can generate new simulations and perform statistical tests on the results, with an accuracy that increases gradually as the results are produced. This tool can be applied to any RePast-based simulation. It deduces variables and parameters used and asks the user to choose the configuration of interest. New simulations are generated, computed and analyzed until all the independence tests between parameters/variables are valid. Finally, the test results and their margins of error are presented to the users.

The aim of these tools is to study several simulations (the parameter space), to compare their result and analyze them. However, none of them aim at studying one complex simulation to describe it. To explore one simulation, the only existing tools are the integrated tools (such as the NetLogo graphs and logs), which are limited to global or userdefined clusters, and classic data mining on logs. We aim to combine the advantage of online and agent-oriented analysis of NetLogo with the flexibility and descriptive potential of Data Mining tools. Some work using data mining tool to identify groups and describe them had been realised with specific applications, for exemple in the SimBogota simulation ([GQPD07][EHT07]). In these simulations, social groups in Bogota city regions where identified by data mining, and the group results were perceived by the agent. The goal was however more a multi-scale simulation than a description of simulation dynamic.

Analysis tools

In this section, we present the two main tools that we use to automatically analyze MAS simulations. To discover the groups of agents we propose to use data clustering, and then value-test evaluations to describe them. These tools are associated in our analysis model (see section 4) in order to automatically describe the evolution of groups and help the modeler to understand what happens in complex simulations.

Finding the groups : clustering

The goal of clustering algorithms is to "find the structure" of a dataset. Most of the time, data represent objects or individuals that are described by a given number of variables or characteristics [START_REF] Lebart | Statistique exploratoire multidimeensionnelle : visualisation et inférence en fouilles de données[END_REF]. The dataset's structure is represented as a partition or a hierarchy of partitions. Every single object is assigned to a given group (cluster) in the partition, or to a several groups (clusters) when considering a hierarchy of partitions. An object is assigned to a given group g if it is more similar (the sense of similar varies with the algorithm) to the objects in g than to objects in other groups. The main hypothesis when clustering a dataset is that the structure exists and the goal is to make it evident. Similarity between objects usually depends on the distance between them (actually between the vector of variables representing them). One of the most used distance (for quantitative variables) is the Euclidean distance.

As the state of an agent is a vector that includes the instantaneous values of the set of variables that describes the agent's behavior, we can consider the set of states of all the agents in a simulation as a dataset and then to cluster it. In that way, groups of agents whose states similarity is maximal can get formed. To conduct the clustering, the "right" algorithm and distance measure have to be chosen. Indeed, different clustering algorithms can lead to different results as their underlying hypothesis and functionning diverge. In our work, it is the responsibility of the modeler to choose those parameters. By including in our model the Weka machine learning library 1 , we provide the modeler with a wide list of clustering algorithms and distances functions. In our experiments we used the X-Means algorithm [PM00], described in the following paragraphs.

One of the most known algorithm of clustering is the K-Means algorithm [START_REF] Stuart | Least squares quantization in pcm[END_REF] whose objective is to find the k prototypes (average characteristic vectors) that represents the best the data. In that algorithm k initial prototypes are defined (usually by random) and at each iteration every object is assigned to its nearest prototype. Every prototype is then updated to the average vector of the characteristics of the objects that were assigned to it. The process is repeated until there is no significant changes of prototypes between two succesive iterations or until a maximum number of iterations is reached. The main pitfall of K-Means is that the number k of clusters must be known. As the idea of clustering is to find the right and unknown structure describing the dataset, there is no reason to know k beforehand.

An improvement of K-Means is proposed with the X-Means algorithm [PM00]. In that algorithm the "right" number k of clusters is determined by successive K-Means executions. The first execution starts with a k min number of clusters (the minimum number of clusters, a parameter of the model), and at each iteration one of the clusters found in the previous iterations is divided into two new clusters. The cluster to be divided is the one whose internal similarity (the similarity between the objects inside the cluster) is the lowest. The process is repeated until k max number of clusters is reached. Then k maxk min + 1 partitions are produced. The chosen partition is the one that maximizes the internal similarity of clusters and maximizes the distance between prototypes.

Interpreting the groups : Value tests

Value test (VT) [START_REF] Morineau | Note sur la caractérisation statistique d'une classe et les valeurs-tests[END_REF][START_REF] Lebart | Statistique exploratoire multidimeensionnelle : visualisation et inférence en fouilles de données[END_REF] is an indicator that allows the automated interpretation of clusters. It determines the more significant factors (for continuous variables) and modalities (for categorical variables) in a given cluster in comparison with the global dataset. The VT compares the deviation between the variables/modalities in clusters and the variables/modalities in the overall dataset. The main hypothesis in the VT calculation is that the variables follow Gaussian distributions. In that condition, for a level of risk of 5% we can consider Figure 1 -Overview of the analysis model that a variable/modality is significant if its VT is greater than 2. The automated description of a group is given by its set of significant variables/modalities. For continous variables the mean value of the variable in the group completes the description.

We present here the calculation of V T for continuous variables, for categorical variables see [START_REF] Morineau | Note sur la caractérisation statistique d'une classe et les valeurs-tests[END_REF]. Given a dataset containing n elements and a cluster k on the dataset containing n k elements. Given a quantitative variable j, its average E(n j ) and its variance S 2 (n j ) in the overall dataset. Given also the average E(n j k ) of j in the cluster k, the V T for the variable j in the cluster k is computed as follows :

V T (n j k ) = (E(n j k ) -E(n j )) ( n-n k n-1 × S 2 (n j ) n k ) (1)

Analysis model

Model overview

Our goal is to describe, online or offline, what happens in a simulation at the cluster level. Our model can be described with several steps as illustrated in Fig. 1 : 1) Model Selection : what do we study ? 2) Data processing : what are the data ? 3) Clustering : can we find homogeneous groups ? 4) Cluster description : how can we describe them ? 5) Cluster evolution : how do they evolve ? 6) Simulation generation : is this reproducible ? In future work, we intend to use the most interesting (or user-selected) agent model (clusters) identified to generate new simulations with similar agents and thus test the clusters behavioral stability.

For a better understanding, we will describe each step with the application of our tool to an illustrative example.

Model selection

The first step is to choose the model to be studied. Our model can be applied both online (with NetLogo) or offline from logs (by simulating an online data stream).

We choose here the Bank Reserves model, provided with Netlogo, where financial agents either save or borrow money via loans (Fig. 2). This is a very simple model illustrating the effect of money creation via savings deposits/loans grants. There is only one Bank and P eople (consumers) agents. Each agent begins with a random amount of money in its wallet (between 0 and a parameter richT hreshold). When an agent has a positive wallet, it deposits its money in the bank, which increases its saving variable (and puts its wallet at 0). At each step, agents move randomly. When they meet someone they make a transaction, which is a simple transfer from one agent to the NetLogo provides some tools to observe an experiment either at an individual or at a global level. For example, on Fig. 2 some global variables are presented to give an overview of an experiment. The global amount of loans and money show an early increase, then a stabilization of the total money when the maximum amount of loans is reached. The Income distribution graph gives an overview of the repartition of wealth between three fixed groups (negative wealth, wealth higher than richT hreshold and the rest). Even if these informations are interesting, a more detailed understanding of the model behavior can not be reached with such global/local observation. For example, it is difficult to answer the following questions : who are the wealthy agents ? Do the rich stay rich ? This would even be truer for more complex models, for which variable interactions are much more difficult to deduce than with such a simple toy simulation.

Data processing

A data matrix is generated every n steps. A line in the matrix represents one agent's state. Raw data from simulations are not the only interesting data for cluster's generation and analysis. Several filters or aggregators can also be used to process the data stream. We use two different aggregators to complete the initial matrix : i) the moving average (for each variable, we add a new variable computed as the average of the last five steps values) ; ii) the initial values for each agent (for each variable, we define a new variable corresponding to the value of the variable for this agent at the starting point of the simulation). These initial values variables are not used in the clustering but used latter for the description of the clusters.

Clustering

Clustering is performed on the final data in order to generate homogeneous agent groups (for now, X-Means is used, but any other clustering algorithm from Weka can be easily selected instead). Clusters are visualized in NetLogo (with colors), and their extension and description are presented. For example, in Fig. 3, three clusters are identified in t=400.

Cluster's description

Once the clusters are identified, it is possible to get an easy-to-read description by using V T (see sec. 3.2). The description with V T (Fig. 4) makes it easy to interpret and describe them. Positively (respectively negatively) significant variables are presented in blue (resp. red) : their average is significantly higher (resp. lower) than the global average. For example, in t = 400, three clusters are identified. Cluster7, with 114 agents, is a "poor" cluster, whose agents have low wealth, savings, wallet and the corresponding moving average variables (M M savings and M M wealth), and a higher amount of loans. Some significant variables are (probably) clustering artifacts (Y Cor) or random effects (T OHeading), and will justify our stability analysis.

Similarly, cluster9 regroups the 66 "wealthy" people, with high wealth and savings and few loans. An interesting result is the significant T 0W allet variable, corresponding to the wallet value of agents at the beginning of the simulation. The wealthy people were significantly richer than the average at the beginning of the simulation.

At the end of the simulation, a description of all the clusters obtained at each time step gives a global overview of the simulation (Fig. 5, with a selection of some results in Table 1). In our experiment, it is always possible to identify a wealthy and a poor cluster, and sometimes (like in t = 400) a middle cluster. From their description, it is already possible to observe that the link between the wealth and the initial wealth (the T 0W allet) is not significant anymore after t = 400. It may be related to the fact (observed with NetLogo global observation) that bank has Figure 4 -Cluster Description at t=400 : listed variables are higher for the member of the cluster (blue positive value) of lower (red negative values) than in the global population. Black variables are not significantly different.

Figure 5 -Global Overview of clusters reached its loan limit (the total amount of money stops to increase at around t = 230).

However, it is difficult to compare clusters at different time steps with this overview since they are different both in intension and in extension. In a more complex model, cluster may have a completely different meaning at different steps.

Tableau 1 -Selection of the "rich" clusters and some variables results from the global overview of a simulation (Fig 5

Clusters evolution

In order to describe the clusters' evolution, we consider two alternative hypothesis : either the extension in every cluster is considered as stable (we keep exactly the same agent population in the cluster), or the intension of every cluster is fixed (we keep the same definition of the cluster).

Evolution by extension/population

Once an interesting cluster is identified (for example the wealthy agents of t = 400, cluster9), it is interesting to follow its evolution. To do it, the first way is to fix the extension (population) of the cluster. Fig. 6 describes the evolution of the cluster9 with fixed extension after t = 400 (see Table 2 for a selection of the most interesting variables). Initial parameters values are stable since the population does not change. Other variables may change (except for example Id since this variable is constant for every agent). This view clearly shows that all wealth-related differences with the other agents decrease : all the (ab-Figure 6 -Cluster evolution of Cluster9 (rich peoples of T=400) by extension (fixed population) solute) V T values for wealth, saving and loans decrease. This mean that, in average, the wealthy people of t = 400 are becoming less and less wealthy. They are still significantly wealthier than the average at t = 1200, but their loans are no more significantly lower in comparison with t = 1000. We can also check on this evolution that the clustering artifacts (like the Y Cor variable) do not stay significant.

Evolution by definition

The second way is to fix the clusters intension (definition). Fig. 7 represents the description of the clusters identified at each step with the intension function of t = 400 (see Table 3 for a selection of interesting variables). All the variables considered in clustering are by definition roughly similar, since the intension of clusters is the same. However, the other variables may evolve (in our example, the initial parameters of the agents). Cluster9, for example, regroup the wealthy agents at each step, but the number and the initial properties of these agents evolve. It is interesting to see that the number of wealthy agent stays approximatively constant (66, 71, 56, 58, 65). But the evolution of the initial parameters confirms the observation made with the global overview : the initial wealth of the agent (T OW allet) is not significant anymore after t = 600.

Cluster evaluation

The evaluation of the cluster is done by combining the size of the cluster and its 'descriptivness' measured by the number of significant V T . The objective is to identify clusters that are both big enough to be generic and descriptive enough to bring some interesting informations. The score of a cluster c at time step t is calculated as the product of the number of significant variables (whose '|V T | is greater than 2) : score(c, t) = |V S c,t | × n c,t where V S c,t is the set of significant variables of c and n c,t is the number of agents in c at the time t.

For example, in Fig. 6, the score of the cluster is initially relatively high (924) because the population is important and the number of significant V T is high. But since the cluster is not stable (the rich people don't stay rich), the number of significant variables decreases. The score decreases, reflecting the fact that the cluster becomes difficult to interpret (except by : "the ones that were poor at t=400").

Experiments

We tested our analysis tool on several simulation models, both online with NetLogo and offline with data logs. Since we have already described a NetLogo application in the previous section to illustrate the description tool, we will describe here an offline analysis application.

Model description

To illustrate how our model deals with more complex simulations, we chose a model following the KIDS approach [START_REF] Edmonds | From kiss to kidsan 'anti-simplistic' modelling approach[END_REF] : the number of parameters and observed variables is kept high to be more descriptive and realistic rather than synthetic. The Rungis Wholesale Market simulation [CCB09a][CCB09b] was developed with the BitBang Franmework [START_REF] Baptista | Bitbang : A model and framework for complexity research[END_REF] and reproduces a Fruit and Vegetable wholesale market. One type of seller agent and 4 types of Buyer agents are considered in the simulation (with many variable parameters, 20 variables in average by agent type). The four type of buyer agents are : Restorators seeking efficiency, TimeFree seeking good opportunities, Barbes seeking low-quality and low price products and Neuilly buyers seeking high quality high prices products. The csv logs record one line for each couple day/agent, with 33 output variables (see below). We analyze here only the 60 Buyer agents during the first 10 days of the simulation.

The main observed variables of this model are the transaction time, the number of sellers per buyer, the quality and quantity of the products and the prices. There are four types of prices, the producer price (price paid by the sellers to the producers), the transaction price(price paid by the buyers to the sellers), the standard price (price given at the beginning of the transaction) and the final price (price paid by the consumers). In general producerprice < transactionprice < standardprice < f inalprice.

Clusters description

Two clusters are identified at each step. They are easy to describe since many variables are significant (see Fig 8 where the high proportion of red or blue numbers illustrates the high number of significant variables, and a selection and description of interesting variables for the cluster 3 in Table 4). For example, at t=0 (and t=1) we identify the "expansive" cluster that is composed by the buyers of high-quality expansive products, and the "cheaps" cluster that is composed by the agents buying low-quality cheap products. Even if this could be anticipated from the buyers definition, one first quick analysis helps to validate the model behaviors : the "expansive" products are fresher (product Age), with higher prices, and the Buyer and Seller profit as well as profit rate higher than for the "cheaps". Other significant variables give some new interesting informations : the buying time (M oyHour) is not significantly different, but the expansive buy significantly more products (Sumof Qty) to more Sellers (N bSeller) than the cheaps agents. These informations were not trivial to identify unless you knew you want to find it before the analysis.

Cluster's evolution

Following the identified clusters brings new interesting insights. For example, in Fig. 9 (and some selected variables in Table 5), we follow the "expansive" cluster identified at t=1. First, it is possible to check that the clusters has a stable behavior. The agents buying expansive and high-quality products still do the same in the following days. It is however possible to describe some cluster evolutions : the Moy-Hour variable, which measures the Average transaction time shows a constant decrease. "expansive" agents buy progressively sooner than the "cheaps" agents. Also, the SumQty variables, which was significant when the cluster was created, becomes quickly not-significant : the Avg Seller margin rate 5,29 expansive buyers do not buy more product thant the orthers, it was just an artifact at the cluster creation.

Conclusions

The framework for the observation of MAS simulation that we present here, provides the modeller with generic tools that allow him/her to get a synthetic descriptive view of MABS. Presently, it can be used to understand the dynamics of simulations and to ease their validation. In future work we intend to develop a mechanism that will allow the modeller to use the definitions of interesting clusters in new simulations as generic "agent models". Indeed, cluster definition and population can be used to define a distribution function to generate agents profiles. It is easy to retrieve the average values and the standard deviation of every variable by using simple statistical analysis tools. Agents generated using the distribution functions could be reintroduced in simulations, clusters can be rebuilt and the global simulation variables can be compared with their previous values. In that way, the clusters stability and their "expressiveness" can be measured over different simulations.

In order to allow the analysis of a wide number of different type of simulations we are currently adapting our framework both to consider qualitative and network variables and facilitate large simulations analysis. The latter will be done by integrating our framework to the Open-Mole engine. That engine provides, among other facilities, the easy use of cluster and grid computing for simulations.

Figure 2 -

 2 Figure 2 -Netlogo Bank Reserves Model

Figure 3 -

 3 Figure 3 -The cluster list panel gives the identity number of the agents included in each cluster.

Figure 7 -

 7 Figure 7 -Cluster evolution of Cluster9 (rich peoples defined at T=400) by intension (fixed definition)

Figure 8 -

 8 Figure 8 -Global Overview of clusters after 10 days of Rungis Market Simulation

Tableau 5 -

 5 Evolution by extension of the expansive group detected at t=1 77 0,63 1,99 -1,12 MoyTransPrice 5,89 4,52 5,14 5,09 5,84 4,42 MoyFinPrice 6,1 5,36 5,38 5,67 6,25 5,5 MoyStdPrice 5,72 3,42 2,87 3,28 4,84 2,8 MoyProdPrice 3,09 0,88 2,32 2,94 3,48 1,96 MoyAge -5,59 -4,73 -4,77 -5,34 -5,57 -4,91 MoyQuality 6,18 4,96 5,23 5,59 6,16 5,14 MoyQty -1,05 0,63 -0,13 1,03 1,79 1,33 MoymargBuyer 3,76 5 4,14 4,91 4,76 4,77 MoymargSeller 5,29 3,95 4,45 3,97 4,83 3,95

Figure 9 -

 9 Figure 9 -Evolution by extension of the expansive group detected at t=1 (cluster3)

http ://weka.wikispaces.com/


Florence Tubach

Philippe Ravaud

Emilio Martin-Mola

Hassane Awada

Nicholas Bellamy

Claire Bombardier

David T Felson

Najia Hajjaj-Hassouni

Marc Hochberg

Isabelle Logeart

PhD ; Ravaud

PhD ; Martin-Mola

PhD ; Awada

PhD ; Bombardier

PhD Felson Dt

MD, MPH ; Hajjaj-Hassouni N

MD ; Matucci-Cerinic

PhD ; Van De Laar

PhD ; Van Der Heijde

PhD ; Dougados

Florence Pr

Tubach

Minimal Clinically Important Improvement and Patient Acceptable Symptomatic State in pain and function in rheumatoid arthritis, ankylosing spondylitis, chronic back pain, hand osteoarthritis, and hip and knee osteoarthritis: the REFLECT multinational study

est destinée au dépôt

Minimum clinically important improvement and patient

acceptable symptom state in pain and function in rheumatoid arthritis, ankylosing spondylitis, chronic back pain, hand osteoarthritis, and hip and knee osteoarthritis: Results from a prospective multinational study.

Significance and Innovations

 Defining a minimal important improvement or an acceptable state to use in reporting the results of trials or cohorts is widely considered important.

 This article reports the minimal clinically important improvement (MCII) and the patient acceptable symptomatic state (PASS) values estimated for the first time for 4 generic outcomes in 5 chronic rheumatic diseasesrheumatic arthritis, ankylosing spondyitis, hand osteoarthritis, hip and/or knee osteoarthritis and chronic back pain --and 7 countries in the same cohort study.

 The estimates of the MCII and PASS were similar across diseases and countries, whatever the outcome measure. This consistency reinforces the robustness of these concepts and allows for promoting the use of values of MCII (15/100 for absolute improvement, 20% for relative improvement) and PASS (40/100) when reporting the results of trials of the 5 involved rheumatic diseases with pain, patient global assessment, physical function or physician global assessment used as outcome criteria.

Introduction

The choice of an outcome measure is a major step in the design of clinical trials. Trials of rheumatology ideally should systematically evaluate pain, functional impairment, and patient global assessment. For most musculoskeletal disorders, such as osteoarthritis (OA), an outcome criterion is recommended for each of these domains (e.g., for OA, a 0-10 numerical rating scale [NRS] for pain and the Western Ontario and MacMaster Universities Osteoarthritis Index [WOMAC] for functional impairment), [START_REF] Altman | Design and conduct of clinical trials in patients with osteoarthritis: recommendations from a task force of the Osteoarthritis Research Society. Results from a workshop[END_REF][START_REF] Bellamy | Recommendations for a core set of outcome measures for future phase III clinical trials in knee, hip, and hand osteoarthritis. Consensus development at OMERACT III[END_REF][START_REF] Grees | Recommendations for the registration of drugs used in the treatment of osteoarthritis. Group for the respect of ethics and excellence in science (GREES): osteoarthritis section[END_REF]. General outcome criteria, such as patient global assessment, are also recommended.

In general, the results of clinical trials are reported at the group level, for example, with the mean change from baseline by treatment group or with the effect size. Because a statistically significant difference is mostly a matter of sample size and standard deviation of change in the variable, the most difficult issue is determining whether an observed difference is clinically important [START_REF] Wright | Interpreting health-related quality of life scores: the simple rule of seven may not be so simple[END_REF]. How does a mean result at the group level (e.g., a mean change in pain of 1 point on a 0-10 point NRS) translate to a proportion of patients with clinically relevant improvement? To take into account the patient's perspective and for easier interpretation of results, considering whether the change from baseline is important at the individual patient level and then reporting the percentage of patients with improved condition by treatment group is an interesting approach. This approach requires that for each patient, the continuous outcome measure (e.g., change in pain from baseline on a 0-10 NRS) be dichotomized into a binary variable (important change in pain from baseline, that is, change greater than a cut-off defining an important improvement in the patient's perspective). Several cut-offs have been proposed to define the Minimal Clinically Important Difference (or Minimal Important Difference) [START_REF] Beaton | Many faces of the minimal clinically important difference (MCID): a literature review and directions for future research[END_REF][START_REF] Crosby | Defining clinically meaningful change in health-related quality of life[END_REF]. Recently, the Minimum Clinically Important Improvement (MCII) [START_REF] Tubach | Evaluation of clinically relevant changes in patient-reported outcomes in knee and hip osteoarthritis: the Minimal Clinically Important Improvement[END_REF] was proposed, because the amount of change patients consider clinically important is not the same with worsened or improved condition [START_REF] Cella | Meaningful change in cancer-specific quality of life scores: differences between improvement and worsening[END_REF][START_REF] Guyatt | Reassessing quality-of-life instruments in the evaluation of new drugs[END_REF][START_REF] Schwartz | Measurement of fatigue. determining minimally important clinical differences[END_REF]. Because only the patient can perceive whether the change in symptom is important, the definition of the MCII relies on an anchor-based method, with the patient's opinion as the external anchor [START_REF] Crosby | Defining clinically meaningful change in health-related quality of life[END_REF][START_REF] Jaeschke | Measurement of health status. Ascertaining the minimal clinically important difference[END_REF][START_REF] Juniper | Determining a minimal important change in a disease-specific Quality of Life Questionnaire[END_REF].

The MCII reflects the concept of improvement ("feeling better"). The Patient Acceptable Symptomatic State (PASS) has been proposed to address the concept of partial symptomatic remission ("feeling good") [START_REF] Tubach | Evaluation of clinically relevant states in patient-reported outcomes in knee and hip osteoarthritis: the Patient Acceptable Symptom State[END_REF]. These concepts are complementary: a patient's condition can be markedly improved with treatment but can still be poor. For example, a score for pain that decreases from 9 to 5 on a 0-10 NRS is considered a success in terms of improved condition, but an absolute value of 5 still reflects a painful and/or disabling condition. Therefore, reporting the results of a trial in terms of both proportion of patients with improved condition and proportion in an acceptable state at the end of the trial is useful. The PASS value is a clinically relevant cut-off, from the patient's perspective, that allows for classifying patients at the end of the trial as being in "an acceptable state" (i.e., with the outcome score ≤ the PASS)

or not (i.e., with the outcome score > the PASS). The PASS was found more relevant than the MCII because the MCII is the amount of change needed to reach the PASS [START_REF] Tubach | Feeling good rather than feeling better matters more to patients[END_REF]. However, when dichotomizing a continuous variable, reporting the between-groups difference in the proportion of patients with improved condition and the proportion in an acceptable state at the end of the trial is less powerful than reporting the effect size [START_REF] Cohen | The cost of dichotomization[END_REF]. Reporting both is considered an interesting approach [START_REF] Altman | The cost of dichotomising continuous variables[END_REF][START_REF] Brozek | How a well-grounded minimal important difference can enhance transparency of labelling claims and improve interpretation of a patient reported outcome measure[END_REF][START_REF] Dworkin | Core outcome measures for chronic pain clinical trials: IMMPACT recommendations[END_REF][START_REF] Tubach | Minimal clinically important improvement and patient acceptable symptom state for subjective outcome measures in rheumatic disorders[END_REF].

The concepts of MCII and PASS are supported by the Outcome Measures in Rheumatology international network (OMERACT), which has focused 2 meeting sessions (2004 and 2006) on this issue [START_REF] Tubach | Minimal clinically important improvement and patient acceptable symptom state for subjective outcome measures in rheumatic disorders[END_REF][START_REF] Tubach | Minimal clinically important difference, low disease activity state, and patient acceptable symptom state: methodological issues[END_REF] and has contributed importantly to understanding in this area.

To date, MCII and PASS cut-off values have been estimated for hip and knee OA in France [START_REF] Tubach | Evaluation of clinically relevant changes in patient-reported outcomes in knee and hip osteoarthritis: the Minimal Clinically Important Improvement[END_REF][START_REF] Tubach | Evaluation of clinically relevant states in patient-reported outcomes in knee and hip osteoarthritis: the Patient Acceptable Symptom State[END_REF] and acute painful shoulder in France [START_REF] Tubach | Feeling good rather than feeling better matters more to patients[END_REF]. PASS values have been estimated for AS, RA and psoriatic arthritis in Norway [START_REF] Kvamme | Identification of cutpoints for acceptable health status and important improvement in patient-reported outcomes, in rheumatoid arthritis, psoriatic arthritis, and ankylosing spondylitis[END_REF][START_REF] Tubach | Stability of the patient acceptable symptomatic state over time in outcome criteria in ankylosing spondylitis[END_REF]; for RA in The Netherlands [START_REF] Klooster | Patient-perceived satisfactory improvement (PPSI): interpreting meaningful change in pain from the patient's perspective[END_REF]; and AS in the United States and Europe [START_REF] Dougados | Evaluation of the patient acceptable symptom state as an outcome measure in patients with ankylosing spondylitis: data from a randomized controlled trial[END_REF][START_REF] Maksymowych | Thresholds of patient-reported outcomes that define the patient acceptable symptom state in ankylosing spondylitis vary over time and by treatment and patient characteristics[END_REF] and Canada [START_REF] Maksymowych | Evaluation and validation of the patient acceptable symptom state (PASS) in patients with ankylosing spondylitis[END_REF]. This prospective multinational cohort study aimed to estimate the MCII and PASS for the main symptomatic criteria used in rheumatology, for 5 different diseases and in 7 different countries.

Materials and methods

Study design

The Rheumatological Evaluation of Facts Leading to Excellent Treatment (REFLECT) study was a prospective, multicenter, multinational, observational, 4-week cohort study.

Study population

Patients were involved from hospital rheumatology department in Australia, Italy, Lebanon, Morocco, Spain and The Netherlands. To be included in the REFLECT study, outpatients had to be > 18 years old, experience pain from musculoskeletal disease ( 3 on a 0-10 NRS), having been prescribed a nonsteroidal anti-inflammatory drug (NSAID) for the next 4 weeks, and be able to understand the objectives of the study and complete questionnaires in the national language of the country where they live (i.e., English in Australia, French in France, Italian in Italy, Arabic in Lebanon and Morocco, Castillian in Spain, and Dutch in The Netherlands). This study involved patients with one of the following musculoskeletal diseases: RA as defined by the American College of Rheumatology (ACR) [START_REF] Arnett | The American Rheumatism Association 1987 revised criteria for the classification of rheumatoid arthritis[END_REF]; AS as defined by the modified New York criteria [START_REF] Van Der Linden | Evaluation of diagnostic criteria for ankylosing spondylitis. A proposal for modification of the New York criteria[END_REF], with a painful axial involvement at baseline; hand OA as defined by the ACR (29); hip and/or knee OA as defined by the ACR [START_REF] Altman | The American College of Rheumatology criteria for the classification and reporting of osteoarthritis of the hip[END_REF][START_REF] Altman | Development of criteria for the classification and reporting of osteoarthritis. Classification of osteoarthritis of the knee. Diagnostic and Therapeutic Criteria Committee of the American Rheumatism Association[END_REF]; or mechanical back pain (BP) defined as back pain for at least 3 months resulting from a mechanical disorder.

Inclusion could begin with the start of a first NSAID or a switch from one NSAID to another.

Patients were excluded if they had prosthesis on the assessed joint or if they had received an intra-articular injection in the 4 weeks before the study began. All patients initially visited the rheumatologist in charge of the patient, and an NSAID was prescribed (the drug and its dosage was chosen by the physician). A final visit to the same rheumatologist was scheduled 4 weeks later.

Data collection

Data were collected at 2 visits (baseline visit and 4 week-follow-up). At the baseline visit, data were collected on demographics (age, gender, weight, height) and disease characteristics (date of symptom onset for each disease and specific characteristics for each evaluated disease).

At both visits, patients assessed their status regarding their musculoskeletal disease by general and disease-specific patient-reported outcomes (PROs).

-General PROs were 1) pain during the 48 hours before the visit, 2) patient global assessment of disease activity, and 3) functional disability. All these 3 outcomes were measured on a 0-10 NRS.

-Disease-specific PROs, according to each disease, were 1) Health Assessment Questionnaire (HAQ) score for RA [START_REF] Fries | Measurement of patient outcome in arthritis[END_REF], 2) Bath Ankylosing Spondylitis Disease Activity Index (BASDAI) [START_REF] Garrett | A new approach to defining disease status in ankylosing spondylitis: the Bath Ankylosing Spondylitis Disease Activity Index[END_REF] and Bath Ankylosing Spondylitis Functional Index (BASFI) [START_REF] Calin | A new approach to defining functional ability in ankylosing spondylitis: the development of the Bath Ankylosing Spondylitis Functional Index[END_REF] for AS, 3) Australian/Canadian Hand Osteoarthritis Index [START_REF] Kelli | Validity, factor structure, and clinical relevance of the AUSCAN osteoarthritis hand index[END_REF] for hand OA, 4) WOMAC [START_REF] Bellamy | Validation of WOMAC: A health status instrument for measuring clinically important patient relevant outcomes to antirheumatic drug therapy in patients with osteoarthitis of the hip or knee[END_REF] for hip or knee OA and 5) Roland and Morris questionnaire [START_REF] Roland | A study of the natural history of back pain. Part I: development of a reliable and sensitive measure of disability in low-back pain[END_REF] for BP.

At the baseline and final visits, the rheumatologist also rated global assessment of the disease activity (physician global assessment) on a 0-10 NRS (a general physician-reported outcome).

In this article we report the results for the 4 general outcome measures: pain, patient global assessment of disease activity, functional disability and physician global assessment of disease activity. Each of these generic outcomes was measured on a 0-10 NRS (appendix 1), 0 being the best score and 10 the worst. For this report, all outcome scores were normalized to a 0-100 score.

A specific Case Report From (CRF) was available in each of the 6 languages involved. The questionnaires were translated to be as conceptually similar as possible by a single group with experience in translation and linguistic validation of patient questionnaires (see appendix 2).

Outcome measures

As previously described [START_REF] Tubach | Evaluation of clinically relevant changes in patient-reported outcomes in knee and hip osteoarthritis: the Minimal Clinically Important Improvement[END_REF][START_REF] Tubach | Evaluation of clinically relevant states in patient-reported outcomes in knee and hip osteoarthritis: the Patient Acceptable Symptom State[END_REF], to determine the MCII and the PASS, we used an external anchoring method based on patient perspective. The external anchors, collected at the final visit, were as follows: for the MCII determination, the patients assessed the change from baseline on a 3-point Likert scale (improved, no change, or worse). If they reported an improvement, they were asked how important this improvement was to them (very important, moderately important, slightly important or not at all important).

The MCII was estimated for both the absolute (final valuebaseline value) and relative ) difference in each patient-reported outcome. The MCII was defined as the 75 th percentile of the distribution of change in PRO scores for patients who considered they had a slightly or moderately important improvement during the study (i.e., this change in symptoms [the MCII] or a lower change was achieved by 75% of patients with a slightly or moderately important improvement) [START_REF] Tubach | Evaluation of clinically relevant changes in patient-reported outcomes in knee and hip osteoarthritis: the Minimal Clinically Important Improvement[END_REF].

For the PASS determination, patients were asked "If you were to remain for the rest of your life as you were during the last 48 hours, would this be acceptable or unacceptable for you?" with a dichotomous response mode: acceptable or unacceptable. The PASS was defined as the 75th percentile of the distribution of the PRO scores at the final visit for patients who considered their state acceptable at the end of the study [START_REF] Tubach | Evaluation of clinically relevant states in patient-reported outcomes in knee and hip osteoarthritis: the Patient Acceptable Symptom State[END_REF].

The wording of the external anchors for determining the MCII and PASS in the different languages is given in the Appendix 2.

Statistical analysis

Analyses were conducted first by disease and then by country.

Logistic regression was used to model the observations and compute MCII and PASS estimates (according to the 75 th percentile of the logistic curve, as described above) and their 95% confidence intervals (95% CIs). Data for patients with missing data in one of the outcomes (general outcome measures or external anchor) were excluded from the analyses concerning this outcome.

If the number of patients for which the MCII or the PASS had to be estimated was too low (<20), the cut-offs were not estimated for that subgroup. Data for patients lost to follow-up were excluded from the analysis. Statistical analysis involved use of SAS 8.2 (SAS Inst., Cary, NC) and Splus 6.2 (TIBCO Software Inc).

Ethics

This study was conducted in compliance with the Good Clinical Practices and Declaration of Helsinki principles.

This study was an observational study and not a clinical trial. The conduct of the study did not interfere with the usual care of patients. If local law required it, the national principal investigator obtained Institutional Review Board approval of the protocol, as well as approval of all subsequent major changes.

Each patient gave their signed informed consent to be in the study if local law required it.

Results

A total of 1,532 patients were enrolled in the study. The number of patients by disease and country is in Table 1; 1,505 (98%) patients completed the final visit.

Clinical and demographic data are in Table 2. In all, 668 (44.5%) reported a slightly or moderately important improvement, and 1,014 (67.4%) reported being in an acceptable state at the end of the study.

Baseline scores on a 0-10 NRS varied between diseases. Patients with hip or knee OA had the highest scores and those with RA the lowest.

Figures 1 and2 show the estimates of MCII and PASS with their 95% CIs for the 4 general outcome measures. Numerical data are given in the supplemental files.

MCII

For pain in the whole sample, the MCII for absolute difference was -17 (95% CI -18, -15) and the MCII for relative difference was -21% (-24%, -19%).

For patient global assessment in the whole sample, the MCII for absolute difference was -15 (-16, -14) and the MCII for relative difference was -20% (-21%, -18%).

For functional disability assessment in the whole sample, the MCII for absolute difference was -12 (-13, -11) and the MCII for relative difference was -15% (-18%, -13%).

For physician global assessment in the whole sample, the MCII for absolute difference was -14 (-15, -14) and the MCII for relative difference was -18% (-20%, -16%).

PASS

In the whole sample, the PASS was 42 [START_REF] Tubach | The variability in MCII and PASS estimates did not have an impact on treatment effect estimates[END_REF]44) for pain, 43 [START_REF] Sprangers | Which chronic conditions are associated with better or poorer quality of life?[END_REF]45) for patient global assessment, 43 [START_REF] Sprangers | Which chronic conditions are associated with better or poorer quality of life?[END_REF]44) for functional disability assessment, and 39 [START_REF] Roland | A study of the natural history of back pain. Part I: development of a reliable and sensitive measure of disability in low-back pain[END_REF][START_REF] Tubach | The variability in MCII and PASS estimates did not have an impact on treatment effect estimates[END_REF] for physician global assessment.

Globally, the MCII and PASS estimates were consistent across diseases and across countries when the samples were sufficient. Australia, Lebanon, Morocco and The Netherlands often had extreme values and wide confidence intervals because of MCII and PASS estimation for small groups of patients.

For each PRO, patients with RA, who had the lowest baseline scores, had the lowest PASS.

Similarly, patients with hip or knee OA, who had the highest baseline scores, had the highest PASS.

The proportion of patients with both change in score higher than the MCII and final score lower than the PASS for the outcome criteria was 44.6% for pain, 39.8% for patient global assessment, 37.9% for functional disability assessment, and 30.1% for physician global assessment.

Discussion

The importance of defining a minimal important improvement or an acceptable state to use in trials reporting is widely accepted [START_REF] Wells | Minimal clinically important difference module: summary, recommendations, and research agenda[END_REF][START_REF] Wells | MCID/Low Disease Activity State Workshop: Low Disease Activity State in Rheumatoid Arthritis[END_REF]. This prospective multinational study allowed for estimating in the same cohort --globally, by disease, and nationally --the MCII and PASS values for 4 generic outcomes used in clinical trials of 5 chronic rheumatic diseases --RA, AS, hand OA, hip and/or knee OA and chronic BP. Furthermore, this is the first determination of the MCII and PASS for hand OA and chronic BP.

Our estimates of the MCII and PASS were similar across diseases and countries, whatever the outcome measure. Results that deviated from the overall findings could be related to low enrolment (with large 95% CIs). Few patients were involved in this estimation of the MCII and PASS from Australia and The Netherlands as compared with France, Italy and Spain. The consistency of the MCII and PASS estimates reinforces the robustness of these concepts and allows for promoting the same estimates of the MCII and PASS, whatever the disease or country, in reporting the results of trials of any of the 5 rheumatic diseases. Indeed, we need to agree on a single value (for each criterion) that could be used for reporting in all trials.

Previously, values chosen for MCII or PASS cut-offs, among a range of relevant MCII and PASS values, had a low impact on the difference in success rate between arms in a trial (i.e., the proportion of patients benefiting from treatment) [START_REF] Tubach | The variability in MCII and PASS estimates did not have an impact on treatment effect estimates[END_REF]. Thus, the slight differences observed across diseases would not result in a distortion of results. We recommend the use of values of MCII (15/100 for absolute improvement, 20% for relative improvement) and PASS (40/100) in reporting the results of trials using pain, patient global assessment, functional disability or physician global assessment as outcome criteria for any of the 5 involved rheumatic diseases (Table 3).

Musculoskeletal conditions have a great impact on quality of life. Indeed, among patients with a wide range of chronic diseases, those with musculoskeletal diseases report the poorest levels of physical functioning and pain [START_REF] Sprangers | Which chronic conditions are associated with better or poorer quality of life?[END_REF]. The consistency of PASS estimates across diseases allows for using the proportion of patients achieving the PASS for any the 4 generic outcome measures to compare the burden of the different diseases.

Only one study has investigated differences in PASS estimates across countries. This study involved a multinational controlled trial database of AS and compared the United States and different European countries. PASS estimates were higher for English-speaking than non-English-speaking countries [START_REF] Maksymowych | Thresholds of patient-reported outcomes that define the patient acceptable symptom state in ankylosing spondylitis vary over time and by treatment and patient characteristics[END_REF]. We did not include enough English-speaking patients to investigate this issue, which may need further investigation.

Our MCII estimates for global pain and patient global assessment are consistent with those from previous studies, ranging from -15 to -20, [START_REF] Tubach | Evaluation of clinically relevant changes in patient-reported outcomes in knee and hip osteoarthritis: the Minimal Clinically Important Improvement[END_REF][START_REF] Tubach | Evaluation of clinically relevant states in patient-reported outcomes in knee and hip osteoarthritis: the Patient Acceptable Symptom State[END_REF][START_REF] Maksymowych | Thresholds of patient-reported outcomes that define the patient acceptable symptom state in ankylosing spondylitis vary over time and by treatment and patient characteristics[END_REF], and the PASS estimates in this study are slightly higher than that found in most studies, about 35 [START_REF] Tubach | Evaluation of clinically relevant changes in patient-reported outcomes in knee and hip osteoarthritis: the Minimal Clinically Important Improvement[END_REF][START_REF] Tubach | Evaluation of clinically relevant states in patient-reported outcomes in knee and hip osteoarthritis: the Patient Acceptable Symptom State[END_REF][START_REF] Maksymowych | Thresholds of patient-reported outcomes that define the patient acceptable symptom state in ankylosing spondylitis vary over time and by treatment and patient characteristics[END_REF], but lower than that found in a Canadian study, about 50 [START_REF] Maksymowych | Evaluation and validation of the patient acceptable symptom state (PASS) in patients with ankylosing spondylitis[END_REF].

For the 4 outcome measures, the PASS cut-off appears to be slightly higher for hip or knee OA than the other diseases. However, hip or knee OA patients had the highest scores for the 4 outcome measures at baseline. This finding agrees with previous results showing that patients with a chronic disease and with the most severe symptoms would consider satisfactory a higher level of symptoms than would patients with few symptoms [START_REF] Tubach | Feeling good rather than feeling better matters more to patients[END_REF].

The main strength of this study is that it is the first estimation of the MCII and PASS in the same cohort. The methodology used was the same for each disease (5 diseases) and each country (7 countries), with the same wording for the external anchor, treatment, follow-up duration, and outcome measures, which allows for comparisons across diseases and countries.

One limitation of the study is that for each of the 4 outcome measures, we used a discrete 0-10 NRS to measure the outcome and then standardized it into a continuous 0-100 scale. This method relies on the assumption that the transformation was valid in this context. Recent work suggests that transformed WOMAC VAS values are clinically acceptable approximations of native NRS values [START_REF] Bellamy | Comparison of transformed visual analogue and native numerical rating scaled patient responses to the WOMAC Index[END_REF]. Thus, for pain, for instance, the overall MCII value will be -17 if pain is measured on a 0-100 visual analog scale (VAS), but the individuallevel MCII will be -2 if pain is measured on a 0-10 NRS, because the NRS is a discrete variable. Similarly, the PASS value for pain will be 42 with a 0-100 VAS and 4 with a 0-10 NRS.

Another limitation of the study is that it involved patients receiving NSAIDs. Whether similar MCII and PASS estimates would be obtained for patients receiving TNF blockers, for example, is an issue. Patients receiving NSAIDs and biotherapy may have different expectations for treatment efficacy, with a more ambitious goal for biologics, which may affect MCII and PASS estimates. PASS values have been previously estimated in a trial of ankylosis spondylitis patients comparing adalimumab to placebo [START_REF] Dougados | Evaluation of the patient acceptable symptom state as an outcome measure in patients with ankylosing spondylitis: data from a randomized controlled trial[END_REF]. PASS values for patient global assessment were 29.5 (95% CI 20.0, 39,0) for the adalimumab group and 53.0 (32.0, 81.0) for the placebo group. Similar values were observed for total back pain. However, these data were from small samples (88 patients in the adalimumab group and 24 in the placebo group) and the large confidence intervals do not allow for concluding any difference between the 2 groups and compared to our results.

This article reports the MCII and PASS estimates for the 4 generic outcome measures recommended for use in trials of rheumatic disorders. The MCII and PASS estimates for disease-specific PROs will be reported in further articles. The consistency of the MCII and PASS estimates for the generic outcome measures allows for using the same estimate across the 5 conditions studied. Addressing the individual level by describing the proportion of patients with an important improvement (as defined by the MCII) or achieving the PASS cutoff provides additional meaningful information from clinical trials and longitudinal observation studies. This approach will aid in the interpretation of results, add useful information for daily practice [START_REF] Wells | MCID/Low Disease Activity State Workshop: Summary, Recommendations, and Research Agenda[END_REF], and provide information that is complementary to the conventional presentation of results at the group level (mean changes in scores). Reporting the proportion of patients achieving the MCII and the PASS aims at complementing, not replacing, information on the effect size, because the effect size remains a more powerful approach. UCL: upper 95% confidence limit LCL: lower 95% confidence limit UCL: upper 95% confidence limit LCL: lower 95% confidence limit

 Castilian

Piense en todas las formas en las que la (nombre de la enfermedad) le ha afectado durante las últimas 48 horas.

1.

En comparación con el momento de su incorporación al estudio, ¿cómo se ha encontrado en las últimas 48 horas? (Marque "X" sólo una de las casillas)

 mejor  sin cambios (pase a la pregunta 3)

 peor (pase a la pregunta 3)

2.

Si ha respondido "mejor" a la pregunta anterior, ¿qué importancia tiene para usted esta mejoría? (Marque "X" sólo una de las casillas)

 mucha importancia Pensando a tutti i sintomi che (nome della malattia) le ha provocato nelle ultime 48 ore.

 importancia moderada  importancia ligera  ninguna importancia

1.

Come si è sentito durante le ultime 48 ore rispetto a come si sentiva quando ha cominciato qiesto studio? (Barrare con una "X" solo una casella)

 meglio  invariato (vada alla domanda 3)

 peggio (vada alla domanda 3)

2.

Se ha risposto "meglio" alla domanda precedente, quanto è migliorata la sua condizione generale? (Barrare con una "X" solo una casella)

 molto migliorata  moderatamente migliorata  migliorata di poco  migliorata quasi per niente

3.

Se la sintomatologia da lei presentata nelle ultime 48 ore, dovesse rimanere invariata per il resto della sua vita, sarebbe accettabile o inaccettabile per lei? (Barrare con una "X" solo una casella)

 accettabile  inaccettabile

  adress: florence.tubach@bch.aphp.fr Subtitle: MCII and PASS in pain and function in RA, AS, chronic back pain, and hand, hip or knee OA Financial support for the work reported on in the manuscript: This study was supported by Merck, Sharp, & Dohme Chibret Laboratories, France. Dr Bombardier holds a Pfizer Research Chair and a Canada Research Chair in Knowledge Transfer for Muskuloskeletal Care. Other benefits from commercial sources for the work reported on in the manuscript, or any other financial interests that any of the authors may have, which could create a potential conflict of interest or the appearance of a conflict of interest with regard to the work: None Word count: 3,242

(

Figure 1 :

 1 Figure 1: Estimates of the Minimum Clinically Important Improvement (MCII; absolute change and relative change) and 95% confidence intervals (95% CIs) for the 4 general outcome measures. Squares represent point estimates for each disease and each country. The overall MCII estimate (for the whole sample) with its 95% CI is at the bottom of the figure, represented as a larger square. Between brackets is the size of the sample for which the MCII was estimated. Australia and the Netherlands are not represented because of insufficient sample size for estimation.

Figure 2 :

 2 Figure 2: Estimates of the Patient Acceptable Symptomatic State (PASS) with 95% CIs for the 4 general outcome measures. Squares represent point estimates for each disease and each country. The overall PASS estimate (for the whole sample) with its 95% CI is at the bottom of the figure, represented as a larger square. Between brackets is the size of the sample for which the PASS was estimated. Australia and the Netherlands are not represented because of insufficient sample size for estimation.

Table 1 . Number of patients included in the study by disease and country.

 1

	Countries	Disease
		AS	BP	Hand OA	Hip/knee OA	RA	Total
	Australia	5	5	5	5	5	25
	France	95	106	102	109	101	513
	Italy	13	50	64	67	76	270
	Lebanon	40	0	0	44	41	125
	Morocco	33	50	0	50	30	163
	The Netherlands	15	0	0	0	23	38
	Spain	82	78	78	78	82	398
	TOTAL	283	289	249	353	358	1532
	* AS = ankylosing spondylitis; BP = back pain; OA = osteoarthritis; RA = rheumatoid
	arthritis

Table 2 . Baseline characteristics of patients and changes from baseline in patient-reported outcomes.

 2

	AS	BP	Hand OA	Knee/Hip OA	RA
	(n=283)	(n=289)	(n=249)	(n=353)	(n=358)

Table 3 . Recommended values for MCII and PASS for pain, patient global assessment, functional disability and physician global assessment. Important improvement in condition is defined as decrease in score higher than or equal to the MCII cut-off. Acceptable state is defined as a score strictly lower than the PASS cut-off.

 3

	Outcome measured on a 0-100 VAS*	Outcome measured on a 0-10 NRS**

Denk aan alle gevolgen van (naam van de ziekte) in de afgelopen 48 uur.

  Als u voor de rest van uw leven zou moeten blijven zoals de afgelopen 48 uur, zou dit dan aanvaardbaar of onaanvaardbaar voor u zijn? (Zet hieronder in slechts één vakje een " X ") If you answer "improved" at the previous question, how important is this improvement to you? (Mark an "X" in only one box below) If you were to remain for the rest of your life as you were during the last 48 hours, would this be acceptable or unacceptable for you? (Mark "X" in only

	3.
		 aanvaardbaar
		 onaanvaardbaar
		 English
	Think about all the ways your (disease name) has affected you during the last
	48 hours.
	1.	Compared to when you started the study, how have you been during the
		last 48 hours? (Mark "X" in only one box below)
		 improved
		 no change (go to question 3)
		 worse (go to question 3)
	2.
	3.	Si tuviera que permanecer el resto de su vida tal y como ha estado las  very important
		últimas 48 horas, ¿lo consideraría aceptable o inaceptable? (Marque "X" sólo una de las casillas)  moderately important
		 aceptable  slightly important
		 inaceptable  not at all important
	3.	 Dutch:
		one box below)
	1.	Hoe ging het met u de afgelopen 48 uur in vergelijking met toen u aan de  acceptable
		studie begon? (Zet hieronder in slechts één vakje een " X ")  unacceptable
		 beter
		 geen verandering (ga naar vraag 3)  French:
		 slechter (ga naar vraag 3)
	2.	Hoe belangrijk is deze verbetering voor u?
		(Zet hieronder in slechts één vakje een " X ")
		 heel belangrijk
		 tamelijk belangrijk
		 een beetje belangrijk
		 helemaal niet belangrijk

Pensez à toutes les façons dont votre (nom de la maladie) vous a affecté(e)

au cours des dernières 48 heures.

  Si vous avez répondu « mieux » à la question précédente, quelle a été pour vous l'importance de cette amélioration ? (ne cocher qu'une seule case ) Si vous deviez rester tout le reste de votre vie dans le même état que celui des dernières 48 heures, cela serait-il acceptable ou inacceptable pour vous ? (ne cocher qu'une seule case)

		 Très importante
		 Modérément importante
		 Légèrement importante
		 Pas du tout importante
	3.
		 acceptable
		 inacceptable
		 Italian:
	1.	Par rapport au début de l'étude, comment vous êtes-vous senti(e) au cours
		des dernières 48 heures
		 mieux

 pas de changement (allez directement à la question 3)

 moins bien (allez directement à la question 3) 2.

ACKNOWLEDGEMENTS

We acknowledge Gabriel Baron, Arthur Allignol and Blandine Pasquet for their technical help in statistical analysis and reporting, and Dr F Nacci and S Guiducci for their efforts in implementing the study in Italy.

Very good Very bad

Functional disability assessment

Circle the number that best describes the difficulty you had in doing daily physical activities due to your rheumatoid arthritis during the last 48 hours. None Extreme

Physician global assessment

Considering all the ways rheumatoid arthritis has affected your patient during the last 48 hours, circle the number that best describes how he/she has been doing:


Virgile Fritsch
email: virgile.fritsch@inria.fr

Gael Varoquaux

Jean-Baptiste Poline

Bertrand Thirion

Non-parametric Density Modeling and Outlier Detection in Medical Imaging Datasets

Keywords: Outlier detection, non-parametric density estimation, One-Class SVM, Parzen, Local Component Analysis, neuroimaging, fMRI

The statistical analysis of medical images is challenging because of the high dimensionality and low signal-to-noise ratio of the data. Simple parametric statistical models, such as Gaussian distributions, are well-suited to high-dimensional settings. In practice, on medical data made of heterogeneous subjects, the Gaussian hypothesis seldom holds. In addition, alternative parametric models of the data tend to break down due to the presence of outliers that are usually removed manually from studies. Here we focus on interactive detection of these outlying observations, to guide the practitioner through the data inclusion process. Our contribution is to use Local Component Analysis as a non-parametric density estimator for this purpose. Experiments on real and simulated data show that our procedure separates well deviant observations from the relevant and representative ones. We show that it outperforms stateof-the-art approaches, in particular those involving a Gaussian assumption.

Introduction

Group studies based on medical images often attempt to extract representative samples from a given dataset in order to summarize the whole population to one or a few data prototypes. A stronger yet standard assumption is to consider that the data are Gaussian distributed around these prototypes. In the case of neuroimaging, this choice seems convenient because of the high dimensionality of the datasets, but a simple univariate Shapiro-Wilk normality test [START_REF] Shapiro | An analysis of variance test for normality (complete samples)[END_REF] demonstrates that the Gaussian hypothesis is not correct (see e.g. [START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF]). Furthermore, a unimodal distribution hypothesis is inconsistent with approaches that emphasize the impact of population stratification such as genome-wise association studies, or diagnosis settings that imply a separation between patients and healthy subjects. More generally, any parametric characterization of the population statistical structure is challenged by the presence of many outliers related to acquisition or processing issues. Outlier detection and subsequent data cleansing is therefore a first step towards a better understanding of the statistical structure (including between-subjects variability) of medical imaging datasets, which in turn would be of broad interest regarding group analyses. Manifold learning is an alternative solution [START_REF] Gerber | On the manifold structure of the space of brain images[END_REF] useful for visualization, but lacking statistical guarantees.

In a recent contribution [START_REF] Fritsch | Detecting outliers in high-dimensional neuroimaging datasets with robust covariance estimators[END_REF], a regularized version of the standard Mahalanobis distances-based outlier detection method was introduced in this context; it relies on the assumption that inliers are Gaussian distributed. Under mild deviations from this assumption, the approach has been shown to be accurate, making it possible to point out outliers in both high-dimensional and highly polluted neuroimaging datasets. However, although Mahalanobis distances-based approaches can rank the observations accurately (a property that we refer to as accuracy), the choice of a threshold on this ranking is strongly related to the actual data distribution.

There has also been interest in non-parametric algorithms such as One-Class Support Vector Machine (SVM) [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF] for subjects versus patients discrimination [START_REF] Kalatzis | Design and implementation of an SVM-based computer classification system for discriminating depressive patients from healthy controls using the P600 component of ERP signals[END_REF][START_REF] Mourao-Miranda | Patient classification as an outlier detection problem: An application of the one-class SVM[END_REF]. One-Class SVM is well suited for medical image models since the algorithm is computationally efficient and does not rely on any prior distribution assumption. However, it defers part of the work to the practitioner, who has to build a training set with already labeled observations, or directly indicate the amount of contamination in the dataset.

The focus of our work is the analysis of the statistical structure of medical imaging datasets through high-dimensional non-parametric density estimation algorithms. We use an algorithm derived from Parzen windows density estimation (or Kernel Density Estimation (KDE)), that estimate the parameter θ of a given kernel K(.; ., θ) so that, given a learning set (x i ) i=1,..,n , the density probability of an observations x ∈ R p can be written p(x) = 1 n n i=1 K(x; x i , θ). As a new contribution, we subsequently embed this density estimator into a modeseeking procedure to build a simple representation of the data and thus discard potential outliers. This framework provides an easy and efficient way of checking data homogeneity, a feature often required when performing further data analyzes or clinical studies [START_REF] Mourao-Miranda | Patient classification as an outlier detection problem: An application of the one-class SVM[END_REF].

In Section 2 and 3, we briefly describe the different tools that we use in this work and our contributions, that adapt these tools to medical imaging settings. In Section 4, we present some experiments on both simulated and real data. We show that the accuracy of density-based outlier detection is greater than the accuracy of state-of-the-art outlier detection methods. We also demonstrate that, starting from the Local Component Analysis density estimator, we can obtain a simple, yet relevant, differentiation between inliers and outliers. Finally, we present the results of our experiments in Section 5 and discuss them in Section 6.

State-of-the-art methods

Regularized Minimum Covariance Determinant. We consider datasets of n observations x 1 , . . . , x n in R p . Based on given mean µ and covariance Σ parameters, we define the Mahalanobis distance d 2 µ,Σ (x i ) of an observation

x i by d 2 µ,Σ (x i ) = (x i -µ) T Σ -1 (x i -µ).
The larger this quantity, the more likely x i is an outlier. We use the Regularized Minimum Covariance Determinant (RMCD) [START_REF] Fritsch | Detecting outliers in high-dimensional neuroimaging datasets with robust covariance estimators[END_REF] as the reference method for outlier detection and also to whiten the data in distance-based approaches (see Section 4).

One-Class SVM. The One-Class SVM novelty detection algorithm [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF] is a supervised clustering algorithm that relies on a thresholded Parzen windows density estimator to define a frontier between two populations. It can be adapted to the unsupervised problem of outlier detection, but remains a descriptive model that only provides a deviation index for each observation.

Local Component Analysis. Local Component Analysis is another extension

of Parzen windows density estimation where the isotropic assumption inherent in most kernels is relaxed to anisotropic covariance parameters. The θ parameter of the kernel hence becomes the local data covariance matrix Σ, which we estimate using a leave-one-out cross-validation scheme as in [START_REF] Roux | Local component analysis[END_REF]:

Σ * = arg min Σ   - n i=1 log   |Σ| -1 2 (2π) -p 2 1 n -1 j =i exp - 1 2 d 2 xj ,Σ (x i )     . (1)

Theoretical contributions

Setting the LCA regularization term λ. In [START_REF] Roux | Local component analysis[END_REF], an internal regularization term is used to ensure LCA computation stability. The proposed default value is set to λ = 10 -4 [START_REF] Roux | Local component analysis[END_REF]. In our work, we choose λ so that it models properly the central mode of the data. Since outliers may have a large influence on the observed variance of the dataset along some dimensions, we use a robust heuristic: we select the 50% most concentrated observations according to a Parzen windows density estimation, compute the Ledoit-Wolf [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] coefficient shrinkage α from this subsample, and set λ = α 1-α .

Building an interactive outlier detection framework. We propose an efficient procedure to summarize the necessary information about the data structure so that the practitioner can find how many observations to discard: Within the LCA computation, proximity measures of each observation from another are computed as

k ij = exp -1 2 (x i -x j ) T Σ * -1 (x i -x j )
, thus providing a kernelbased representation of the data as a symmetric positive definite matrix K = (k ij ) i,j∈ [1..n] 2 , that summarizes the whole data set structure. Let K = U DU T , where U and D are the matrix of the eigenvectors and diagonal matrix of eigenvalues (σ i , i ∈ [1..n]) of K. Let D δ be the diagonal matrix obtained by shrinking the elements of D by a factor δ ≥ 0 like in [START_REF] Wang | Structural similarity and distance in learning[END_REF]:

D δ (i, i) = 1 -δ σi if σ i > δ, D δ (i, i) = 0 otherwise. D δ yields a shrunk density estimate at each observa- tion g δ (x i ) = e T i U D δ U T e
, where e i is a vector whose entries are 0 except its i-th element which is 1 and e is a vector of ones; note that the normalization constant is omitted as it plays no role in our analysis. We finally define ∆(x i ) = min δ {δ : g δ (x i ) < 0.5}, which associates each observation with the minimal shrinkage value δ that -almost-cancels it. ∆ can be further used to identify different levels of homogeneity amongst the data. Typically, outliers would correspond to a group of observations that vanish with the smallest values of δ, whereas larger δ also trim off regular observations as in Fig. 2.

We define the disappearance function, a ranked version of ∆, as: ∆ rank (i) = ∆(x i ) i:n , where ∆(x i ) i:n is the i-th order value of ∆(x i ). Working with simulated datasets and various values of p, p/n ratios and contamination amount γ, we show that the first knee in the variation of ∆ sort provides a reliable estimation of the number of outliers in the dataset, while no such estimation can be made from the LCA's ranked density function g rank (i) = g(x i ) i:n , where g(x i ) = 1 Z e T i Ke and Z is a normalization factor. As we will show in Section 4, ∆ rank better characterizes data structure than g rank .

Experiments

In our experiments, we first compare One-Class SVM, Parzen density estimation and LCA outlier detection accuracy on both simulated and real data. Parzen density estimation and One-Class SVM are also applied to whitened data (see Section 2), which we refer to as One-Class SVM w and Parzen w . In a second set of experiments, we demonstrate that the subsequent interactive outlier detection framework described in Section 3 provides a usable representation of the data distribution that helps the user to isolate a set of homogeneous samples.

Data description. We generate a γ-polluted n × p dataset by drawing (1γ)n observations from a N (µ, Σ) distribution and γ n observations from a N (µ, αΣ) distribution, α > 1. We also consider Student distributed datasets. In both cases, µ can be set to zero without loss of generality. In all our simulations, we generated the outlier observations so that they can be distinguished from the inliers data.

Our real data were functional MRI contrast maps that were acquired with 3T scanners from multiple manufacturers. Each contrast was available for a number of subject n comprised between 1700 and 2000. BOLD time series was recorded using Echo-Planar Imaging, with TR = 2200 ms, TE = 30 ms, flip angle = 75 • and spatial resolution 3mm isotropic. Standard pre-processing, including slice timing correction, spike and motion correction, temporal detrending, Gaussian smoothing at 5mm-FWHM, and spatial normalization, were performed on the data using the SPM8 software and its default parameters. Functional contrast maps were obtained by fitting a General Linear Model to this data using SPM8 also. For each available subject, we extract the mean signal of every region of interest defined by an anatomical atlas [START_REF] Perrot | Joint bayesian cortical sulci recognition and spatial normalization[END_REF] and concatenate the extracted values so that every subject is represented by a p-dimensional vector in an n× p matrix.

Outlier detection accuracy measures. With simulated data, for various values of p/n, γ, Σ, and α, we build Receiver-Operating Characteristic (ROC) curves [START_REF] Zweig | Receiver-operating characteristic (ROC) plots: A fundamental evaluation tool in clinical medicine[END_REF] of outlier detection accuracy for each method. ROC curves were averaged over 10 runs for each method and each experiment. We compute Area Under Curve values (AUC) that reflect the general outlier detection accuracy of a method for a realistic range of p/n ratios.

Our real datasets are composed of n ∼ 1900 (the exact value depends on the contrast) observations, each described by 113 features (p = 113). With such a p/n ratio, we can construct a fair approximation of the ground truth with covariance-based outlier detection. Computing AUCs as we did with simulated data, and sub-sampling the original dataset, we assess the ability of the different methods to accurately rank the observations by their degree of abnormality.

Results

Quality of outlier detection with density estimators

Gaussian and Student distributed data. On both Gaussian and Student distributed data, the accuracy of outlier detection with LCA dominates the accuracy of the other methods. Generally, all methods perform well with an AUC above 0.9, except when the condition number of the covariance matrix κ(Σ) is above 100. Table 1 illustrates this phenomenon. In the latter case, one has to whiten the data previously to using One-Class SVM and Parzen. The main advantage of LCA is that such a transformation is part of the algorithm.

Real dataset. Fig. 1 shows the accuracy of the different methods on a real neuroimaging dataset, using a contrast related to the perception of angry versus neutral faces. All methods perform well with an AUC above 0.8. LCA achieves the highest accuracy yet, which remains above 0.95 for all p/n ratios. Whitening the data prior to outlier detection with One-Class SVM or Parzen density estimation is relevant since it increases the accuracy of the latter methods by roughly 0.1. Similar results were obtained in five other functional contrasts. 1. AUC values of the different outlier detection methods confronted with variance outliers (Gaussian distributed data, p = 100, γ = 0.4, κ(Σ) = 1000, α = 1.15). Fig. 1. Outlier detection accuracy of nonparametric density estimation algorithms, represented by their AUC (real data). LCA outperforms both Parzen density estimation and One-Class SVM, even applied on whitened data. RMCD parametric method has the same accuracy than the latter. LCA seems to be sensitive to the p/n ratio as its performance decreases with this ratio.

Relevance of an interactive outlier detection procedure

Finding outliers on simulated datasets. We verified with extensive simulations that the first knee of the disappearance function directly provides an estimate of the number of outliers. Fig. 2 illustrates this statement. This result holds for various p, p/n and α values, even though the decision showed to be a bit conservative. This behavior is yet required to guarantee a low false detections rate on heavy tailed distributions such as the Student distribution. Fig. 3 shows that our procedure does not encourage discarding observations when applied to pure Student distributed data.

Investigating the statistical structure of real data. Fig. 4 gives the spectrum of real neuroimaging datasets as obtained from LCA-learned density transformations. Knees can be easily identified in this curve, indicating that two or more relevant groups of observations are present. This observation rank property could not be inferred from the standard decision function. It is noticeable that many observations (about half of the dataset) seem to be suggested as outliers, while looking at a standard bidimensional PCA plot (not shown for the sake of place) would have suggested a much lower number.

Discussion

Statistical modelling of medical images is challenging because of the dimensionality of the data. Most current approaches rely on a Gaussian assumption. Here we use density estimation to rank medical images according to their degree of abnormality in a group study. We demonstrate that outlier detection with Local Component Analysis (LCA) achieves higher accuracy than state-of-the-art methods. This was shown for various p/n settings on both Gaussian and Student distributed data contaminated with up to 40% outliers. Real data experiments showed that LCA accuracy is generally above 0.9, although it seems to slightly decrease in high-dimension. Our choice of the LCA regularization parameter seemed to be optimal in that regard and our experiments demonstrated that LCA should be preferred to other non-parametric methods. Non-parametric methods do not rely on distributional assumptions, but at the expense of explicit statistical control. We propose a simple way to perform outlier detection in an unsupervised framework that does not require any prior knowledge and guides the user in his final decision about how many observations to discard. Because it uses an internal cross-validation scheme, LCA adapts to the data local structure and comes with a natural kernel-based representation of the data. We apply a trace norm penalization to capture the information carried in the kernel matrix which reveals important features on the structure of the data [START_REF] Wang | Structural similarity and distance in learning[END_REF]. As this penalization is the convex relaxation of principal components analysis-based truncation of the kernel matrix, it results in a stable criterion. We used it to characterize the difference between outliers and inliers in a robust procedure. This is meant to provide practitioners a faithful representation of possible inhomogeneities in the population under study. We verified on sev-eral simulated and real functional neuroimaging datasets that this heuristic to chose the regularization of LCA does not yield spurious outlier detections. An attractive generalization of the LCA approach for high-dimensional settings is a mixed model, in which some dimensions are simply modeled as a Gaussian, while others are modeled through equation (1) [START_REF] Roux | Local component analysis[END_REF]. Conclusion. Local Component Analysis was shown to have a good outlier detection accuracy under more general settings than parametric approaches. The interactive outlier detection framework presented in this contribution is of broad interest in medical imaging where manual data screening is impossible because of the high-dimensionality and sample size of the data, and yet essential due to the poor quality of the datasets used in many clinical studies.

This work was supported by a Digiteo DIM-Lsc grant (HiDiNim project, N o 2010-42D). JBP was partly funded by the IMAGEN project, which receives research funding from the E.U. Community's FP6, LSHM-CT-2007-037286. This manuscript reflects only the author's views and the Community is not liable for any use that may be made of the information contained therein.

4 Fig. 2 .Fig. 3 .Fig. 4 .

 4234 Fig. 2. Functions summarizing the data structure from their density. Difference between outliers (red dots) density and inliers (black dots) density only appears in the disappearance function. Choosing δ ≃ 5 yields an outlier detection corresponding to estimating γ ≃ 35% for a real value of γ = 40%.