A Gibon

Á S Ladet

Á G Balent

D Sheeren

Á C Monteil

Agricultural land use Á Integrated landscape assessment Á Landscape scenario Á Agent-based model Á Transdisciplinarity Á Rural development

Keywords: European mountains, Land use and cover change, Landscape multifunctionality

Modelling and simulating change in reforesting mountain landscapes using a social-ecological framework Annick Gibon • David Sheeren • Claude Monteil • Sylvie Ladet • Ge ´rard Balent Abstract Natural reforestation of European mountain landscapes raises major environmental and societal issues. With local stakeholders in the Pyrenees National Park area (France), we studied agricultural landscape colonisation by ash (Fraxinus excelsior) to enlighten its impacts on biodiversity and other landscape functions of importance for the valley socio-economics. The study comprised an integrated assessment of land-use and land-cover change (LUCC) since the 1950s, and a scenario analysis of alternative future policy. We combined knowledge and methods from landscape ecology, land change and agricultural sciences, and a set of coordinated field studies to capture interactions and feedback in the local landscape/land-use system. Our results elicited the hierarchically-nested relationships between social and ecological processes. Agricultural change played a preeminent role in the spatial and temporal patterns of LUCC. Landscape colonisation by ash at the parcel level of organisation was merely controlled by grassland management, and in fact depended on the farmer's land management at the whole-farm level. LUCC patterns at the landscape level depended to a great extent on interactions between farm household behaviours and the spatial arrangement of landholdings within the landscape mosaic. Our results stressed the need to represent the local SES function at a fine scale to adequately capture scenarios of change in landscape functions. These findings orientated our modelling choices in the building an agent-based model for LUCC simulation (SMASH-Spatialized Multi-Agent System of landscape colonization by ASH). We discuss our method and results with reference to topical issues Introduction Natural reforestation of mountain agricultural landscapes is of special concern for society and policy decision-makers in Europe from the local to the regional level [START_REF] Soliva | Envisioning upland futures. Stakeholder responses to scenarios for Europe's mountain landscapes[END_REF]. Many studies have stressed its detrimental impacts on the various environmental and socio-economical services provided by these landscapes of emblematic character and rich biodiversity (Rey [START_REF] Benayas | Abandonment of agricultural land. An overview of drivers and consequences[END_REF]). The role of agricultural abandonment as the major proximate driver of natural reforestation has been well established for several years. However recent works suggest that spatial patterns of agricultural abandonment do not in fact depend as much as previously believed on biophysical environmental conditions but also to a significant extent on social factors [START_REF] Mottet | Agricultural land-use change, mountain landscape dynamics since 1950. A case study in the Pyrenees[END_REF][START_REF] Gellrich | Agricultural land abandonment and natural forest re-growth in the Swiss mountains: a spatially explicit economic analysis[END_REF]. One can therefore assume that studies of land-use and land-cover change (LUCC) that include a spatially-explicit analysis of the relationships between social and ecological systems can provide further insights into changes in mountain landscapes and their drivers, and consequently improve our understanding of the conditions required for their sustainability.

Agent-based models (ABM) are considered to be of major interest for supporting progress in the modelling of interactions between human actors and their environment on a spatially-explicit basis (e.g. [START_REF] Parker | Multi-agent systems for the simulation of land-use and land-cover change. A Review[END_REF][START_REF] Bousquet | Multi-agent simulation and ecosystem management: a review[END_REF][START_REF] Grimm | Pattern-oriented modeling of agent-based complex systems: lessons from ecology[END_REF]. But representing human decisions and biophysical processes in an integrated way in ABM for LUCC simulation (ABM/LUCC) is still a challenge [START_REF] Janssen | Empirically based agent-based models[END_REF][START_REF] Matthews | Agent-based land-use models: a review of applications[END_REF][START_REF] Robinson | Comparison of empirical methods for building agent-based models in land use science[END_REF]). To facilitate this integration in a study of landscape reforestation processes and impacts in the Pyrenees National Park (PNP) area, we relied on the social-ecological system (SES) concept coined by Holling (2001) and we coupled an integrated assessment of past landscape change with the building of an ABM/LUCC. Here we present and discuss the method and current results of the study, which was conducted in close cooperation with landscape stakeholders.

Material and methodology

Research context and orientations

Participatory study on landscape multifunctionality at the initiative of stakeholders In the early 2000s, the PNP scientific service, the agricultural services and rural planning agencies of the Hautes-Pyrenees district considered the continual reforestation of agricultural landscapes with ash (Fraxinus Excelsior) as a serious threat to the future of the natural environment and to valley socioeconomics. They asked our research group to be involved in a study of landscape reforestation and its impacts on biodiversity, agricultural production and the landscape's cultural heritage and aesthetics. The resulting participatory research group (PRG) included scientists in landscape ecology, agricultural sciences, geographical information sciences and computer modelling. The PRG met in plenary workshops once a year, additional meetings being organized as required.

Sustainability sciences recommend making a clear distinction between sustainability of coupled humannatural systems, i.e. their capacity to persist despite uncertainty and change in their environment, and their sustainable development, i.e. the political vision of the desired function of the systems in the future [START_REF] Kates | What is sustainable development? Goals, indicators, values, and practice[END_REF]. The relevance of this view for landscape ecology is acutely recognized [START_REF] Wu | Landscape ecology, cross-disciplinarity, and sustainability science[END_REF][START_REF] Naveh | Landscape ecology and sustainability[END_REF]. It leads to consider landscape sustainability not as an endpoint but as an idea supporting the choice of development directions that will have a positive impact on landscape functions. Material and immaterial landscape functions can be regarded as the output of intended and non-intended impacts of concrete land use [START_REF] Brandt | Multifunctional Landscapes-motives, concepts and perceptions[END_REF][START_REF] Helming | Sustainable development of multifunctional landscapes. Springler, Berlin Holling CS[END_REF]. The operative objectives of our study were thus (1) to gain an integrated understanding of the processes of change in the local landscape/ land-use system regarded as a SES, i.e. a selforganised complex system, and (2) to build a ABM/ LUCC to describe plausible behaviour of the system under a range of alternative scenarios, in order to assist management and policy decision making.

Overview of the study method

Landscape ecology supports an integrated understanding of spatial and temporal changes in ecological systems, their relation to human activities, and their role in the sustainability of landscape systems (e.g. [START_REF] Turner | Landscape ecology: what is the state of the science?[END_REF][START_REF] Otte | Biodiversity at the landscape level. Recent concepts and perspectives for multifunctional land use[END_REF]. Land change science (LCS) addresses interactions among social systems and environmental factors regarding the behaviour of households or firms as the major proximate drivers of land change [START_REF] Turner | emergence of land change science for global environmental change and sustainability[END_REF]. Agricultural sciences consider change in agricultural land-use in relation to farming systems viewed as adaptive socio-technological systems driven by land condition, household behaviour, technology, socio-economic environment and public policy [START_REF] Deffontaines | Agricultural systems and land use patterns: how can we build a relationship?[END_REF][START_REF] Dent | Simulation of ecological, social and economic factors in agricultural systems[END_REF][START_REF] Gibon | Etudier la diversite ´des exploitations agricoles pour appre ´hender les transformations locales de l'utilisation de l'espace: exemple d'une valle ´e des Pyre ´ne ´es centrales[END_REF][START_REF] Bontkes | Modelling the dynamics of agricultural development at farm and regional level[END_REF]. We combined concepts and methods from these scientific domains to address landscape changes in the Pyrenees by means of a case study.

Study area

The study site consisted in four neighbouring villages (42°57 0 N, 0°3 0 W) located about 20 km south of Lourdes, selected as representative of local dynamics [START_REF] Mottet | Agricultural land-use change, mountain landscape dynamics since 1950. A case study in the Pyrenees[END_REF]. The local climate is typical of Atlantic mountains (average annual temperature of 12.5°C and precipitation of 900 mm). Village territories range from 450 to more than 2,000 m asl. Private agricultural land is located between 450 and 1,300-1,400 m asl., common pastures and forests occupying the remaining altitudes. Soils in the cultivated area are relatively deep brunisols [START_REF] Julien | Patterns of ash (Fraxinus excelsior L.) colonization in mountain grasslands: the importance of management practices[END_REF]). Local agriculture is specialised in raising beef cattle and/or sheep. Almost all the usable agricultural area (UAA) is dedicated to grassland, the rest being cultivated for fodder and cereal crops for animal feed. In 1955, there were 120 farms in the four villages; in 2000, only 42 farms remained.

As is the case in most European mountains (Rey [START_REF] Benayas | Abandonment of agricultural land. An overview of drivers and consequences[END_REF][START_REF] Soliva | Envisioning upland futures. Stakeholder responses to scenarios for Europe's mountain landscapes[END_REF], Pyrenean agricultural landscapes were shaped by a very ancient agro-pastoral tradition, and experienced a slow evolution for a long period up to the first half of the twentieth century [START_REF] Gibon | Landscapes on the french side of the Western and Central Pyrenees[END_REF]. The village was (and to a large extent still is) the basic institutional level of their organisation. Ancient communities organised agricultural systems according to a 'house-based' social system [START_REF] Augustins | Les transmissions entre ge ´ne ´rations dans les socie ´te ´s paysannes europe ´ennes[END_REF], including farm holding ('house') transfer to a single heir, who was bound to care it for next generation. The house comprised private land and rights for use of the common land. Every household in the community was allocated a parcel of private land in each agricultural landscape unit (ALU) to enable it to benefit from the full range of village natural resources. Every parcel in a given ALU was subject to similar management decided by the community council. This resulted in close co-organisation between ecological and social systems at respectively three main levels of spatial organisation [START_REF] Balent | Organisation collective et individuelle dans la gestion des ressources pastorales. Conse ´quences sur la durabilite ´agro-e ´cologique des ressources[END_REF]: (1) the ecological plot, the landscape unit, and the whole landscape for ecosystems; (2) the parcel, the individual farmholding, and all village private land for land use. Traditional landscapes in the study area comprised two broad types of units (Fig. 1). Bottom ALUs, located in valley bottoms and on the lowest slopes, included village settlements. Intermediary agricultural landscape units (IALUs) located higher up on valley slopes included scattered barns in which hay was stored and consumed by herds on the spot. The system remained mostly unchanged until the 1950s, when agricultural modernisation and decline in the number of farms led to diversification of farming strategies, and to land-use heterogeneity in the same ALU.

Research methodology

Integrated assessment of landscape and land-use change since the 1950s

Ecological dimension of reforestation processes: The ecological processes of landscape reforestation were studied in Villelongue agricultural landscape (360 ha), which was regarded as a good illustration of regional trends. To assess LUCC at the landscape scale, we built land-cover maps from aerial photographs taken at three dates (1948, 1971 and 2001; French National Geographic Institute) using the method of [START_REF] Muraz | Comment valoriser les photos ae ´riennes dans les SIG[END_REF].

To assess the susceptibility of grasslands to colonisation within the landscape, we mapped all the old trees in the study landscape using aerial photographs taken in 1948 and applied a 100 m circular buffer around each tree according to the literature. [START_REF] Wardle | Fraxinus excelsior[END_REF] observed ash seed dispersal over 125 m from the tree and [START_REF] Wagner | Ein Modell zur Fruchtausbreitung der Esche (Fraxinus excelsior L.) unter Beru ¨cksichtigung von Richtungseffekten[END_REF] showed that the density of the seed rain was high even at long distances (1.33 seeds/m 2 at a distance of 90 m).

The role of land management in ash colonisation was assessed at the parcel level using an empirical model of the relationships between botanical composition, and grassland biomass production and utilization [START_REF] Balent | Construction of a reference frame for studying the changes in species composition in grassland[END_REF]. The model is based on a Correspondence Analysis ordination (CA-Model) of the botanical composition of a reference set of Pyrenean grasslands. The scores of the grasslands on the F1 and F2 axes of the CA model were calibrated with respect to soil fertility using the concentration of N, P and K in plants, which has been shown to be closely correlated to dry matter production [START_REF] Duru | Mineral nutritional status and botanical composition of pasture. I. Effect on herbage accumulation[END_REF], and to biomass utilization from direct measurement of grazing intake and surveys of hay crop yields [START_REF] Balent | Construction of a reference frame for studying the changes in species composition in grassland[END_REF]. The model enables the intensity of use of a parcel to be estimated as the ratio between biomass utilization and biomass production i.e. F2 score/F1 score in the model. We calculated the ordination scores of 96 grasslands randomly selected in the Villelongue territory in the CA-model by averaging the scores of the species they share with the CA-model (passive ordination sensu Økland 1990). Additionally, the presence of ash seedlings was visually observed along two walking transects covering the longest and shortest diagonals of each of the 96 parcels.

Human dimension of land-use change: In order to assess land use (LU) in 2003 and land use change (LUC) since 1950 at the parcel, farm and landscape levels, we combined proven methods based on LCC maps and census data [START_REF] Tasser | An integrative approach for analysing landscape dynamics in diverse cultivated and natural mountain areas[END_REF]) with a field survey of every household that used some agricultural land in the study area, whatever their official profession and landholding area. We built a spatially-explicit survey method using a model of Pyrenean farms as complex SESs [START_REF] Mottet | Agricultural land-use change, mountain landscape dynamics since 1950. A case study in the Pyrenees[END_REF]; Fig. 2). The data collated concerned: (1) the characteristics of the farmholding and the household, the production and land-management systems, (2) land tenure, land cover (LC) and LU practice on each parcel in 2003, (3) change in land tenure and LU of each parcel dating back to 1950, and (4) the history of the farm since 1950. After assessing LUC and its drivers at the landscape level [START_REF] Mottet | Agricultural land-use change, mountain landscape dynamics since 1950. A case study in the Pyrenees[END_REF], we analysed LU and LUC at the farm level in relation to land holding condition and household behaviour.

The diversity among farm households and holdings was analysed from a series of typologies built using selected sets of indicators and manual and multivariate statistical analyses according to [START_REF] Gibon | Etudier la diversite ´des exploitations agricoles pour appre ´hender les transformations locales de l'utilisation de l'espace: exemple d'une valle ´e des Pyre ´ne ´es centrales[END_REF]. Typologies concerned [START_REF] Mottet | Transformations des syste `mes d'e ´levage depuis 1950 et conse ´quences pour la dynamique des paysages dans les Pyre ´ne ´es[END_REF]: (1) the farm structure, the household and its livelihood strategy, (2) the farmholding territory viewed both as a land estate and as a component of the agricultural landscape, (3) how the farming system operated with a focus on its management and on the land-use system. Changes in farm structure and management since 1950 were assessed using the concept of a 'farm trajectory' coined by [START_REF] Capillon | Participatory simulation of land-use changes in the northern mountains of Vietnam: the combined use of an agent-based model, a role-playing game, and a geographic information system[END_REF] form of socio-technological archetypes concerning organisation and change in the farmholding structure and land condition, the land-use system and farm management, and the short-and long-term behaviour of the household.

Simulation of future LUCC scenarios and their impacts on landscape functions

General method for the development of the AMB/ LUCC: The methods used in ABM/LUCC simulation to model individual behaviour and interactions between social actors, and their interactions with the natural environment in the determination of LCC, mainly depend on the orientation of the study [START_REF] Janssen | Empirically based agent-based models[END_REF][START_REF] Robinson | Comparison of empirical methods for building agent-based models in land use science[END_REF]). In simulations for communication and social learning in common property natural resource management (e.g. the 'Companion Modelling' (ComMod) method; [START_REF] Barreteau | Our companion modeling approach[END_REF], the social behaviour and interactions of individual land users are mostly modelled with their cooperation as heuristic rule-based models, under the real-world configuration of landholdings and households (e.g. [START_REF] Capillon | Participatory simulation of land-use changes in the northern mountains of Vietnam: the combined use of an agent-based model, a role-playing game, and a geographic information system[END_REF]. At the other extreme, in simulations conducted for research or policy purposes, population demographics and change in landholdings are represented, and short and long-term decisions made by individual households are mostly modelled using an economic model from the literature and utility functions calibrated from census or/and survey data (e.g. [START_REF] An | Exploring complexity in a human-environment system: an agent-based spatial model for multidisciplinary and multiscale integration[END_REF][START_REF] Millington | An agent-based model of mediterranean agricultural land-use/cover change for examining wildfire risk[END_REF]. Hybrid approaches are now being developed to enhance the account of the variety and mechanisms of household decision-making at the short and long term, in ABMs including modules that allow nesting utility functions, heuristic rule-based techniques, etc. (Manson and Evans 2007; [START_REF] Le | Land-use dynamic simulator (LUDAS): A multi-agent system model for simulating spatio-temporal dynamics of coupled humanlandscape system. I. Structure and theoretical specification[END_REF][START_REF] Fontaine | An agent-based approach to model future residential pressure on a regional landscape[END_REF][START_REF] Valbuena | An agentbased approach to model land-use change at a regional scale[END_REF]. Because of our emphasis on the social dimension of behaviours, we chose the ComMod method to build an ABM/ LUCC that would in the end combine the two approaches.

Framework approach for the development of the conceptual model of SMASH: Relying on hierarchy theory [START_REF] Allen | Hierarchy: perspective for ecological complexity[END_REF][START_REF] Wu | A spatially explicit hierarchical approach to modeling complex ecological systems: Theory and applications[END_REF][START_REF] Moreira | Dynamical coupling of multiscale land change models[END_REF], we assumed that the local landuse/landscape system could be represented as nested modules according to three levels of spatio-temporal organization:

• Module 1 deals with the long-term simulation of the farms' development until the scenario horizon as trajectories of change (e.g. [START_REF] An | Exploring complexity in a human-environment system: an agent-based spatial model for multidisciplinary and multiscale integration[END_REF][START_REF] Le | Land-use dynamic simulator (LUDAS): A multi-agent system model for simulating spatio-temporal dynamics of coupled humanlandscape system. I. Structure and theoretical specification[END_REF][START_REF] Fontaine | An agent-based approach to model future residential pressure on a regional landscape[END_REF]. It is run at a 5-year time step at which changes in individual farm structure and land allocation are simulated from a farm-development behaviour model according to household type and demographic stage, computerised change in the LU and LC condition of parcels (as an output Module 3) and results of interactions between landholders on the land market (as an output of Module 2); • Module 2 deals with the simulation of interactions and transactions between actors on the land sale and rental market (e.g. [START_REF] Parker | A conceptual design for a bilateral agent-based land market with heterogeneous economic agents[END_REF][START_REF] Valbuena | An agentbased approach to model land-use change at a regional scale[END_REF]). It assumes they occur simultaneously at the same 5-year time steps. The interactions and transactions are simulated using policy-delineated areas for land-use allocation, landowner and landholder behaviour models, and contextual opportunities for land sale and rental (in relation with Module 1); • Module 3 deals with the simulation of LCC at the parcel level as a result of the interactions between ecological dynamics and farmers' behaviour in land management. Yearly periods are used to simulate the agricultural production cycle. The conditions of individual farmholding and land management strategies are assumed to remain unchanged during the 5-year period, until their update as an output of Modules 1 and 2.

Five years was selected as the longest time step as an acceptable compromise between ABM simplicity and scenario realism. Module 3 addresses not only land allocation, but also the interactions and feedback between ecosystem dynamics and land-use behaviours. The PRG considered it to be the most critical factor in the future ABM capability to adequately simulate future spatio-temporal patterns of LCC. It was therefore given priority in the development of ABM. Up to now, it is the only module that is fully developed and implemented.

Scenario setting and assessment of landscape functions: The PRG selected the prospective ('whatif') scenario method and decided on the macro-drivers of LUCC to be considered as agricultural and environmental policy, and urbanisation pressure from tourism and populations relocating from towns. Policy stakeholders expressed a preference for extreme assumptions as regards hypothesized changes and their impact on landholder behaviours, but also wanted a realistic simulation at the fine scale of the resulting spatial and temporal patterns of change in the landscape. Landscape change was analysed using 2D maps, 3D virtual images and a set of metric indicators according to a method for multifunctional-landscape assessment that we elaborated with other teams in the European VisuLands project (Miller et al. 2009).

A GIS-based scenario study [START_REF] Gibon | Landscape and land-use change[END_REF] enabled the PRG to refine the view of the building and function of the ABM. Scenarios were worked out manually and mapped in GIS using policy-delineated areas, a basic model of landowner behaviour and a set of broad rule-based behaviour models of farm households built on the field-survey results. The Villelongue agricultural landscape and its condition in 2003 were selected as the scenario baseline, and their horizon was fixed at the year 2030.

ABM/LUCC building and implementation: We used the CORMAS platform [START_REF] Bousquet | Cormas: commonpool resources and multi-agent systems[END_REF] to develop and implement the ABM. The PRG combined (1) the use of the ARDI method (i.e. Actors, Resources, Dynamics and Interactions) developed by [START_REF] Etienne | ARDI: a co-construction method for participatory modelling in natural resources management[END_REF], (2) the knowledge and data provided by the historical assessment, and (3) results of previous transdisciplinary studies conducted on Pyrenean agropastoral systems.

Results of the historical assessment of LUCC

Ecological processes of landscape reforestation with ash

Forested land increased over 20% in the Villelongue landscape from 1948 to 2001, while agricultural land decreased by 24% and built up areas increased by around 9%. The impacts of reforestation were particularly visible in the IALUs (Fig. 3).

The mapping of ash seed rain showed that the density and spatial layout of old ash trees in the landscape (more than 3,000 for about 360 ha) made nearly all the grassland parcels susceptible to natural reforestation (see Supplementary material). This result was confirmed by in situ observation of the presence of ash seedlings in most of the 96 grasslands studied. Therefore, we considered landscape to be a neutral variable regarding the susceptibility of grasslands to colonization by ash. The passive ordination of the 96 grasslands in the CA Model (Fig. 4) showed that the presence of established ash trees was negatively linked to the biomass removal axis (n = 96, t = 5.79, p \ 0.001) and independent of the biomass production level (n = 96, t = 0.86, NS). The colonization process was prevented by hay cutting [START_REF] Julien | Patterns of ash (Fraxinus excelsior L.) colonization in mountain grasslands: the importance of management practices[END_REF]. In grazed-only grasslands and for a given level of grassland production, established ash trees were preferentially found with a low utilization/production ratio. A threshold of intensity of use separating the parcels with and without established ash trees could be eye-adjusted. It corresponded to a 0.5 ratio line, i.e. when half the grass produced annually is used.

In fact, due to the spatial arrangement of ash trees in the landscape, which was inherited from the past, the way the farmers manage the land was found to be the major proximate driver -and current control-of grassland encroachment.

Agricultural LU and its change

The cumulated agricultural area occupied in 2003 by the 40 farms surveyed was about 1,000 ha, i.e. 75% of the traditional area. While farm households owned almost all their farmland in 1950, in 2003, parcels rented under registered and oral agreements amounted respectively to 20 and 27% of their area [START_REF] Mottet | Transformations des syste `mes d'e ´levage depuis 1950 et conse ´quences pour la dynamique des paysages dans les Pyre ´ne ´es[END_REF]. During the course of the study period, agricultural abandonment concerned less than 10% of their parcels, and a return to agricultural use occurred in 70% of their area. With the previous findings of [START_REF] Mottet | Agricultural land-use change, mountain landscape dynamics since 1950. A case study in the Pyrenees[END_REF], these results showed that agricultural abandonment in the study area since the 1950s mainly resulted from the fact that land belonging to collapsed farms had not been taken over when the farmer retired, rather than the intentional abandonment of parcels with high natural constraints in operational farms.

Contemporary farms were found to mostly result from the continuation of an ancient family-household farming tradition (31/40) or the taking over of a farmholding belonging to a retiring neighbour (5/40). After 1950, only four (small) farms were created by families of non-agricultural and/or non-local origin.

In 2003, about 65% of the farmers were more than 50 years old, and only third of them had a certain successor. The most noticeable change in farm characteristics since the 1950s was a large but uneven increase in their size (Table 1). The LU strategies of the households appeared to vary among four archetypes (Table 2). The 'patrimonial strategy' was the most widespread. This involves the patrimonial use of each parcel including scythe mowing on some steep slopes, maintenance of irrigation channels, etc. The 'Retreat strategy' was found among ageing farmers who basically had a patrimonial strategy but who, due to the absence of a successor, were progressively restricting their farming activity until full retirement. The 'selective strategy', which relies on farm technical and economic efficiency rationales, was encountered on part of the biggest farms. The 'niche strategy' was found on very small landholdings kept for social and cultural reasons.

These results revealed the deep-rooted character of traditional cultural values of the farmer population who were attached to patrimonial land management in both the short and long term. The PRG concluded that it has been a positive factor in the preservation of local agroecosystems and landscapes up to now.

Simulation of future LUCC scenarios

Results of the GIS-based study

Three 'what-if' scenarios were built: a 'trend scenario' (Scenario 1) assuming continuation of the current trend; a 'CAP reform' scenario (Scenario 2), to address the impacts of the forthcoming application of the 2003 CAP reform, assuming that farmers will no longer be permitted to use parcels without an official rental agreement; and an 'urbanisation scenario' (Scenario 3), assuming a municipal policy favouring new population settlement (enlargement of the area delineated for urbanisation). The behaviours of individual farm households were modelled according to archetypes as sets of decision rules to comply with hypothesized LU restrictions and to profit from opportunities on the land rental and sale market. The LU changes considered included sale of building land and agricultural abandonment, the 2003 LU being maintained in other cases. Reforestation of abandoned land was simulated from LC transition rules (Cf. Fig. 7). Scenarios resulted in very distinct land-cover patterns in 2030 (Fig. 5) and values of landscape functions [START_REF] Gibon | Landscape and land-use change[END_REF]. Under the assumptions of the simulation and the set of indicators used, the 'trend' scenario was revealed to be the most favourable of the three for the desired future landscape multifunctionality.

Conceptual model and implementation of the first SMASH version

The first ABM version (Module 3 of the whole SMASH) focuses on the simulation of LUCC in grasslands at the parcel level, under constant configuration of the farms. Other LU types (cropland, abandoned land, etc.) are assumed to be unchanged.

In the whole SMASH model, the module will operate for a period of 5 years, at the end of which changes computed in the other ABM modules will initialise the configuration of farms and landscape for the following 5-year period.

The theory-driven view of interactions between grassland management and reforestation processes

The PRG had a common understanding of annual grazing pressure on a parcel being the result of interactions and feedback between (1) the condition of the grassland and LU at the parcel level, (2) the structure and operation of the grassland system at the farm level and (3) the climate conditions of the year concerned. In mountain conditions, fine scale heterogeneity in bio-physical conditions and grassland vegetation communities mean these interactions are very complex [START_REF] Tasser | An integrative approach for analysing landscape dynamics in diverse cultivated and natural mountain areas[END_REF]). This makes them peculiarly difficult to model and simulate in agricultural sciences [START_REF] Balent | Activite ´s de pa ˆturage, paysages et biodiversite ´[END_REF][START_REF] Andrieu | Simulation of forage management strategies considering farm-level land diversity: example of dairy farms in the Auvergne[END_REF]).

To address the relationships between grassland and farmer management from the parcel to the farm level, we used a conceptual model of the Pyrenean grassland system at the farm level built by [START_REF] Gibon | The heterogeneity of grassland fields as a limiting factor in the organization of forage systems. Development of a simulation tool of harvests management in the Central Pyrenees[END_REF]. In this model, the grassland system is regarded as an adaptive SES coupling a set of grasslands and a farmer's decision and action system aimed at meeting year-round and period-specific objectives at various time-steps at the farm and the parcel level. It is represented as a hierarchical system of decision rules for the application of the set of 'technical operations' (TO, e.g. herd grazing, hay cutting and harvesting), that shape the spatial and temporal patterns of production and use of the farm grasslands year round (Fig. 6).

Relationships between grassland production and management were modelled using an empirical model for Pyrenean grasslands built by [START_REF] Duru | Fonctionnement et dynamique des prairies permanentes. Exemple des Pyre ´ne ´es centrales[END_REF] in explicit relation to [START_REF] Balent | Construction of a reference frame for studying the changes in species composition in grassland[END_REF] model. This model enables estimation of production from grass growth curves as a function of daily climate data and the time of the prior TO application, parameterized according to growth cycle and season. It includes rules to account for the impact of altitude on grassland production, and synthesizes the diversity of grass production patterns according to the condition of grassland in the form of calibrated growth curves for high, medium and low productivity grasslands (HP, MP and LP grasslands respectively).

ABM structure and operation

The baseline used to simulate the individual farm's grassland system consisted in (1) the set of grassland

Decision-rules for the spatio-temporal coordination of all the TO management practices at all the LMUs year-round

Year

LMU j provisional TO sequence Fig. 6 Framework conceptual model of a farmer's management strategy of the grassland system at the farm level (adapted from [START_REF] Gibon | The heterogeneity of grassland fields as a limiting factor in the organization of forage systems. Development of a simulation tool of harvests management in the Central Pyrenees[END_REF]. TO i : technical operation of category i (e.g. early spring grazing, hay harvest). LMU j : land management unit number j, , , hierarchical levels in the organisation of a farmer's decision rules from the parcel/ day ( ) up to the whole farm/year ( ) parcels in real-world farm, (2) a virtual cattle herd whose size in 2003 was fitted by converting the number of sheep and cattle into Livestock Units and adding them up, and (3) an attributed rule-based model of management. Social entities (farmers viewed as computer agents) and spatial entities are described at different scales (see Supplementary material for details on the ABM structure). At the farm scale, the farm territory is represented by the set of its farming parcels, i.e. its land-management units, and the set of its cadastral parcels, i.e. the basic units for LU rights and land transactions. An additional category of parcels called elementary parcels was created to enable the model to mimic LU practice on heterogeneous parcels used both for hay cutting and grazing, in which the parts of the parcel that farmers consider to be too steep are usually not mown. Both the cadastral parcels and the farming parcels were modelled as aggregates of elementary parcels. The latter were initialised in CORMAS from a cell grid with a set of attributes imported from a rasterized vector layer built in GIS [START_REF] Monteil | Participatory modelling of social and ecological dynamics in mountain landscapes subjected to spontaneous ash reforestation[END_REF]; Table 3).

The basic time step used to represent the socioecological dynamics in the course of the year is 15 days. Its length results from a long-debated compromise in PRG workshops between an easy-to-run ABM simulation and a realistic representation of the interactions and feedback between agroecological processes and land management decisions.

ABM representation of interactions between LUC and LCC from the farm level to the parcel level: The grassland management practices surveyed in 2003 were narrowed down to three behaviour types: (1) 'close-to-tradition' grassland management attributed to households with a 'traditional' or 'retreat' strategy, (2) 'modernised' grassland management attributed to households with a 'selective' strategy, and (3) 'conservative' grassland management attributed to households with a 'niche strategy'. For the latter, the LU and LC baseline at each parcel were assumed to continue over time. Other management types were modelled as sets of parcel attributes and farmer management rules according to [START_REF] Gibon | The heterogeneity of grassland fields as a limiting factor in the organization of forage systems. Development of a simulation tool of harvests management in the Central Pyrenees[END_REF] model. Only 'close-totradition' grassland management is currently implemented and simulated in the first version of SMASH. The management practice for each TO type is represented in the form of a fixed calendar (Table 4), for now neglecting most of the within-year adaptive rules. We represented the variety of the coupled grassland conditions-projected TO sequences as 'operational LU' categories which consider grassland productivity (HP, MP, LP), altitude and the TO sequence together Classes of accessibility 0 (no information), 1 (road), 2 (track suitable for tractor), 3 (track not suitable for tractor), 4 (no direct access)

Classes of altitude 1 (450-600 m), 2 (600-800 m), 3 (800-1,000 m), 4 (1,000-1,350 m)

Classes of land use 1(cropland), 2 (meadow), 3 (pasture), 4 (abandoned land), 5 (woodland)

Classes of landcover 1(cropland), 2 (grassland), 3 (encroached grassland), 4 (young forest), 5 (mature forest), 6 (building), 7 (other)

Classes of slope 1 (0-10%), 2 (10-30%), 3 (30-50%), 4 ([50%)

Classes * Operational land use combine the type of grassland productivity (HP, MP, LP), the altitude (Alt 1: \ 600 m; Alt 2: 600-1,000 m; Alt 3: [1,000 m), and the projected TO sequence (as a succession of mowing (M), early spring grazing (g) and grazing (G) operations; operations in brackets are optional) (Table 3). Categories are initialised in the ABM by attributing to each parcel an operational LU category based on its condition and TO sequence as collated in the field survey. Grass production at the parcel level is computed in the form of average daily grass growth for each 15 day-time step calibrated according to [START_REF] Duru | Fonctionnement et dynamique des prairies permanentes. Exemple des Pyre ´ne ´es centrales[END_REF] model from a local series of meteorological data (see Supplementary material). In the current ABM version, grass growth simulation does not account for grassland condition and altitude. It is simulated by applying the HP grassland model and a low altitude category. Annual grass consumption through grazing at the parcel level is calculated by cumulating simulated herd consumption during each grazing TO computerized in the course of the year, assuming a daily requirement of 13 Kg Dry Matter (DM) per Livestock Unit. The assessment of undergrazing for each grazed-only parcel is computed at the end of each year by calculating the simulated grass consumption/production ratio and comparing it with the empirical 0.5 threshold. The LUCC processes at the parcel level are modelled as a set of LC transition rules (Fig. 7). The PRG designed them from practical knowledge and the field data collated.

ABM simulation of the grassland system during a one-year cycle: The operation of the grassland system is simulated at a 15-day time step from the time the herd is turned out to graze in spring to when the herd is brought back into winter indoors. TOs are carried out at the 'farming parcel' level according to the farm calendar.

The organization of hay harvest (mowing) TOs is simulated as follows:(1) Selection of the parcels that have to be harvested during a given 15-day time step based on their projected TO sequences over the year;

(2) A mowing objective for the time step is computed to deal with all the meadows concerned during the mowing period; (3) Mowing priority is given to meadows with the highest biomass and easiest access; (4) Mowing implies removing all available fodder from each elementary parcel that makes up the farming parcels concerned, except in heterogeneous parcels where the slope is greater than 30%.

The grazing operations are simulated as follows (see Supplementary material for details): (1) Parcels that can be grazed during a given time step are selected according to projected TO sequences; (2) Grazing priority is given to parcels with the smallest number of TOs already completed; (3) Grazing of additional Table 4 Representation of the projected annual calendar of the grassland management practices for the 'close-to-tradition' management strategy in SMASH parcels is continued until the herd's needs for the 15-day period are satisfied, as long as there are still parcels to graze. At the end of the simulated year, the calculation of annual pressure for each parcel and the updating of reforestation indicators are computed.

Results of the first simulation experiment

The simulation experiment aimed to test the ABM module, with reference to two main directions for ABM validation, i.e. the assessment of their computational and conceptual validity [START_REF] Rykiel | Testing ecological models. The meaning of validation[END_REF]. The ABM was therefore run for 30 years, regardless of the reconfiguration of the farms that will be computed every 5 years once the whole SMASH model is complete.

The results reported here concern two farms (Farms 1 and 3) whose attributes are listed in Table 5. The 'close-to-tradition' model matches their real-world grassland management. The simulated impacts of reforestation at the end of the 30-year period on annual grass production at the farm level, i.e. hay harvested and grazing, are illustrated in Fig. 8a. For Farm 3, stable values for cumulated annual grass production and grass removal were obtained from the 13th year on. From this date, the computed grazing intensity of the grazed-only parcels became systematically higher than the threshold enabling encroachment, thus stopping their reforestation. For Farm 1, where computed annual grassland production largely exceeded herd needs, the process continued throughout the simulation period. Simulated LU and LC of the parcels resulted in 4.3 ha (i.e. 20.7% of the farm UAA) and 23.3 ha (i.e. 32.3% of the farm UAA) of encroached grasslands for Farm 3 and Farm 1 respectively. The general ABM behaviour revealed was consistent with the conceptual model of the reforestation process, but the large values obtained for the encroached areas underlined discrepancies in the simulation of the grassland system for both farms. These discrepancies appear to result mainly from the oversimplification applied in the simulation of grassland production at the parcel level, which will be corrected in the near future. But the higher value obtained for Farm 1 than for Farm 3 could also be due to the conversion of sheep into cattle in the ABM. The impacts of neglecting management differences between species therefore require further examination to determine whether there is a real need for separate representation of the species.

In terms of spatial distribution, simulated reforestation mainly concerned the grazed parcels located at the highest altitudes and the steepest elementary mown parcels, as expected from the conceptual model (Fig. 8b). However the detailed results revealed some gaps in simulated TO sequences for parcels used both for grazing and hay cutting (see Supplementary material). This questions the interest of attempting to simulate LCC at a scale smaller than the cadastral parcel; i.e. the elementary parcel, in the next version of the SMASH model.

Discussion and conclusion

Changes in land-use systems and their impact on landscape functions

The character of agricultural landscapes in the PNP area can be understood as the output of the smooth co-evolution of social and natural systems in the traditional land-use system, and the changes they have undergone since the 1950s as the symbol of its collapse. Our results illustrate the major role played by social institutions in the sustainability of traditional coupled human-natural systems [START_REF] Berkes | Building resilience in lagoon social-ecological systems: a local-level perspective[END_REF]. Until the second half of the twentieth century, the local organisation of landscape/land-use relationships relied on a set of properties, i.e. hierarchical organisation, coordination, complementarity and integration between components, which are known to form the basis of adaptiveness and persistence in complex self-organised systems [START_REF] Kolasa | Ecological systems and the concept of biological organisation[END_REF][START_REF] Balent | Organisation collective et individuelle dans la gestion des ressources pastorales. Conse ´quences sur la durabilite ´agro-e ´cologique des ressources[END_REF]. The acceleration of global dependency and the decrease in local autonomy after the Second World War to a visible break in the continuity with the past in European landscapes [START_REF] Antrop | Why landscapes of the past are important for the future?[END_REF]. It resulted in an increasingly wide array of uncoordinated social actors and bodies (up to the European level) who have direct or indirect authority over rural land use [START_REF] Gibon | Managing grassland for production, the environment and the landscape. Challenges at the farm and landscape level[END_REF]. From the point of view of a SES, the request of the Pyrenean stakeholders in charge of nature conservation and rural policy at the origin of our study can be viewed as an attempt to reintroduce some organization into local land-use systems with respect to the currently expected functions of landscape.

Culture as a pivotal driver of change in the mountain agricultural landscapes of Europe: Our results provide evidence for the interest of a finescale assessment of past organisation in land-use practice and the changes they have undergone since World War II to enhance the understanding of LUCC and its impacts on landscape functions. The spatiallyexplicit survey method we fined-tuned to assess the socio-ecological dimensions of LU change in the study landscape enabled us to overcome the basic difficulty linked to the limited availability of direct data on past land use [START_REF] Tasser | An integrative approach for analysing landscape dynamics in diverse cultivated and natural mountain areas[END_REF]. Our method could be useful in other locations where the actors' recollection of LUC at the small scale is still alive. Insights gained into the social-ecological processes that shape LUCC in local landscapes confirm culture as a driver of landscape change [START_REF] Burgi | Driving forces of landscape change. Current and new directions[END_REF]) in the Pyrenees. Farmers' behaviours strongly imprinted with traditional values have up to now restricted the magnitude of landscape reforestation. However, the results of our GIS-based scenarios indicate that the resilience of local agriculture is being increasingly eroded by economic liberalisation and rural urbanisation, as observed in other European mountains, where agriculture is already vanishing in some locations [START_REF] Lundstrom | Linking models of land use, resources, and economy to simulate the development of mountain regions (ALPSCAPE)[END_REF][START_REF] Streifeneder | Eclairage sur les transformations des structures agricoles dans les Alpes. Comparaison des indicateurs agro-structurels a l'e ´chelle locale[END_REF]. Therefore, unless there is a shift in the direction of socio-economic change and policy, Pyrenean agricultural landscapes will likely undergo a much more drastic change in the near future than they have in recent decades. Drawing on [START_REF] Plieninger | Traditional land-use and nature conservation in European rural landscapes[END_REF], we consider support for the continuation of remaining traditional land-use systems to be a major requirement for the future multifunctionality of European mountain landscapes.

Traps and tracks in the modelling of socialecological processes for the ABM/LUCC simulation of landscape change Our SES framework approach combining a range of methods helped us to smooth over recognised basic difficulties in the representation and parameterisation of a concrete landscape/land-use system [START_REF] Berger | Multi-agent simulation for the targeting of development policies in less-favored areas[END_REF][START_REF] Verburg | Simulating feedbacks in land use and land cover change models[END_REF][START_REF] Robinson | Comparison of empirical methods for building agent-based models in land use science[END_REF]). The transdisciplinary assessment of landscape change produced essential conceptual and empirical knowledge and data to link the land and the people on spatially-explicit bases. The GIS-based scenarios supported a clear definition of the objectives of the ABM simulation and the functionalities expected from the ABM/LUCC. The concentration of ecological and socio-technical field studies in a small geographic area facilitated the parameterization and calibration of the model.

Modelling land-use behaviours and their interactions with the natural environment: The limitations of conceptual models and agent computational abilities to represent the behaviour of LU actors, and the interactions among them and with their environment have hampered the simulation of interactions and feedback in coupled human-natural systems and the capacity of ABM/LUCC simulation to answer realworld questions [START_REF] Janssen | Empirically based agent-based models[END_REF][START_REF] Matthews | Agent-based land-use models: a review of applications[END_REF]). Our study is part of an effort to bridge this gap. Recently published ABM/LUCC works also focus on enhancing the conceptual and empirical modelling of land-users' decision-making processes at the local or regional scale, in the simulation of agricultural-land and/or forest change [START_REF] Le | Land-use dynamic simulator (LUDAS): A multi-agent system model for simulating spatio-temporal dynamics of coupled humanlandscape system. I. Structure and theoretical specification[END_REF][START_REF] Manson | Agent-based modeling of deforestation in southern Yucatan, Mexico, and reforestation in the Midwest United States[END_REF][START_REF] Valbuena | An agentbased approach to model land-use change at a regional scale[END_REF] or urban sprawl [START_REF] Brown | Exurbia from the bottom-up. Confronting empirical challenges to characterizing a complex system[END_REF][START_REF] Fontaine | An agent-based approach to model future residential pressure on a regional landscape[END_REF]. The approaches applied are very similar to each other and to our approach. All rely on complexity science and combine typologies of individual land-user households and trajectories of change in the households/landholdings (in our case family farms) to represent the heterogeneity in social behaviours and the adaptive character of LU decision-making processes in the short and long term.

One main specificity of our work is addressing relationships and feedback between ecosystems dynamics and land management in ABM/LUCC simulation of landscape change.

Modelling cross-scale interactions and feedback between grassland ecosystems and their management at the landscape scale: a tricky challenge: The representation of interactions and feedback between ecosystems and their management in ABM was not only supported by the view of the landscape/landuse system as a SES but also dictated by its contextual conditions and objectives. The account of the interrelationships between grasslands and their management in the grassland-dominated landscape studied was shown to be crucial for simulating reforestation processes and their impacts on landscape functions in a plausible way, which was the stakeholders main request.

To model the interactions and feedback between biophysical features, local grassland ecosystems and their management, from the parcel to the landscape level and from the short to the long term, we used the results of the integrated landscape assessment and also relied on conceptual and empirical models built in prior socioecological and technological research into sustainable development of Pyrenean agriculture. We spent a lot of time and effort in searching for the simplest representation and calibration of these interlinkages. But their complexity in mountain conditions led us to use a very fine basic time step (15 days) in the ABM/LUCC, which is-as far as we know-novel in ABM simulation of landscape change at a time scale of 2-3 decades. Indeed, simplifications are known to be a major stumbling block for the integrated simulation of human decision making and biophysical systems [START_REF] Matthews | The people and landscape model (PALM): towards full integration of human decision making and biophysical simulation models[END_REF]. Our simulation results show that oversimplification could render the ABM incapable of adequately describing the spatiotemporal patterns of reforestation in the local context, and consequently hamper the design of plausible scenarios of change in landscape functions. However, at the current stage of the development of the model, when one only of its three modules is available, it is still too early to fully assess the strengths and weaknesses of our modelling choices.

Conclusion

This study of landscape change was orientated by our strong commitment to the stakeholders of nature conservation and rural development in the Pyrenees National Park area, and our aim to help them 'navigate' (cf. [START_REF] Berkes | Navigating social-ecological systems[END_REF] the landscape/ land-use system. As stressed by [START_REF] Fu | Expanding the bridging capability of landscape ecology[END_REF], transdisciplinarity is an important stimulus for enhancing the bridging capability of research in support of the sustainability of multifunctional landscapes. This led us to invest a significant amount of work in integrated landscape assessment before building an ABM/LUCC to simulate of landscape change in our study area. European mountains are well-known to be 'complicated places' for landscape studies [START_REF] Soliva | Envisioning upland futures. Stakeholder responses to scenarios for Europe's mountain landscapes[END_REF]. We developed a novel method combining knowledge from landscape ecology, land-change and agricultural sciences, and a coordinated set of field studies to elicit the cascade of short-to long-term interactions and feedback in the coupled human-natural system, from the parcel to the whole landscape, that locally shape landscape reforestation processes and patterns and modify landscape functions. Our results provide additional evidence of the special role of mountain agriculture in the creation and maintenance of rural landscapes of high ecological, social and economic values.

Like other recent studies in the literature, our study illustrates the advantages of approaches that rely on complexity sciences and hybrid mixes of methods to improve the representation of land-use decisionmaking processes and their interactions with the natural environment in ABM/LUCC simulation. Typologies of individual (household) behaviours and models of their trajectories of change appear as a general way of accounting for the social heterogeneity and adaptive character of human behaviours in the modelling of interactions and feedback in coupled human-natural systems.

Our study also provides evidence of the need to go beyond changes in LU and CC categories and to model interactions and feedback between the condition and management of ecosystems, the ABM/ LUCC simulation being used to handle the multifunctionality of grassland-dominated landscapes. Results indicate that one way of meeting this special challenge is making available local typologies of grassland ecosystems and their short-and long-term management by farmers, from the parcel to the landscape scale, together with a conceptual model of their interrelationships. We therefore believe that agroecology and farming systems research can make a crucial contribution to enhancing the capacity of landscape ecology to support the multifunctionality and sustainability of mountain and other grasslanddominated landscapes.

  . The diversity encountered among farms was synthesized in the Valley upper slopes -Hayfields and pastures on steeper slopes -Barns for storing hay -Winter feeding of herds in barns -Spring and autumn grazing Altitude commons -Summer grazing Valley bottom and lower slopes -Crops, hayfields and village buildings -Winter grazing of herds and indoor feed supplementation

Fig. 1

 1 Fig. 1 Diagrammatic organisation of agricultural landscape units and land use practice in a traditional Pyrenean agro-pastoral system (adapted from a DDAF-65 drawing)

Fig. 2

 2 Fig. 2 Methodological framework used for the integrated study of family-farm dynamics (from Mottet 2005)

Fig. 3

 3 Fig. 3 Agricultural landscape change between 1948 and 2001 in Villelongue village (Hautes Pyre ´ne ´es) (from Gibon and Fidalgo 2009). Valley commons are in grey. Colour codes for land-cover categories in the agricultural landscape: see legend in the map

Fig. 5

 5 Fig. 5 Prospective LCC of the agricultural landscape in 2030 in Villelongue village according to the scenarios of the GISbased study (from Gibon and Fidalgo 2009). Scenario 1: trend

  of operational land use* P1: HP, Alt 1, MMM(G), P2: MP or LP, Alt 1, MG; P3: MP, Alt 1, MMG; P4: HP, Alt 2, gMMG(G); P5: HP, Alt 2, gM(M)G(G); P6: MP, Alt 2, gMG(G): P7: MP or LP,Alt 2, G(G)G; P8: MP or LP, Alt 3, MG(G); P9: MP or LP, Alt 3

Fig. 7

 7 Fig.7Transition rules of natural resources with and without interaction with land-use practices (adapted from[START_REF] Monteil | Participatory modelling of social and ecological dynamics in mountain landscapes subjected to spontaneous ash reforestation[END_REF]

Fig. 8

 8 Fig. 8 Simulated results of the interactions between land use and land cover change over the 30-year simulation period for grasslands on Farm 1 and Farm 3. a Cumulated year-round grass production (solid lines) and grass removal by hay cutting and herd grazing (dotted lines) at the farm level for Farm 1 (filled square) and Farm 3 (filled triangle), in tonnes DM. b

Table 1

 1 Categories of family farms in the study area in 2003 (adapted from[START_REF] Mottet | Transformations des syste `mes d'e ´levage depuis 1950 et conse ´quences pour la dynamique des paysages dans les Pyre ´ne ´es[END_REF]

	N°Types of family farms

Table 2

 2 Archetypes of farmers' behaviour in land management and farm development in the study area

	Land use strategy	Patrimonial strategy	Selective strategy	Retreat strategy	Niche strategy
	Main drivers	Cultural	Economic	Social	Cultural
	Family socio-	Local or other origin	Local or other origin	Local origin	Local or other
	economics	Part or full-time farming	Part or full-time farming	Late in lifecycle	origin
		Long-term perspective	Long-term perspective	No successor
	Main objectives	Sustainability of	Limitation of labour	Coping with a	Maintenance of a
		agricultural land	input to increase	decrease in	small family-
		resource	livestock production	available labour	farm unit and/
			efficiency		or ''recreation''
					farming
	Farm size and	Medium to large	Medium to large	Small to medium	Small
	farmland spatial pattern	Spread over valley slope-with barns	Large land units and/or few access constraints	Spread over valley slope-with barns	Land units around the
					farmstead
	Animal production	Cattle only or cattle &	Cattle only or cattle &	Cattle & sheep or	Alternative
	system	sheep	sheep	sheep only
	Herd management	More or less close to	'Modernised'	Close to tradition
		tradition
	Land use practice	Maintenance	Abandonment of small	Progressive	Maintenance
		agricultural use of	or constrained land	abandonment of	agricultural use
		each land unit	units	remote land units
	Farmland change	Acquisition of entire	Acquisition of large	Renting or	None
		farmland of retiring	parcels and/or parcels	abandonment of
		farmers	adjacent to existing	excess land units
			ones
	Occurrence	????	??	?	?

Table 3

 3 Imported attributes of the ABM/LUCC cells from which the farming, cadastral and elementary parcels' attributes are derived

	Attribute	Possible values
	Farmer's number	1, 2, 3, 4, 6, 7, 501, 502
	Farming parcel number	Number attributed with the farmer in the survey (Integer [ 0)
	Cadastral parcel number	Number in the official cadastral registrer (Integer [ 0)

Table 5

 5 Characteristics of grassland parcels and herd of Farms 1 and 3FARMOperational land-use categories of grassland parcels (area in ha) Herd size (Livestock units*)

		P1	P2	P3	P4	P5	P6	P7	P8	P9	Total	Cattle	Sheep	Total
	1	6.65	0.86	-	10.94	2.05	0.57	6.35	31.78	12.96	72.16	75	23	98
	3	4.34	1.27	3.63	-	-	2.04	0.46	8.91	-	20,65	34	0	34
	* Livestock Unit = 1 adult cow = 7 adult ewes

Acknowledgments We are indebted to all the farmers and municipalities in the study area for their kind cooperation, to our partners from the agencies and services for and rural development (DDAF-65; CRPGE; CDA-65), particularly to M. Fily, R. Fauge `re, A. Sallent, and P. Caperaa. We warmly thank A. Cabanettes and A. Gavaland (INRA-DYNAFOR) for advising us on ash growth dynamics, M. Etienne (INRA Avignon) for methodological support on the development of the ABM model, C. Simon, A. Mottet, M.P. Julien, and R. Gonzalez for collaboration on this work. We also warmly thank the anonymous referees for their comments which helped us to substantially improve our manuscript. Our project benefited from financial support from INRA, the National Park of the Pyrenees (convention INRA-PNP n°2002-45-S, 2003-2005), the European Union (project QLRT-2001-01017 ''VisuLands'', 2003-2005), the French National Research Agency (projects ADD 'ComMod' and 'Trans', 2005-2008), and the INRA-Region Midi-Pyrenees PSDR programme (project 'Chapay ', 2008-2011).


Based on these observations, numerical models of hydrothermal convection including an inclined basement-sediments interface, a local horst structure, and realistic petrophysical properties have been carried out. The depth-decrease of permeability, the inclination of the interface and the fixed heat flow condition at the base of the model, explain why only a few upwellings can be triggered. Thermal anomalies and a measured temperature profile can be reproduced when fault permeability equals 10 -14 m². Interestingly, structure of convective patterns also exhibits a second hotter upwelling, 8 km east of the Soultz-sous-Forêts upwelling zone, where another geothermal exploration project is now underway. The understanding of hydrothermal convection with realistic fluid and rock properties clearly appears as a predictive tool for geothermal exploration strategies.

Abstract accepted for 2012 AGU Fall meeting, San Francisco.

  Structure of hydrothermal convection in the Upper Rhine Graben as inferred from corrected temperature data and basin-scale numerical models L. Guillou-Frottier 1 , C. Carré 1 , B. Bourgine 1 , V. Bouchot 1 , A. Genter 2 1. BRGM, GeoResources Division, Orléans, France 2. GEIE Exploitation Minière de la Chaleur, Kutzenhausen, FranceGeothermal anomalies in sedimentary basins are strongly controlled by fluid circulation within permeable zones. This study presents a new compilation of newly corrected bottom-hole temperature data in the French part of the Upper Rhine Graben, where the Soultz-sous-Forêts temperature anomaly is locked beneath a local horst structure. After a geostatistically constrained interpolation procedure, maps and cross-sections are extracted from the 3D thermal block, together with the associated standard deviations. Thermal anomalies are preferentially associated with the thickest zones of the permeable fractured Buntsandstein (sandstones) formation, in apparent contradiction with previous models where two major fault zones were suggested to control fluid flow. The underlying fractured granitic basement hosts fluid circulation patterns which are apparently controlled at large-scale by the inclined basement-sediments interface.


Petra Wolf
email: petra.wolf@tum.de

Georg Schmidt

Kurt Ulm

The use of ROC for defining the validity of the prognostic index in censored data

Keywords: Discrimination, Nelson-Aalen estimator, Sensitivity, Specificity

The validity of a diagnostic marker can be summarized using statistical measures either for the goodness of the fit like the deviance, measures of the explained variation like R 2 or the misclassification rate. Other intuitive measures are sensitivity and specificity in the case of binary response. In the absence of censored data the calculation of these measures is widely used. In the presence of censoring the estimation of timedependent sensitivity and specificity is not well known. In this article we propose a new method of calculating ROC curves with censored data using the observed number of events and calculating the additional number of expected events for censored observations. The new method is illustrated with data for predicting mortality in patients surviving a myocardial infarction.

Introduction

The identification of clinical markers to predict a defined event is an important topic in medicine. The event can be the diagnosis of a certain disease, the recurrences of a disease, the response to a certain therapy or mortality. The outcome of these markers may influence further diagnostic procedures or the selection of the appropriate therapy. Good examples can be found in oncology and cardiology, e.g. the PSA-testing for the diagnosis of prostate cancer (which is currently under debate) or the use of leftventricular ejection fraction (LVEF) in predicting mortality after surviving an acute myocardial infarction. The validity of these markers can be assessed with statistical measures either for the goodness of the fit, like the deviance, with measures of the explained variation like R 2 [START_REF] Schemper | Predictive accuracy and explained variation in Cox regression[END_REF][START_REF] O'quigley | Explained variation in proportional hazards regression[END_REF] or the miss-classification rate. However the clinicians prefer measures like sensitivity and specificity, which can be summarized in ROC curves. If the defined event can be assessed within a short period of time, the use of ROC curves in order to present the sensitivity and specificity of a continuous marker is very popular. ROC curves can also be estimated using a combination of markers based on a statistical model. With these curves one can select a cut point in order to divide the marker or the prognostic index into two regions (e.g. positive and negative) based on the choice for sensitivity and specificity depending on the consequences. Therefore ROC curves provide an index of a test' s ability to discriminate between alternative states of health.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

2 However in some studies there is a time lag between the measurement of the markers and the occurrence of the event. This is obvious when the event of interest is death. Some of the patients will die within the follow-up period, but there are always censored observations (patients still alive at the end of the observation period). In such a situation the ROC curves cannot be computed directly and it is difficult to demonstrate the validity of a prognostic marker. However there are some proposals in the literature how to calculate sensitivity and specificity in order to obtain ROC curves even in a study with failure time data [START_REF] Heagerty | Time-dependent ROC curves for censored survival data and a diagnostic marker[END_REF][START_REF] Heagerty | Survival model predictive accuracy and ROC curves[END_REF][START_REF] Cai | The sensitivity and specificity of markers for event times[END_REF]. Unfortunately these extensions are not well-known in the medical community, e.g. in the December issue of the New England Journal of Medicine in 2006 the following statement can be found: "Because standard methods do not exist for deriving ROC curves for time-to-event data ..." [START_REF] Wang | Multiple biomarkers for the prediction of first major cardiovascular events and death[END_REF]. Besides the approaches by [START_REF] Heagerty | Time-dependent ROC curves for censored survival data and a diagnostic marker[END_REF], [START_REF] Heagerty | Survival model predictive accuracy and ROC curves[END_REF] and [START_REF] Cai | The sensitivity and specificity of markers for event times[END_REF] calculating ROC curves, several related methods exist to evaluate the predictive accuracy: an often used and well-known measure is the concordance probability (c-index) [START_REF] Harrell | Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors[END_REF][START_REF] Pencina | Overall C as a measure of discrimination in survival analysis: model specific population value and confidence interval estimation[END_REF]. Another possibility is the analytic calculation of the area under the ROC curve (AUC) [START_REF] Chambless | Estimation of time-dependent area under the ROC curve for long-term risk prediction[END_REF]. If there is no censoring, c-index and AUC are identical and equate the Mann-Whitney statistics P(X>Y). The Mann-Whitney statistic measures the probability that an observation drawn at random from population X (e.g. patients with disease) will exceed an observation drawn at random from population Y (e.g. patients without disease). However under random censorship there are quite differences between these measures: it is shown that Harrell' s concordance index does not meet the Mann-Whitney statistic under censorship [START_REF] Koziol | The concordance index C and the Mann-Whitney parameter Pr(X>Y) with randomly censored data[END_REF]. There are many extensions of Harrell' s c-index, like methods using the inverse probability weighting [START_REF] Liu | A Non-Parametric Approach to Scale Reduction for Uni-Dimensional Screening Scales[END_REF][START_REF] Uno | On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data[END_REF] or with estimating the probability as a function of the Cox model [START_REF] Gönen | Concordance probability and discriminatory power in proportional hazards regression[END_REF]. These approaches got interesting characteristics and forms of application, but in order to use them properly one needs to know the objects they are measuring: the concordance index uses survival time as response whereas the AUC uses the binary endpoint event/no event at or up to a given time point. The Mann-Whitney parameter on the other hand compares two survival curves instead of comparing subjects with events and subjects without event.

Regarding the concordance probability also different definitions are possible: let T be the event time and Z a possible risk score. The conception in [START_REF] Harrell | Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors[END_REF] and [START_REF] Pencina | Overall C as a measure of discrimination in survival analysis: model specific population value and confidence interval estimation[END_REF] is the evaluation of the probability P(Z 1 >Z 2 |T 2 >T 1 ). This definition can be linked to the AUC proposed by [START_REF] Heagerty | Survival model predictive accuracy and ROC curves[END_REF] as a weighted area under the ROC curve. [START_REF] Gönen | Concordance probability and discriminatory power in proportional hazards regression[END_REF] in contrast define a slightly different concordance probability: P(T 2 >T 1 |Z 2 Z 1 ). With this definition it is possible to estimate the probability as a simple function of the Cox model. Another distinction must be made concerning the observation time: some measures quantify a probability at or until a given time, others are not interested in a particular time point and give an overall index for the whole course of time.

For clinical purpose three main aspects can be defined considering the prognostic ability of a marker: first of all the prognostic ability of a marker as addressed with the AUC or the c-index is of interest. The next step is the search for an adequate cut point to differentiate between subjects with high risk, who e.g. should get a treatment and patients with low risk for which a treatment is not essential. To enable a clinician to make such a decision, measures like sensitivity and specificity are useful tools. The information given by a c-index or an AUC alone is not sufficient. Third another useful information can be found in the time-dependence of the AUC. Is the prognostic ability the same over the complete observation time or is the marker only a good measure for some extent of time? To answer these questions, we will present a new and simple method for calculating ROC curves for failure time data. Section 2 describes the estimation of time-dependent sensitivity and specificity in detail. In section 3 we compare the proposed method with the existing method of calculating timedependent ROC curves by [START_REF] Heagerty | Time-dependent ROC curves for censored survival data and a diagnostic marker[END_REF]. Section 4 demonstrates the method on data predicting mortality in patients surviving a myocardial infarction.

Methods

Notation

Assume a continuous marker Z has to be divided into low and high in order to use it for diagnostic or therapeutic purpose. If the disease status D is known for all subjects (D = 0 for no disease, D = 1 with disease) than sensitivity and specificity for a cut point z are defined as

Se (z) = P (Z > z | D = 1) and Sp (z) = P (Z ≤ z | D = 0).
For this definition we assume that the risk of the disease is positively correlated with the marker Z, otherwise the two regions Z > z and Z < z have to be exchanged. In the situation of failure time data, the observation time t has to be taken into account. We use the following notation: T i is the survival time for subject i and we are able to observe only the minimum of T i and C i where C i is the independent censoring time. The follow-up time X i is defined as the minimum of T i and C i (X i = min (T i , C i )) and δ i is the censoring indicator; δ i = 1 represents an event at time

X i (X i = T i ≤ C i ) whereas δ i = 0 indicates a censored observation (X i = C i < T i ).
To denote failure (disease) status of subject i at any time t, we use D

i (t) = 1 if T i ≤ t and D i (t) = 0 if T i > t.
There are two proposals for the calculation of the sensitivity and specificity for time-to-event data [START_REF] Heagerty | Survival model predictive accuracy and ROC curves[END_REF], the cumulative/dynamic and the incident/dynamic version: a) cumulative/dynamic (=C):

The sensitivity for a given cut point z at time t is defined as:

) | ( P ) , ( Se C t T z Z t z ≤ > =
The specificity is given by:

) | ( P ) , ( Sp C t T z Z t z > ≤ =
With this definition all subjects will be used at any fixed time t. Subject i with survival time T i will be counted as control up to time T i and as a case afterwards. b) incident/dynamic (=I):

) | ( P ) , ( Se I t T z Z t z = > = ) | ( P ) , ( Sp I t T z Z t z > ≤ =
Using this definition, the number of subjects used will be decreasing with increasing observation time t. Subjects with a survival time T i play the role of a control for time t < T i and then, for T i = t count as cases. But all subjects with T i < t do not contribute to the calculation of sensitivity/specificity at time t. This definition is similar to the contribution of the observations available at time t to estimate certain parameters in the Cox-model using the partial likelihood method. Both definitions have interesting properties. The cumulative/dynamic definition is useful if all the markers of interest are measured at baseline and one is interested in the prognostic properties of the markers up to time t. The use of the incident/dynamic definition allows time-varying markers [START_REF] Heagerty | Survival model predictive accuracy and ROC curves[END_REF].

In the following we are concentrating on the cumulative/dynamic definition. In the example considered the aim is to estimate sensitivity and specificity up to time t based on measurements obtained at baseline. Therefore we omit the superscript C for the cumulative/dynamic version of sensitivity and specificity.

Estimation

Using Bayes theorem we get

) ( P ) ( P ) | ( P ) | ( P ) , ( Se t T z Z z Z t T t T z Z t z ≤ > ⋅ > ≤ = ≤ > =
Consequently the specificity Sp(z, t) can be estimated as

) ( P ) ( P ) | ( P ) | ( P ) , ( Sp t T z Z z Z t T t T z Z t z > ≤ ⋅ ≤ > = > ≤ =
The problem in both formulas is related to the estimation of the probabilities P(T t | Z > z) and P(T > t | Z z) resp. P(T t) and P(T > t). [START_REF] Heagerty | Time-dependent ROC curves for censored survival data and a diagnostic marker[END_REF] propose to use the Kaplan-Meier method to estimate theses probabilities. However the proportions of patients wit Z > z are varying over time due to censorship and events. Therefore one can not use the result of the Kaplan-Meier curve for the whole sample. For censored observations it does not hold

) ( ) | ( ) ( ) | ( ) ( z Z P z Z t T P z Z P z Z t T P t T P > ⋅ > ≤ + ≤ ⋅ ≤ ≤ = ≤
with fixed proportions P(Z z) resp. P(Z > z), as should be valid under the law of total probability.

In order to obtain the ROC curve at time t one can estimate the expected number of events (e 0 (z,t) and e 1 (z,t) ) with n 0 (z) the number of observations wit Z z and n 1 (z) those with Z > z (n 0 (z)

+ n 1 (z) = n): Z z: ) ( ) | ( ) , ( 0 0 z n z Z t T P t z e ⋅ ≤ ≤ = Z > z: ) ( ) | ( ) , ( 1 1 z n z Z t T P t z e ⋅ > ≤ =
The sensitivity Se(z,t) can be estimated by The other option is to estimate the expected number of events at time t in using all observed events up to time t and calculate the additional number of events in using those who are censored before t.

For censored observations one has to estimate the probability for having an event between the censoring time c and time t (t > c) where the ROC curve should be calculated. This probability P(T t | T > c) can be estimated in using the hazard rates between c and t.

If the censoring occurs between t k-1 and t k , the probability for having an event up to time t is given as

) ( ) ( ) ( ) | ( ) | ( ) | ( ) | ( c S t S t S c T t T c P c T t T t P c T t T c P c T t T P e k k k k - + > ≤ < = > ≤ < + > ≤ < = > ≤ = .
There are some possibilities to estimate the probability

) | ( c T t T c P k > ≤ <
: the exact calculation would be using the definition of a poisson process

1 ) ( ) | ( - - ≈ > ≤ < k k k c t c T t T c P λ . With this definition a
person which is censored shortly after t k-1 has a higher probability of having an event than a person how is censored just before t k . A slightly simplification is to assume that the event has occured in the middle between t k-1 and t k , so that the probability

) | ( c T t T c P k > ≤ < can be estimated by 1 1 2 - - - k k k t t λ .
Under the assumption that the probability for a censored observation having an event is the same if it is censored at the beginning of the interval (t k-1 ; t k ] or at the end of this interval or the interval is short, the probability for having an event reduces to

) ( ) ( ))) ( ) ( ( exp( 1 ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( ) ( ) | ( ) | ( c t c t c S t S c S t S t S c S t S c S c T t T t P c T t T c P e k k k k Λ - Λ ≈ Λ - Λ - - = - = - + - = > ≤ < + > ≤ < =
The approximation results as a first order Taylor approximation and is valid under the condition .

1 )) ( ) ( ( << Λ - Λ c t
For this calculation, as an alternative to the Kalpan-Meier estimator, one can use the Nelson-Aalen estimator for the cumulative hazard rate

= = = Λ ≤ j j j t t j n d t t j λ λ ) ( ) ( ) (
with j d the number of events at time ) ( j t and j n the risk set at )

( j t .
This calculation has to be done for all censored observations before t and has to be performed in both groups separately (Z > z and Z z). With these probabilities together with the observed number of events one gets the expected number of events up to time t and one is able to construct the corresponding table to calculate sensitivity and specificity (Table 1).

Z > z E 1 (t) n 1 -E 1 (t) n 1 (z) Z < z E 0 (t) n o -E 0 (t) n 0 (z) E 1 (t) + E 0 (t) n -(E 1 (t) + E 0 (t)) n ≤ + = t t j j j e d t
) ( E ) ( E ) ( E ) , ( Se 0 1 1 t t t t z + = ; )) ( E ) ( E ( n ) ( E n ) , ( Sp 0 1 0 0 t t t t z + - - =
One characteristic of a ' proper' ROC curve [START_REF] Egan | Signal Detection Theory and ROC Analysis[END_REF]) is its monotonicity. The ROC curve is a monotone increasing function starting at (0,0) and ending at (1,1) [START_REF] Pepe | The Statistical Evaluation of Medical Tests for Classification and Prediction[END_REF][START_REF] Egan | Signal Detection Theory and ROC Analysis[END_REF]. Like with the proposed Kaplan-Meier estimator by [START_REF] Heagerty | Time-dependent ROC curves for censored survival data and a diagnostic marker[END_REF] there can be some situations where this assumption may be violated. There is a proposal in the literature how to solve this problem in using the nearest neighbour estimator (NNE). This method was first described by [START_REF] Akritas | Nearest neighbour estimation of a bivariate distribution under random censoring[END_REF] and adopted for the given situation by [START_REF] Heagerty | Time-dependent ROC curves for censored survival data and a diagnostic marker[END_REF]. All Kaplan-Meier curves near the cut point z are investigated. This approach guarantees monotonicity but it highly depends on the choice of the smoothing parameter, which can result in quite different estimates of sensitivity and specificity. Another approach could be to calculate the ROC curves and apply afterwards a smoothing procedure to assure monotonicity.

One method available could be isotonic regression [START_REF] Salanti | A non-parametric framework for estimating threshold limit values[END_REF]. If the monotonicity is failed than PAVA (Pooling Adjacent Violator Algorithm) has to be applied. With this algorithm the two adjacent points with (Se(z, t), Sp(z, t)) and (Se(z+1, t), Sp(z+1, t)) are pooled together ( ) ( ) A usual measure of the performance of the prognostic power of the marker Z is the AUC, the area under the ROC curve. In studies without censored data the AUC can be estimated in several ways, the Mann-Whitney-U statistic or the trapezoidal rule for integrating a function [START_REF] Pepe | The Statistical Evaluation of Medical Tests for Classification and Prediction[END_REF]. With censored data the trapezoidal rule seems to be appropriate.

+ + + + = + + = 2 ) , 1 ( S )
With the AUC one is able to address two interesting questions. a) Is there a difference between two markers (are two ROC curves different)? b) Is there a change in the prognostic power of a marker over time? For the comparison of two ROC curves their variances and covariances have to be calculated. Several methods are available for the comparison of two ROC curves if the event of interest is known [START_REF] Pepe | The Statistical Evaluation of Medical Tests for Classification and Prediction[END_REF]. In the case of censored data one can use bootstrap methods for calculating variances and confidence intervals. All other methods depend on the disease status which is not known for censored observations.

Comparison with existing methods

The existing method to calculate survival ROC curves via the Kaplan-Meier estimator [START_REF] Heagerty | Time-dependent ROC curves for censored survival data and a diagnostic marker[END_REF] has one basic drawback. Despite the missing guarantee for monotonicity the major problem concerns the fact, that it doesn' t satisfy the condition 0 Se(z,t), Sp(z,t) 1. There can be situations where this condition is violated. Figure 1 presents data that illustrate this problem. These data describe the time to death of 863 kidney transplant patients with age as prognostic factor. A detailed description of the data can be found in [START_REF] Klein | Survival Analysis: techniques for censored and truncated data[END_REF]. In Figure 1 the ROC curve for age at time t=9 years is shown.

< insert figure 1 here >

For this data a cut point beneath an age of 25 leads to a sensitivity above 1. This peculiarity can be explained through the varying probabilities P(Z > z) and P(Z z) as is shown in the Appendix. Due to this inequality the simple approach with the Kaplan-Meier estimator together with the Bayes theorem doesn' t provide a valid estimation of ROC curves in two respects. Besides the common problem with monotonicity particularly the problem with holding 0 Se(z,t), Sp(z,t) 1 necessitates the complex calculation of ROC curves with an approach like the nearest neighbour estimation. A simple smoothing technique like the isotonic regression would fail in this setting because of the missing theoretical justification of this approach. This approach only works without censored observations. To overcome this problem the nearest neighbour estimator was proposed. With this estimator monotonicity and the restriction 0 Se(z,t), Sp(z,t) 1 is guaranteed. Though this causes another problem: how to choose the appropriate bandwidth for smoothing and what is the best bandwidth? There are no clear instructions in literature how to choose the bandwidth. In examining the course of the AUC over time there can be different best bandwidths at different time points. Simulation studies showed two different outcomes (results not shown): in some situations there is a huge difference in AUC for different bandwidth, but in other situations the resulting AUC for different bandwidths is not substantial different. But although the AUC was similar there were partly large differences in the values for sensitivity and specificity for a defined cut point z. A small bandwidth results in few steps in the ROC curve, which is equivalent with few different values of sensitivity and specificity and a ragged curve with big steps. If one chooses a large bandwidth the curve is very smooth with many little steps and therefore many different values for sensitivity and specificity. Using the data of kidney transplant patients there results with a span for the nearest neighbor estimation between 0 and 0.2 with increments of 0.01 an AUC for predicting 9-years-survival between 0.663 and 0.695. An even more different result arises if one only uses the female: there results an AUC between 0.619 and 0.685 for predicting 9-years-survival.

The new proposed method for calculating ROC curves in the context of censored data overcomes these problems. Based on its definition the restriction 0 Se(z,t), Sp(z,t) 1 is naturally achieved. In this example the new ROC curve is similar to the one achieved with the Kaplan-Meier method but it can be seen in Figure 1 that there is a clear restriction to 0 Se(z,t), Sp(z,t) 1. A problem with monotonicity further persist but this can easy be solved with the introduced isotonic regression which is independent of any smoothing parameter and leads to an unique solution for sensitivity, specificity and the AUC.

Application

To demonstrate the new method, data from a study on patients who survived an acute myocardial infarction (MI) are analyzed. This study has been used to derive a new prognostic marker (= DC) in order to predict mortality [START_REF] Bauer | Deceleration capacity of heart rate as a predictor of mortality after myocardial infarction: cohort study[END_REF]. In this publication ROC curves were calculated to compare the prognostic impact of the new markers DC and AC with some of the established markers, e.g. LVEF and SDNN. In the publication the vital status at the end of an observation period of 2 years was used. There are only a limited number of patients with a follow-up less than two years. Therefore the effect of ignoring the censoring is fairly small. In the meantime the follow-up period was extended. We want to demonstrate the prognostic impact of two markers LVEF and DC on the whole sample of n = 2343 patients. The construction of the ROC curve will be first demonstrated using the established marker LVEF. The cut point used for defining the high risk population is z = 30%. The patients with a LVEF of 30% or less are at increased risk. All analyses for this paper were done using R 2.9.2 (R Foundation for Statistical Computing, Vienna, Austria). For the calculation of the expected number of events we used the exact formula via the poisson process as well as the simplification via 1-exp(-(Λ(t)-Λ(c))). The AUC between these methods for LVEF and DC at 2 and 5 years only differed about 0.1%. This is due to the fact that in this study the event times are measured almost daily. Therefore the reported results correspond to the simpler formula which was more feasible to deal with for a bootstrap analysis. Figure 2 shows the Kaplan-Meier curves for both subgroups (LVEF < 30% resp. LVEF > 30%).

< insert figure 2 here > Within five years 181 patients died (39 with LVEF < 30%, 142 with LVEF > 30%) and 1022 observations were censored (40 with LVEF < 30%, 982 with LVEF > 30%). The mortality rate at 5 years among patients with LVEF < 30% was 37.91% and 7.64% among patients with LVEF > 30%.

The expected number of events at t=5 years, calculated as the sum of observed events and additional expected events for censored observations in both groups, results in 46.2 and 170.2 expected events compared to 39 and 142 observed events. These values lead to the following estimates for the sensitivity and specificity of z = 30% at 5 years: the sensitivity is 21.3% and the specificity is 96.5%. In Figure 3 the ROC curve for LVEF at time point 5 years is shown.

< insert figure 3 here >

In order to investigate the change of the impact of LVEF over time the areas under the ROC curves are calculated and plotted against time t (Figure 4).

< insert figure 4 here > It is easy to see that the impact of LVEF for predicting total mortality is very high within the first 6 months after the MI with a value of about 0.80. Afterwards the AUC drops down to values between 0.70 and 0.75.

The new marker DC is also a very powerful predictor. The ROC curve for DC at 5 years is shown in Figure 5. In order to compare both markers the area under each ROC curve is calculated and plotted over time. In Figure 6 the AUCs for both markers are shown. In this graph one can see the change of the impact of both markers over time. LVEF seems to be more powerful within the first 1.5 years.

Afterwards there seems to be a slight advantage for DC. After 2.5 years both markers yield approximately the same AUC with values between 70% and 75%.

< insert figure 5 here > < insert figure 6 here >

We used bootstrap methods to test the difference in AUC between the two markers at time points 2 and 5 years. The mean AUC difference over the 100 bootstrap samples was calculated, and the central 95% of the differences used as a confidence interval. The 95% confidence intervals for the differences were [-0.097; 0.0267] at 2 years and [-0.089; 0.0179] at 5 years. By the bootstrap confidence interval the AUCs from LVEF and DC at 2 and 5 years were not significantly different.

Discussion

With this article a new method of calculating sensitivity and specificity for censored data using the Nelson-Aalen estimator is introduced. When model development is intended to find important prognostic factors the key aspect is how to evaluate the predictive validity. A significant variable in a statistical model alone implies no reliable information about the additive prognostic value. An established discrimination measure is the area under the ROC curve. The AUC has a straightforward interpretation and in studies without censored data the ROC is a standard method. However, as mentioned in the introduction, there is no standard approach for deriving ROC curves for time-to-event data or more generally for the evaluation of the discriminating ability of prognostic markers. There a some proposals in literature, but no one has gained the standard approach.

The method proposed provides a simple straightforward calculation of sensitivity and specificity and the AUC if censoring occurs. For estimating sensitivity and specificity the unknown vital status for censored data at time t is substituted by the expected number of events at time t, using the observed events and calculating the additional expected events for censored observations. To achieve monotonicity for the ROC curve isotonic regression was used after calculating the raw sensitivities and specificities. This method to assure monotonicity displays an intuitive and easy way and is independent of the choice of any smoothing parameter.

Recently a new index was proposed for evaluating the added predictive ability of a new marker, the integrated discrimination index (IDI) [START_REF] Pencina | Evaluating the added predicitve ability of a new marker: From area under the ROC curve to reclassification and beyond[END_REF][START_REF] Chambless | Several methods to assess improvement in risk prediction models: Extension to survival analysis[END_REF]. This index is based on integrated sensitivity and specificity and has a close relationship to the logistic regression R-square [START_REF] Pepe | Comments on 'Evaluating the added predicitve ability of a new marker: From area under the ROC curve to reclassification and beyond[END_REF]. Like the AUC, it can be interpreted as a corrected average sensitivity. The idea for this approach is to overcome the drawback of the AUC in evaluating the added predictive ability of a new marker: to substantially improve the AUC a new marker must have a very large independent influence on the outcome. For this purpose the IDI is a more sensitive index. The increase in AUC is often small between a model with and without the new marker. In these cases the IDI can help to decide on the usefulness of a new marker.

Though an attractive approach, which was presented in this article, there are some more topics which warrant additional research. For our example we used bootstrap to obtain confidence intervals and to test the ability of discrimination of the two markers. But to make strong inference it would be desirable to have an analytic estimator for the variance of the ROC curve and the area under the curve. For evaluating long time periods the use of bootstrap is computationally intensive and it provides no rapid solution for practical application.

Figure 1

 1 Figure 1 ROC curves for 9-years-survival after kidney transplantation with age as prognostic marker. The solid line indicates the curve calculated with the Kaplan-Meier estimator and the dashed line indicates the new estimation with the Nelson-Aalen estimator.

Figure 2

 2 Figure 2 Kaplan-Meier curve of mortality stratified by risk, according to left-ventricular ejection fraction (LVEF).

Figure 3

 3 Figure3ROC curve for prediction of total mortality by LVEF at 5 years.

Figure 4

 4 Figure4AUC(t) for prediction of total mortality based on LVEF with pointwise 95% confidence intervals.

Figure 5

 5 Figure5ROC curve for prediction of total mortality by DC at 5 years.

Figure 6

 6 Figure 6 AUC(t) for prediction of total mortality based on the two markers LVEF and DC.

Figure2A

  Figure2

Figure3A

  Figure3

Table 1 .

 1 Data for calculating sensitivity and specificity at time t

	time t	with event	without event

Acknowledgements

We would like to thank the Editor and a reviewer for helpful comments and suggestions that improved the content of the paper.

Conflict of Interests Statement

The authors declare no conflict of interest.

Appendix

The law of total probability states: Let B 1 , B 2 , …, B n be mutually exclusive events with P(B i )=1 for all i=1,…,n.

In the context of ROC curves with censored data the law of total probability can be formulated as:

. But with censored observations and fixed probabilities P(Z z) and P(Z > z) this equation does not hold, as can be shown with the following data of kidney transplant patients. As an example we use as cut point the median of age (43 years). An extract of the data is shown in the following