Prof Philippe Bousquet

Prof Martyn Plummer

Dr Joël Chadeouf

Delphine Derrien

Anna Bilkova

Nadia Ben Abdallah

Ha Trung

Laurine Claire

Kristan

Laurie

Marion Naoufal

Romain Alison

Mirian Annouck

etc... Valente

Henri

Joe Robert

Bayesian selection model applied to soil carbon dynamics 155

Keywords: Bayesian selection approach, SSVS, spike and slab prior, hierarchical Bayesian model, latent variables, organic carbon dynamics, radiocarbon méthode bayésienne de sélection de variables, recherche stochastique de sélection de covariables, modèle hiérarchique bayésien, variables latentes, dynamique du carbone organique, radiocarbone subject classifications: 35L05, 35L70

First, I would like to express my sincere gratitude to my advisors, Christine HATTÉ and Eric Parent. The first day of our meeting, in the "train bleu" restaurant at Lyon station, stayed etched in my memory. The atmosphere was very relaxed and confortable. In this day, i did not know that the Lyon station will be the departure station to take the train to start my research and have a final destination to defense my thesis successfully. I thank Eric for sharing with me his expertise in Bayesian inference and for all the scientific discussions addressed. These discussions allowed me to gain more experience and perspective in statistics. On the other hand, I mostly appreciated the multidisciplinary aspect of my thesis. I was so lucky to have Christine Hatté as my thesis supervisor. Thank you for helping me discover the new field of soil carbon with finest pedagogy. You surprised me by your reflection on statistics methods which remarkably

The yellow stars highlight a strong collinearity between covariates. The more the correlation tends to -1 (dark blue) or +1 (dark red), the stronger the relationship between the two covariates. Both minimum and maximum temperatures are strongly correlated with the mean annual temperature, with Pearson Correlations (PC) of 0.94 and 0.82 respectively. The maximum precipitation and the seasonal shift precipitation are highly correlated with one another (PC = 0.99) and both are correlated with the mean annual precipitation (PC = 0.86 and 0.81 respectively). . . . . . . . . . . 48 2.9 The F 14 C (upper panels) and carbon content (lower panels) profiles for three profiles from the database: discontinuously sampled (panels a), continuously sampled (panels b and c) with a different (panels b) and the same (panels c) sampling step for both F 14 C and carbon content (panels c). The black lines refer to the soil layer from which the radiocarbon and carbon content were sampled. The blue points indicate the mean level corresponding to the sampled soil layers. . . . . 57 2.10 The carbon content variations (g/kg) according to the soil type for the top (between 0 and 5 cm) (top panel) and deep soil (greater than 80 cm) layers (bottom panel) obtained on the 125 profiles selected for modeling carbon content dynamics. These Box plots underline quartiles: the extreme of the lower whisker, the lower hinge, the median, the upper hinge and the extreme of the upper whisker. 58 5.7 A comparison between the Posterior Inclusion Probability (given in %) for levels within the significant categorical predictors versus the real variation of radiocarbon. The green bars correspond to a Posterior Inclusion Probabilities at least equal to 0.5 (active levels) while the red ones underline the non active levels (PIP < 0.5). The yellow stars in the box-plots indicate the level detected as active by the Bayesian Sparse Group Selection. The boxplots for deep radiocarbon (φ 1 ) are obtained based on the profiles of the database where the depth is greater than 100 cm. In contrast, the boxplots corresponding to the topsoil radiocarbon (φ 2 ) are plotted based on the profiles of the database where the depth is smaller than 10 cm. Luvisol and the natural forest are the baseline levels for soil type and ecosystem type respectively. . . . . . . . . . . . . . . . . . . . . . . . .

5.8

The fusion of levels for ecosystem type for the deep radiocarbon (φ 1 ), topsoil F 14 C (φ 2 ) and the 5.11 The barcharts (right panels) illustrate the Posterior Inclusion Probability (PIP) for soil and ecosystem type in the ω 1 and ω 2 latent linear models. The box-plots (left panels) illustrate the real variation of deep and topsoil soil carbon content according to soil and ecosystem types respectively. The yellow stars indicate the levels detected as active by the Bayesian Sparse Group Selection. Luvisol and natural-forest are the baseline for the soil type and ecosystem type respectively. . . . . . . . . . .

5.12

The results of the ecosystem type levels fusion for each of the three latent variables ω 1 (deep carbon content), ω 2 (topsoil carbon content) and ω 3 (carbon content incorporation). The fusion of levels is based on the Posterior Fusion probability (PFP) for pairs of levels and the posterior median for regression effects. Two levels are fused together if their PFP is at least equal to 50%. . . . . . . . . . . . . . . . . . . . . . . . 5.13 The results of the soil type levels fusion for each of the three latent variables ω 1 (deep carbon content), ω 2 (topsoil carbon content) and ω 3 (carbon content incorporation). The fusion of levels is based on the Posterior Fusion probability (PFP) for pairs of levels and the posterior median for regression effects. Two levels are fused together if their PFP is at least equal to 50%. . . . . . . . 5.14 Distribution of the topsoil (less than 10 cm depth) F 14 C from the database versus the atmospheric F 14 C. Colors highlight sampling years grouped into four periods.: before 1980, [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990], [1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000] 5.20 Distribution of topsoil radiocarbon from the 131 database profiles where the depth does not exceed 10 cm according to the Aridity Index (AI). The boundaries that define various degrees of aridity and the approximate areas involved are given in the table on the right. The more AI tends to 0, the more arid the area is. In contrast an AI higher than 0.65 refers to a humid zone. . . . . . . . . . . . . . 131 5.21 Distribution of the topsoil (less than 10 cm depth) F 14 C from the database profiles versus the seasonal shift (temperature difference between the hottest and the coldest month of the year). The green circle highlights the specific UK sites from the Moor House Nature Reserve "British profiles". 132 LIST OF FIGURES 5.31 Synthetic view of the dependence of soil F 14 C and carbon content on soil-climate-biome. Ten sites from the database were selected as representative of 10 major biomes, taking into account only the explanatory covariates detected as significant for soil F 14 C and soil carbon dynamics respectively.

The orange (the green) band corresponds to the confidence in the local (within site) estimate of F 14 C (soil carbon content), and the gray band corresponds to the between-site variability of soil with similar environmental variables. Dark lines represent the sampled horizon of observed data and the blue points the real F 14 C (soil carbon content) measurements. . . . . . . . . . . . . . . . show the ratio of 14 C to 12 C as a fraction of the atmospheric ratio [START_REF] Bowman | Radiocarbon calibration-current issues[END_REF]. . . .

7.2

The natural variation of radiocarbon production in the last 50,000-year period [START_REF] Hughen | 14c activity and global carbon cycle changes over the past 50,000 years[END_REF] 7.3 IntCal 13 by Reimer et al.(the shift between the 1:1 line and the trend of the calibration is due to the false half-time, the variations around the trend result from changes in production and changes in the

LIST OF TABLES

5

.17 Climatic and environmental conditions for nine F 14 C Brazilian profiles under a "natural forest" from the database. Vegetation cover is reported according to the authors descriptions, "transition" is for the vegetation type at the transition between cerrado and natural forest. Köppen-Geiger subgroup is calculated according to [START_REF] Kottek | World map of the köppen-geiger climate classification updated[END_REF] rule (see Chapter 2") . . . . . . . . . . . . . . . . . .

5.18

The variation of the Net Primary Productivity (NPP) per unit area according to the land type [START_REF] Jackson | A global budget for fine root biomass, surface area, and nutrient contents[END_REF] Plusieurs modèles mécanistes ont été proposés pour exprimer la variation du carbone du sol en fonction de la profondeur. Cependant, un effort pour mieux représenter les simulations de la dynamique du carbone du sol est nécessaire, d'autant plus que la vision des processus physiques est incomplète. Les modèles de carbone du sol les plus utilisés sont:

1. le modèle de RothC [START_REF] Coleman | Simulating trends in soil organic carbon in long-term experiments using rothc-26.3[END_REF] qui intègre également un modèle simple eau/sol, 2. le modèle CENTURY [START_REF] Parton | Analysis of factors controlling soil organic matter levels in great plains grasslands[END_REF] qui inclut des modules plus complexes pour la croissance des plantes, la gestion des opérations, etc.

Le modèle RothC simule que les 30 premiers centimètres du sol et le modèle CENTURY simule les premiers 20 centimètres [START_REF] Falloon | Modelling soil carbon dynamics[END_REF] [START_REF] Houghton | Climate change 1994: radiative forcing of climate change and an evaluation[END_REF].

Les isotopes du carbone permettent de valider la représentation de la dynamique du carbone dans le sol La meilleure façon d'évaluer la performance des modèles mécanistes de la dynamique du carbone du sol est de les comparer avec les données empiriques. Ainsi, une comparaison directe entre les résultats du terrain, du laboratoire, des données et des sorties du modèle mécaniste peut être établie. Pour représenter la matière organique, la spécifier, la suivre et donner une cinétique aux processus, des mesures de la matière organique à plusieurs profondeurs sont nécessaires. Tout d'abord, la quantité du carbone dans le sol peut être définie par les données de teneur en carbone produites par le laboratoire analysant des échantillons prélevés du terrain. Ensuite, il existe des méthodes de traçage isotopique comme les traceurs 13 C et 14 C pour quantifier le temps de résidence de la matière organique du sol, allant de quelques jours jusqu'à plusieurs milliers d'années. La première technique de traçage est fondée sur la surveillance d'abondance du 13 C en cas de changement de végétation (des plantes de type C3 en C4 ou vice versa). Malheureusement, les données disponibles à partir de cette technique ne sont pas en nombre suffisant pour l'évaluation du modèle parce qu'un changement de type de photosynthèse de la végétation est exigé. La deuxième technique, la datation au radiocarbone, est plus puissante. Effectivement, le sol est un témoin des variations des concentrations du radiocarbone de l'atmosphère, en particulier la variation due aux essais nucléaires des profils échantillonnés dans les années 1990) implique que la basse de données n'est pas représentative de la plage de variation du F 14 C atmosphérique surtout par rapport aux variations dues aux essais nucléaires atmosphériques.

Le nombre de mesures de radiocarbone et de la teneur en carbone varie entre 3 et 73 pour les profils échantillonnés.

La majorité des profils de la base de données ont des enregistrements qui varient entre 4 et 10 observations. Ce nombre est satisfaisant pour avoir une idée sur la structure des profils du radiocarbone et de la teneur en carbone.

Une large variation naturelle est observée pour les teneurs en carbone en surface. Cette large variation est en accord avec la base de données mondiale sur le carbone du sol "SoilGrids".

Recherche stochastique de sélection de variables mixtes: application aux variables latentes du modèle hiérarchique de la dynamique du carbone des sols Le deuxième chapitre du manuscrit est présenté sous forme d'un article publié le 13 septembre 2018 dans le Journal de la Société Française de Statistique (SFDS). Dans cet article, nous proposons une approche statistique bayésienne de sélection de variables pour mieux cerner la dynamique du carbone des sols en examinant la variation en profondeur du radiocarbone pour 159 profils sous différentes conditions de climat (température annuelle moyenne, précipitation annuelle moyenne, indice d'aridité, latitude, décalage saisonnier de température, F 14 C atmosphérique) et d'environnement (type de sol, type d'écosystème). Le modèle statistique utilisé dans cet article est inspiré du modèle statistique proposé par Mathieu et al. (2015).

Le modèle hiérarchique non linéaire à variance homogène d'occurrence locale des mesures (Figure 6), s'écrit de la façon hiérarchique suivante: pour un site s ∈ [1 : S], et pour une mesure m s ∈ [1 : m s ], on modélise l'évolution du F 14 C du sol noté par y(s, x) en fonction de la profondeur x par:

y(s, x) = φ 1 (s) + (φ 2 (s) -φ 1 (s)) exp - x φ 3 (s) φ 4 (s) + ε(s, x) ε(s, x) ∼ N(0, σ 2 )
• φ 1 : F 14 C en grande profondeur.

• φ 2 : F 14 C en surface.

• φ 3 : distance relative au point d'inflexion de la courbe.

• φ 4 : décroissance plus ou moins forte.

Les variables latentes φ 1 , φ 2 , φ 3 et φ 4 (s) sont reliées linéairement aux variables potentiellement explicatives de la dynamique du radiocarbone: température annuelle moyenne, précipitation annuelle moyenne, F 14 C atmosphérique, indice d'aridité, décalage saisonnier de température, latitude, type de sol et type d'écosystème. [START_REF] Katsuno | A study of the carbon dynamics of Japanese grassland and forest using 14C and 13C[END_REF]3,[START_REF] Scharpenseel | University of bonn natural radiocarbon measurements vi[END_REF], est le vecteur des effets de régression relative à la variable latente i, E i ∈ R P représente l'effet aléatoire désignant la variabilité inter-sites et X ∈ M S,P (R)) est la matrice de design construite en considérant un contraste traitement.

φ i = Xβ i + E i E i ∼ N(0, σ 2 i I) i = 1, 2 log(φ i ) = Xβ i + E i E i ∼ N(0, σ 2 i I) i = 3, 4
β i = (β i1 , . . . , β iP ) ′ ∈ R P , où i = 1,
d'inclusion a posteriori au moins égales à 0.5 (Figure 7). De plus, pour être sûr que les variables catégorielles non détectées par le SSVS ne sont pas des fausses négatives, on a ajouté les variables catégorielles non détectées d'une manière successive afin de voir si une amélioration du critère DIC (Deviance Information Criterion) peut être établie. Posterior inclusion probability the deep F14C carbon Posterior inclusion probability F14C in topsoil Posterior inclusion probability F14C incorporation

• • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • •
Y = 1µ + G ∑ g=1 X g β g + ε, ε ∼ N(0, σ 2 ) (1)

Spécification des priors:

• Bayesian Group Lasso with Spike and Slab (BGL-SS):

β g |τ 2 g , σ 2 ∼ (1 -π g )N m g (0, σ 2 τ 2 g I m g ) + π g δ 0 (β g ) τ 2 g ∼ G( m g + 1 2 , λ 2
2 )

τ g ∼ Ber(p g ) λ ∼ G(a, b) (2)
Le BGL-SS est une technique qui permet d'estimer et de sélectionner les effets de régression simultanément.

La formulation bayésienne du Lasso a été justifiée par [START_REF] Kyung | Penalized regression, standard errors, and bayesian lassos[END_REF]. [START_REF] Kyung | Penalized regression, standard errors, and bayesian lassos[END_REF] ont montré que le prior double exponentielle proposé par [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] peut être écrit sous forme de la convolution d'une distribution gaussienne sur β g avec un prior Gamma sur son paramètre d'échelle τ g . Pour rendre le modèle plus sparse, [START_REF] Xu | Bayesian variable selection and estimation for group Lasso[END_REF] ont considéré une loi de mélange entre une gaussienne et une masse de Dirac à 0 (pour avoir des effets qui valent exactement 0). Le résultat de la sélection est sensible au choix du paramètre de pénalité. Une petite valeur de λ tend souvent a préféré le modèle nul. Une valeur de 0.5 pour p g g = 1, . . . , G est un choix raisonnable pour faire de la sélection. La règle de décision est basée sur l'estimation a posteriori de la médiane.

• Bayesian Sparse Group Selection (BSGS):

ν g ∼ Ber(p g ) λ lg |ν g ∼ (1 -ν g )δ 0 + ν g Ber(p lg ) β lg |λ lg ν g ∼ (1 -λ lg ν g )δ 0 + λ lg ν lg N(0, τ 2 lg ) (3)
Le BSGS consiste à définir deux indicateurs binaires imbriqués ν g (1: facteur g est sélectionné, 0: sinon) et λ lg (1: la modalité l du facteur g est sélectionnée, 0: sinon). De plus, si le facteur g n'est pas sélectionné dans le modèle (ν g = 0), on affecte des indicateurs nuls pour toutes ces modalités. On peut poser p g = p lg = 0.5, ainsi, tous les sous-modèles sont équiprobables. δ 0 est une masse de Dirac en 0. Le choix de τ lg a un effet important sur la sélection: par exemple une grande valeur de τ lg pour l = 1, . . . ,C g diminue a posteriori la probabilité que le facteur g soit séléctionné. Cette technique nous permet d'avoir un jugement probabiliste non pas seulement sur l'inclusion des variables catégorielles dans le modèle mais aussi sur les effets des modalités au sein d'un même groupe.

• Bayesian Effet Fusion using model-based clustering (BEF):

P(β gl ) = C g ∑ l=0 ν cl N(β gl |µ l , ψ g ) ν l ∼ Dir C g +1 (e 0 ) pour l = 0, . . . ,C g µ 0 = 0 µ l ∼ N(m g , M g ) pour l = 1, . . . ,C g (4)
Cette approche est originale du fait qu'elle permet, non seulement de sélectionner les variables catégorielles significatives pour le modèle, mais aussi de fusionner les modalités au sein du même groupe ayant le même effet sur la variable réponse. Des détails supplémentaires sur le choix des hyperparamètres sont donnés dans l'article de Malsiner (2017). La règle de décision consiste à fusionner les modalités appartenant au même groupe de classification. Une variable catégorielle est éliminée du modèle si toutes ses modalités sont fusionnées avec la modalité de référence.

Les performances de sélection et l'analyse de sensibilité du réglage des hyperparamètres pour la spécification des priors ont été testées pour les trois approches de sélection dans une étude de simulation. Dans cet article, nous présentons également, en détail, la mise en oeuvre des codes sous JAGS pour les trois méthodes de sélection bayésiennes. Pour les profils F 14 C, le meilleur sous-ensemble de facteurs climatiques et environnementaux est identifié par le Bayesian Group Lasso en se basant sur la probabilité d'inclusion a posteriori. En ce qui concerne les profils des teneurs en carbone, le meilleur sous-modèle est identifié à partir du Bayesian Effect Fusion, en se basant sur la probabilité de fusion a posteriori et la médiane des effets de régression.

Interprétations physiques des facteurs climatiques et environnementaux détectés comme significatifs

La deuxième partie de ce chapitre apporte une interprétation physique des facteurs climatiques et environnementaux détectés comme significatifs pour la dynamique du radiocarbone et des teneurs en carbone séparément.

Le radiocarbone atmosphérique, qui aurait dû ressortir parmi les variables explicatives du radiocarbone en surface selon les avis des experts, n'a pas été détecté comme significatif. Ce résultat peut être lié à la sur-représentation des profils échantillonnés en 1990 dans la base de données. D'autre part, la température annuelle moyenne est détectée comme significative pour la teneur en carbone en surface, le F 14 C profond et l'incorporation du F 14 C en profondeur. Notre étude est en accord avec le résultat de [START_REF] Fang | Similar response of labile and resistant soil organic matter pools to changes in temperature[END_REF] qui montre que la matière organique non labile est plus sensible à la température que la matière organique labile.

Les précipitations annuelles moyennes influencent la signature en 14 C et la teneur en carbone en surface ainsi que l'incorporation de ces deux quantités en profondeur. Ces résultats peuvent être liés à la fois à la dilution des composantes des couches superficielles des sols par les composés organiques récemment introduits dans le sol (augmentation de la production primaire résultant de l'augmentation de la MAP) et au priming effect qui entraîne une perte des anciens composés organiques du sol. D'autres interprétations physiques portent sur l'indice d'aridité et le décalage saisonnier de température .

Une surestimation des variables latentes correspondant à l'incorporation du radiocarbone et des teneurs en carbone en profondeur est identifiée. Cette mauvaise estimation, loin des valeurs qu'on peut avoir en réalité, peut être expliquée par la complexité du modèle et au lien non linéaire entre ces variables latentes et les réponses du radiocarbone et de la teneur en carbone.

Le résultat de la fusion des types de sol pour le radiocarbone en surface et pour celui en profondeur souligne que le profil de 14 C est davantage dominé par le climat/la végétation et la texture du sol pour les premiers centimètres du sol et par la teneur en argile pour les couches les plus profondes.

Prédictions du modèle statistique dans un contexte de changements globaux L'avantage du modèle statistique est d'être utilisé pour prédire des profils en F 14 C et des teneurs en carbone dans des endroits où aucune donnée n'est disponible. Ici, en particulier on a essayé de prédire les profils du radiocarbone et de la teneur en carbone dans le cas de la conversion des forêt équatoriales en terres cultivées. Cette étude se base sur neuf profils localisés au Brésil. Une augmentation significative du radiocarbone en profondeur de 0.45 à 0.58, est observée pour les couches profondes. Autrement dit, le temps de résidence du carbone dans les couches profondes est plus long pour la forêt tropicale humide que pour les terres cultivées. Mais aucun changement n'est révélé pour la teneur en carbone en profondeur. Ces résultats sont conformes à ceux de Balesdent et al. (2018) qui montrent que l'utilisation des sols pour les cultures réduit l'incorporation de carbone dans la couche superficielle du sol, mais pas dans les couches plus profondes. Nos résultats ajoutent d'autres éléments à la discussion de Balesdent et al. (2018), allant au-delà du fait que la matière organique de nos sols actuels est l'héritage de sa gestion par plusieurs générations d'agriculteurs.

Par ailleurs, cette partie contient aussi une étude sur le reboisement des terres cultivées et prairies dans les régions tempérés. Les résultats montrent que le reboisement des terres cultivées et des prairies tempérées entraîne une augmentation des stocks de carbone à court et à long terme.

Une étude aussi à été faite pour étudier l'augmentation de la température de 1, 1.5 et 2 • C sur la dynamique des profils du radiocarbone et celle de la teneur en carbone.

Conclusion et perspectives

Le dernier chapitre de ma thèse se divise en deux parties: la première résume brièvement les points principaux de la thèse partant du défi scientifique du départ et la mauvaise compréhension de la dynamique du carbone dans le sol jusqu'à l'avantage de la modélisation statistique. L'extrapolation du modèle statistique bayésien développé pour la teneur en carbone des sols est utile pour obtenir une estimation globale (ou régionale) du stock de carbone des sols. Les modèles statistiques bayésiens pour la dynamique de la teneur en carbone et du radiocarbone nous permettent de prédire les profils du contenu en carbone et en carbone d'un nouveau site, en connaissant les informations climatiques et environnementales correspondantes. Comme aucune mesure n'est fournie pour ce site, les intervalles de crédibilité des paramètres inconnus du modèle seront plus larges que ceux observés pour les sites échantillonnés. En premier lieu, les profils prédits, lors du changement d'utilisation des sols ou des conditions climatiques, sont obtenus sans tenir compte des mesures observées (chapitre 5.5.2, section 5.5). Cela signifie que le site est considéré comme, un nouveau uptake, have the advantage to better ensuring food security by preserving soil fertility. In that respect, a better understanding of the significance of the soil carbon pool was reached after the 2013 IPCC report as for the first time soils were considered as one of the resources for climate change mitigation. It is also worth mentioning the "Four per Thousand" initiative (https://www.4p1000.org) which aims at increasing the world soil carbon sequestration to a 40 cm depth at the rate of 0.4% per year in order to mitigate the global issues of climate change, food insecurity, and environmental pollution [START_REF] Lal | Beyond cop 21: potential and challenges of the "4 per thousand" initiative[END_REF]. To distinguish between sequestration and storage, it is commonly established that sequestration should be sustainable (at least 100 years, as recommended by the Kyoto protocol), whereas storage may be either short-term or long-term.

The large capacity of carbon exchanges with the atmosphere, the huge uncertainties about the response in soil carbon to global changes in climate and land use practices (positive or negative feedback) and lastly the fact that soil carbon is the only pool that humans can manage. All these factors show the crucial global interest of better understanding the fate of soil carbon.

1.2-The current representation of soil carbon dynamics in Land Surface Models is not entirely satisfactory Current representations of mechanistic models for soil carbon dynamics Several mechanistic models have been proposed to express the variation of soil carbon with depth. In these models, the representation of the physical processes at work is incomplete, however. Further research is therefore required to improve the simulation of soil carbon dynamics. The most widely used soil carbon models are included in the RothC Model [START_REF] Coleman | Simulating trends in soil organic carbon in long-term experiments using rothc-26.3[END_REF] which also incorporates a simple soil water model, and the CENTURY model system [START_REF] Parton | Analysis of factors controlling soil organic matter levels in great plains grasslands[END_REF] which includes more complex models for plant growth, management operations, etc.

The RothC model simulates profiles for the top 30cm of soil while the CENTURY model simulates profiles for the top 20cm [START_REF] Falloon | Modelling soil carbon dynamics[END_REF]. Most mechanistic models are calibrated for the top centimeters but there is a need to include deep soil carbon into models, particularly since this stable carbon can be reintegrated into the global carbon cycle thanks to changes in climate and land use practices. In addition, a study done by [START_REF] Todd-Brown | Causes of variation in soil carbon simulations from cmip5 earth system models and comparison with observations[END_REF] shows that the majority of Earth System Models (ESMs) cannot reproduce grid-scale variation in soil carbon and may be missing key processes. Differences across soil carbon models included in ESMs are primarily due to differences in the estimation of Net Primary Product (NPP) and the parametrization of soil decomposition sub-models. The weakness and the limitation of soil carbon dynamics models comes also from the fact that these models are parametrized under specific management and climatic conditions. Furthermore, ESMs seldom consider depth carbon and even when they do so, discretization does not consider the changes in physical conditions and superposition results of box model layers. This overview of mechanistic models for soil carbon dynamics points out the importance of considering the total carbon of soils and of extending conceptualizations of processes to all scales of time and space.

Newly revealed processes and deep carbon are missing In addition, there are large uncertainties about the processes that slow down mineralization and protect the organic matter in soil. Among these processes, one can distinguish: spatial inaccessibility to microorganisms and enzymes, hydrophobicity, encapsulation in organic macromolecules, litter resistance, organic matter-mineral interactions, etc. Until now, a major challenge has been to prioritize the role and impact of the stabilization process on soil carbon models [START_REF] Paul | The nature and dynamics of soil organic matter: Plant inputs, microbial transformations, and organic matter stabilization[END_REF]. It will be a great challenge to express the new concepts of soil carbon stabilization/destabilization by differential equations in order to incorporate them into mechanistic modeling. Furthermore, the majority of soil carbon mechanistic models underestimate the amount of soil carbon since deep carbon is not considered in the 1.3-A statistical approach is an alternative to better represent uncertainties on soil carbon dynamics In fact, integrating the new concepts of soil carbon protection into mechanistic models, considering what takes place in deep soil layers, and extending the conceptualization of processes to all scales of time and space, requires intensive development work. For that reason, a statistical approach based on observed empirical soil carbon data is useful in order to understand soil carbon dynamics, represent the different sources of uncertainties and provide answers about soil carbon behavior for the near future. The first attempt at the statistical modeling of soil carbon dynamics was made by Mathieu et al. (2015). Their proposed statistical approach was based on a world wide meta-analysis of 122 soil radiocarbon profiles collected from 87 articles in the soil science and archeology/paleoclimatology literature under different climatic and environmental conditions. The unknown model parameters were estimated relying on frequentist statistical inference. The study also sought to identify the effect of climatic and environmental factors on soil carbon dynamics. The analysis done by Mathieu et al. (2015) showed that the age of topsoil carbon was primarily affected by the climate and vegetation. In contrast, the results obtained on deep soil proved that the carbon content was impacted more by soil taxa than by climate. Moreover, they argued that the dependence on soil type points out the effect of other pedologic traits such as clay content and mineralogy. However, their interpretation was based on an expert analysis of the predictive results obtained without considering any statistical selection procedure to assess confidence about these judgements.

1.4-Contribution of my research work

The crucial aim of my research is to improve the statistical model proposed by Mathieu et al. (2015) in order to better express the soil carbon dynamics, using Bayesian inference for estimation. This inference has the advantage of taking into account the uncertainties on the unknowns and made it possible to integrate into the statistical model the knowledge on soil carbon dynamics given by soil science experts (see Appendix 7.2 for further details). A subsequent goal was to put into practice a Bayesian selection approach in order to assign a probabilistic judgment and numerically quantify the respective contributions of climatic and environmental factors such as: land use, soil type, temperature, precipitation, aridity index, etc. on soil carbon dynamics. A particular concern is to predict the gain or loss of soil carbon by computing the carbon stock and residence time when changes in temperature or land use occur. Moreover, it will be useful to know which type of land use conversion can sequestrate more soil carbon and predict the soil carbon response if the temperature increases by 1.5°C or 2°C. Finally, we propose a Gaussian Bayesian model that considers jointly the soil carbon content and radiocarbon activity. This model takes into account, on the one hand, the correlation between soil carbon content and radiocarbon, and on the other, hand the correlation between depth measurements. This model is constructed in such a way as to provide information on soil carbon at deep layers. Our study will be useful to have an overview of the behavior of soil carbon dynamics in a context of global warming and will help make some decisions concerning land use practices.

The statistical modeling of soil carbon dynamics has several important advantages: a better representation of uncertainties on soil carbon dynamics, the presence of various tools that numerically quantify quantities of interest for soil scientists, and faster responses to the issues of today and the near future.

The manuscript is organized as follows: in the first chapter we detail the soil carbon database used for the study.

We illustrate the heterogeneity of data sources, the available climatic and environmental information and the variety of sample sizes between different sites. The second chapter is an article published in the Société Française de Statistique SFdS journal in which we discuss the statistical model used on radiocarbon data. This chapter underlines the performance of Stochastic Search Variable Selection (SSVS) which is a Bayesian selection approach used as a first attempt to numerically quantify the climatic and environmental factors. The results obtained on artificial data

Introduction

show that SSVS can misrepresent some significant categorical explanatory covariates such as soil type and land use. However, a cross validation test on real radiocarbon data, conducted on the statistical model including all climatic and environmental factors and the SSVS model, showed that the latter achieves a better prediction and adjustment level. Chapter 3 is motivated by the results obtained on the SSVS approach. It gives an overview of three recent Bayesian selection methods appropriate for categorical potential predictors: Bayesian Sparse Group Selection, Bayesian Group Lasso based on spike and slab priors and Bayesian effect fusion using model-based clustering.

In this chapter, these three methods are applied on a simple regression model in order to better understand the functions and the characteristic features of each of the prior specifications. This chapter also includes a tutorial on these three Bayesian selection methods using Just Another Gibbs Sampler (JAGS) for Markov chain Monte

Chapter 2

• Vegetation is heterogeneously reported: from vague information (e.g. "natural vegetation"), to very detailed information (including precision at the species level). In most cases, the main type of ecosystem (forest, field, grassland, savanna, desert) is also indicated.

• Soil levels: they are reported as upper and lower depths (cm) of the sampling slice.

• The soil fraction on which 14 C and organic carbon content were measured. It includes bulk, density fraction, particle size fraction, molecular fraction, even the specific molecule.

• Soil horizon (i.e. L for litter, O horizon, A, B, etc.) if available.

• Based on the soil horizon designation or designated as such in the article, paleosols are also specified.

• Radiocarbon activity provided for different sampling levels. Various units are used: yr BP, pMC, F 14 C , ∆ 14 C (see Appendix 7.1). All values are reported as such and translated into F 14 C.

• Soil organic carbon content provided for different sampling levels, as organic carbon concentration and/or as stock depending on what is available in the article. Soil organic carbon concentration is given as %wt or g/kg or derived unit. Soil organic carbon stock is expressed as kg/ha and derived unit.

• Bulk density is seldom available and if so, is given either for different depths or as the mean value whatever the depth.

• Other information such as clay content, granulometry, pH, soil texture are also reported when available (a few cases only).

Processed database

In order to focus on soil organic carbon in general and not on the specific aspect of dynamics, we removed the sites with the following features:

1. Soil levels corresponding to a "paleosol" (244 levels from 51 profiles) were removed since they no longer have any carbon exchange with the atmosphere.

2. Levels above the horizon O (soil litter) are not considered.

3. Some studies were carried out on specific molecules or granulometric fractions that are not representative of all soil organic matter. We only kept data obtained on a "bulk", "bulk after HCl" and a "bulk after concentrated HCl". These three supports provide a more correct overall picture of the total soil organic matter (152 profiles were removed).

4. Sites with overlapping layers were removed (this concerned two sites).

5. Thirty-four profiles with less than three observations were removed to ensure a good estimation and prediction for statistical inference for both soil radiocarbon and soil carbon content dynamics. In fact, four and three unknown parameters have to be estimated for the 14 C and organic carbon statistical models, respectively. Three observations is thus the minimum required number.

6. Profiles with an unknown soil type, land use or vegetation cover were removed: this concerned 8 profiles with unknown soil type and 32 profiles with unknown land use type.

7. Six profiles showing odd patterns of organic carbon distribution were removed from the carbon content dynamics modeling (but kept for 14 C modeling).

Chapter 2

After data cleaning, only 131 profiles of soil radiocarbon and 125 profiles of soil carbon content remained for the statistical modeling.

Note on carbon content

To predict soil carbon profiles, we decided to work on soil carbon concentration (e.g. %wt) rather the soil carbon stock (e.g. g/m 2 ), since concentration is measured, while the soil carbon stock is calculated based on soil carbon concentration and on soil bulk density. However, bulk density is rarely or never provided in the articles used to build the database. Therefore, a pedotransfer function was used to predict the bulk density in the database, in order to complete the datasets. [START_REF] Alexander | Bulk densities of california soils in relation to other soil properties[END_REF] provided the most generic equation, where bulk density was derived from carbon concentration. But Alexander's equation is much too generic and does not account for soil type nor agro-pedo-climatic conditions, thus resulting in major uncertainties [START_REF] Tifafi | Large differences in global and regional total soil carbon stock estimates based on SoilGrids, HWSD, and NCSCD: Intercomparison and evaluation based on field data from Usa, England, Wales, and France[END_REF]. Hence, although carbon stock is more relevant for agronomical and climatic purposes, in view of the fact that it would greatly increase the uncertainties compared to carbon concentration, it was decided to establish the statistical model on the carbon concentration profiles. A second step will be to extend to carbon stock, from the modeled profile.

Potential explanatory covariates affecting soil carbon dynamics

The behavior of soil carbon was investigated by modeling the dynamics of soil radiocarbon and of carbon content.

Numerical (temperature, precipitation, etc.) and categorical (soil type, ecosystem type) predictors were considered for the meta-analysis. Explanatory covariates such as clay content, pH and granulometric information are not considered in this study since this information was seldom available.

The geographical information such as latitude, longitude and altitude are not considered in the statistical study since they do not impact the soil carbon dynamics. They are reflected in climatic parameters, such as temperature and to a lesser extent in "ecosystem".

Potential climatic numerical predictors

Climatic information is of prime importance to specify soil carbon dynamics. Statistically, taking all the monthly records of temperature and precipitation parameters (33 variables) into consideration would decrease the estimation and prediction performances of the linear model by increasing the variance of the estimated coefficients and making the model very sensitive to minor changes. In addition, it may enhance multicollinearity problems (Figure 2.2). For these reasons, in a first step, the number of predictors was reduced from 33 to 9. To select the potential climate predictors in this first step, we summarize information given by the monthly temperature and precipitation by considering:

• the extremes of temperature and precipitation regimes: minimum and maximum monthly precipitation (min_P and max_P, respectively), minimum monthly temperature (min_T and max_T, respectively),

• the mean annual temperature and precipitation (MAT and MAP, respectively),

• the seasonal shift between the warmest and coldest months (Dif_T),

• the seasonal shift between the wettest and the driest months (Dif_P),

• the aridity index (AI).

Chapter 2

• Regosol-Arenosol, Arenosol and Leptosol were grouped together, as they share a lack of significant soil horizon.

• Kastozem and Phaeozem were pooled into Chernozem, as they are all humus-rich soils, at least for their surface layers.

• Fluvisol was merged with Cambisol, as they are often found in conjunction.

• Nitisol was merged with Ferralsol as they both contain a high amount of iron oxides that interact with organic compounds.

• Lastly, Plinthosol and Planosol were grouped with Gleysol as they are all susceptible to waterlogging and drought/frost stress. Table 2.1: Merging of WRB soil type groups for soil radiocarbon profiles according to expert advice. For ease of reference, we will hereafter use the soil type group "short name" (e.g. Chernozem) to refer to the concatenation of the merged groups (e.g. Chernozem, Kastanozem, Phaeozem).

We grouped "land use" and "vegetation" into a single term, "ecosystem" that combines the two types of information. We identified 9 categories that we further merged into 6 groups (Table 2.2) : field, forest (forest, natural-forest), cultivated-forest, natural-grassland, cultivated-grassland and undefined natural (natural + naturaldesert + natural-savanna) (see Table 2.2,Figure 2.4). The aggregation of ecosystem type was done in order to include the anthropogenic impact. In order to avoid categories with a small number of observations and to increase the prediction power of the statistical model, we created a group called "others" (Table 2

Chapter 2

Köppen published his first climate classification in 1884 [START_REF] Köppen | Die wärmezonen der Erde, nach der Dauer der heissen, gemässigten und kalten Zeit und nach der Wirkung der Wärme auf die organische Welt betrachtet[END_REF] and later improved it (e.g. [START_REF] Koppen | Das geographische system der klimat[END_REF][START_REF] Geiger | Klassifikation der klimate nach w. köppen. Landolt-Börnstein-Zahlenwerte und Funktionen aus Physik[END_REF]) to achieve the Köppen-Geiger climate classification. It has since been regularly updated.

Kottek provided the latest version in 2006 [START_REF] Kottek | World map of the köppen-geiger climate classification updated[END_REF] for the second part of the 20 th century. The classification is based on the division of earth climates into 5 main climate groups: equatorial (Figure 2.5: acronym beginning with A, red-derived zones), arid (B, yellowish), warm temperate (C, green), snow (D, blue), polar (E, gray). The main groups are then subdivided into subgroups (29 in total). The definition of climatic groups is provided in [START_REF] Kottek | World map of the köppen-geiger climate classification updated[END_REF], table 1 and the first subdivision is reproduced as follows (Table 2.5

T max ≤ 10 • C 2 ET Tundra climate 0°C ≤ T max ≤ 10 • C 2 EF Frost climate T max ≤ 0 • C 0 Table 2
.5: Description of the Köppen-Geiger classification (1st and 2nd letter description only) and number of radiocarbon profiles selected from the database that correspond to the different subgroups (last column). P min (P max ) and T min (T max ) are for the minimum (maximum) monthly precipitation and temperature, P ann is for the MAP, S and W subscripts are for summer and winter respectively. P th = 2 * MAT + a, with a = 0 if at least 2/3 of MAP occurs in winter, a = 28 if at least 2/3 of MAP occurs in summer and a = 14 otherwise. The calculation key implies that the polar climates (E) have to be determined first, followed by the arid climates (B) and subsequent differentiations into the equatorial climates (A) and the warm temperate and snow climates (C) and (D), respectively.

Examination of the database shows that 20 of the selected profiles belong to "equatorial climates", 1 to "arid climates", 70 to "warm temperate climates", 38 to "snow climates" and 2 to "polar climates (Table 2.5). At the first order, this results in a homogeneous representativeness of intermediate climate types, i.e. tropical, warm temperate and snow climates, leaving out extreme climates. "Arid climate" is represented by only one profile from the Sonora Desert, AZ, USA) and "polar climate" by two Italian mountain profiles. Warm temperate climates are overrepresented, and this tendency is even stronger when compared with the land surface ratio they occupy (Figure 2.5). This is due to the fact that most agronomical studies have traditionally been performed in temperate regions, while investigations in other regions are a recent phenomenon. At the second order, however, not all sub climates are present in the database. Some sub-climates are overrepresented. So, whereas equatorial climates are well balanced between "monsoon" and "savannah climates", "rainforest climate" is not represented. The high weight of "fully humid warm temperate climates" is in line with the respective weight of the "dry season" and "fully humid" within this type of climate. The imbalance is rather between Cs and Cw where one would have expected an Chapter 2 tiveness with most of the profiles regularly distributed in the [0.15;2] range, i.e. from very arid to humid zones, and including 7 profiles above 2.

Temperature seasonal shifts vary from a few • C to more than 35 • C (Figure 2.6, panel b), covering the wide range of climates from tropical to continental. The profile distribution within this range is less continuous than for annual parameters, and rather patchy. Nevertheless, the [2; 38 • C] range of seasonal shift is homogeneously swept.

Briefly, our database will allow the statistical model to predict soil profiles under equatorial, warm temperate and snow climates but precludes using the model in extreme conditions such as arid climates, polar climates and very wet regions.

Soil type diversity and distribution of profiles

The database offers a good representativeness of the land coverage diversity. As shown in table 6, the 17 soil types from the database represent about 71% of the land coverage. However, the distribution between the 17 categories (or between the 9 merged groups of soil types) is unbalanced. So with 12 profiles Andosol, which only represents 1% of total land area on Earth will be much better constrained than areno-regosol represented by only 7 profiles but covering 22% of land area. The difference in representativeness is not that large for the other soil types.

Introduction

A significant current issue when trying to predict our planet's future is to understand the feedback effects between climate evolution and the future soil carbon balance. Soil constitutes the largest carbon pool in interaction with atmospheric carbon, containing 2000 to 2400 Gt of organic carbon in the first meter, i.e. at least the equivalent of 250 years of current fossil carbon emissions that are estimated at 10±0.6 Gt/year [START_REF] Stocker | Climate change 2013: the physical science basis: Working Group I contribution to the Fifth assessment report of the Intergovernmental Panel on Climate Change[END_REF].

The stock of soil organic matter (SOM) has been defined as a balance between input of organic matter through vegetation and loss through microbial decomposition. A large variation in the soil organic carbon (SOC) stock amongst soil types and land use has been shown, ranging from 2 kg/m 2 for arenosols to more than 10 kg/m 2 for podzols (Batjes, 1996). Regarding land use, [START_REF] Martin | Spatial distribution of soil organic carbon stocks in france[END_REF] show that relationships between soil organic carbon stocks and pedo-climate depend on the type of land use and that they differ between forest and cultivated soil.

The global analyses carried out by [START_REF] Carvalhais | Global covariation of carbon turnover times with climate in terrestrial ecosystems[END_REF] and [START_REF] He | Radiocarbon constraints imply reduced carbon uptake by soils during the 21st century[END_REF] point out the lack of knowledge of carbon residence time in soil and an increasing concern about the importance of climate factors in the variability of carbon storage. For instance, a temperature increase may clearly impact the activity of soil microorganisms and the subsequent organic carbon sequestration by soils. Moreover SOM evolution plays a key role in the CO 2 atmospheric content since the soil is a crucial pool for CO 2 emission or sequestration. No consensus has been reached, however, on the relative importance of the various climatic factors that affect SOM dynamics, such as temperature, precipitation, aridity, moisture, etc.

In fact, several questions remain unclear for soil scientists: Could soil capacity be durably increased to sequestrate more carbon by changing land use? What quantitative changes in SOM occur when modifying agricultural practices? Will that change the soil carbon stock/the organic matter residence time? What is the contribution of each climatic or environmental factor to soil carbon? Is the potential increment of the soil carbon stock to be considered as sustainable ? These questions highlight the importance of assessing the uncertainties as well as understanding the complex mechanisms of soil carbon dynamics. To investigate this point through data collection, in addition to soil carbon concentration, F 14 C measurements are also taken into account to describe SOM dynamics on the grounds that radiocarbon content can be considered as a clock that registers SOC residence time [START_REF] Scharpenseel | Radiocarbon dating of soils-problems, troubles, hopes. Paleopedology: Origin, Nature and Dating of Paleosols[END_REF].

A worldwide meta-analysis of radiocarbon profiles is described in Mathieu et al. (2015). In their study, a hierarchical non linear model is designed under the frequentist paradigm with inference performed by the "Expectation-Maximization" algorithm. The radiocarbon dynamics is parameterized as a smooth function of depth with random effects taking into account potentially explanatory climatic and environmental factors. Once calibrated, the model is used for statistical prediction along various typical scenarios of (modified) forcing conditions; according to an expert interpretation of their predictive results, deep soil carbon dynamics is driven more by soil type than by climate. Although such a result was based on a statistical model with unknown parameters, there was no direct probabilistic judgment to assess the strength of their claim. Our aim in this article is to scrutinize this claim more closely and check the robustness of the statistical model in view of the many uncertainties: how confident can we be in the effective roles of environmental covariates and climatic factors for the phenomenon under study? What are the respective contributions of signal and noise in what we see? In this paper, we revisit Mathieu's approach under the Bayesian paradigm since Bayesian inference has the advantage of expressing the uncertainties on the unknowns throughout the statistical analysis. We re-parametrize the model to obtain more directly interpretable parameters, change the error term structure to clarify the different sources of uncertainties, and weight the influence of the climatic and environmental drivers for prediction. A Bayesian selection approach is hereby used in order to quantify the contribution of climatic and environmental factors to soil carbon dynamics. Several Bayesian selection approaches for linear models have been developed in the literature such as: Variable Selection for Regression Models (VSRM) (Kuo and Mallick, 1998), Gibbs Variable Selection (GVS) (Dellaportas and Ntzoufras, 1997) and Stochastic Search Variable Selection (SSVS) (George and McCulloch, 1993). These methods were applied within the framework of the linear model, where y i is the outcome response for individual i (i = 1, . . . , n) predicted by p potential explanatory covariates x i j for j = 1, . . . , p. The intercept is expressed by α and the measurement error by e i .

y i = α + p ∑ j=1 θ j x i, j + e i e i ∼ N(0, σ 2 ),
with N(µ, σ 2 ) referring to the Normal distribution with mean µ and variance σ 2 . In frequentist selection methods, each variable combination corresponds to a different model, so the variable selection chooses among all possible models the best sub-model based on criteria for model selection such as: AIC, BIC and Mallows's C p . For a large number of covariates p, it is not computationally achievable to consider all 2 p possible sub-models. The idea of Bayesian variable selection is to define a binary variable I j which indicates whether a covariate x j is influential (I j = 0) or not influential (I j = 0) for the response y. I j is generated from a Bernoulli prior. The VSRM and GVS selection methods set θ j = I j × β j . For VSRM, I j and β j are considered as independent and β j is sampled from a vague normal prior (Kuo and Mallick, 1998). For GVS, β j is sampled from a conditional prior that depends on I j such as a Gaussian mixture prior: P(β j |I j ) = (1 -I j )N(µ, S 2 ) + I j N(0, τ 2 ), where µ, S 2 and τ 2 are hyperparameters chosen to ensure good mixing of the Monte Carlo Markov Chains (MCMC) (Dellaportas and Ntzoufras, 1997). Therefore, these two Bayesian selection methods enable the best sub-model to be selected by affecting null regression coefficients (I j = 0 ⇒ θ j = 0) for the non influential predictors. SSVS considers a "slab and spike" prior which depends on I j for the regression coefficients β j , with a spike around 0, and a flat slab elsewhere. Then if I j is null, we assign a value close to 0 for θ j , which means that the corresponding covariate x j has no effect on response y. This method was chosen for the present study. The major difference between the scope of the original a soil layer characterized by the depths of its top and bottom levels. Such a preliminary data cleaning was based on the following criteria: i-) the radiocarbon data must have been acquired on bulk organic carbon (not on specific fraction, nor specific molecule), ii-) sites must contain more than 3 observations. Figure 1 shows the site locations where radiocarbon data at various depths were collected. The number of observations varied from one site to another (from 3 to 88 measurements per site). For each of the 159 profiles, the following information of interest is provided: sampling year, location, climate, soil type, land use, organic carbon content and radiocarbon. Soil texture is not considered as it is poorly recorded in many articles from the literature. More details on the database can be found in (Mathieu et al., 2015). In this study, the potential climatic and environmental explanatory covariates are as follows:

• Mean annual precipitation (MAP), mean annual temperature (MAT), aridity index (AI), and absolute shift between July and January temperatures (∆T) are included as representative of the average climate and seasonality of the site. The aridity index, defined by UNEP as the ratio of annual precipitation to annual potential evaporation, was obtained from the FAO 10-minute mean climate grids for global land areas for the period 1950-2000(Trabacco and Zomer, 2009).

• Latitude (Lat).

• The atmospheric radiocarbon of the sampling year ( 14 Catm).

• Soil type with 13 different categories ordered alphabetically: andosol, arenosol, cambisol, chernozem, ferralsol, fluvisol, gleysol, kastanozem, luvisol, nitisol, phaeozem, podzol, vertisol. We pooled phaeozem and kastanozem soil types into chernozem due to similar characteristics, as they are poorly present in the database. Hereafter, soil type will be considered as a categorical variable with 11 levels.

• Vegetation and land use were combined to form a new factor dubbed "ecosystem", with originally 9 categories distinguished as follows: cultivated-field, cultivated-forest, cultivatedgrassland, forest, natural, natural-desert, natural-forest, natural-grassland and natural-savanna. We pooled natural-desert into the "natural" ecosystem. Ecosystem will therefore be considered as a categorical variable with 8 levels. Among the 159 profiles collected, 55 with missing climatic or environmental covariates were removed from the database. After previous data cleaning, the dataset finally includes 104 sites TABLE 1. Contingency table of pairwise combinations of levels between soil type and ecosystem. Abbreviation "C" in column names refers to Cultivated and "N" to Natural. and 951 records. The dataset results from an observational study, which may lead to some confusion due to the spurious association between the correlated and/or the poorly contrasted covariates. The very small number of observations for pairwise-combinations of factors (even a null number for many of them) rules out the possibility of including interactions between soil type and ecosystem in the model (see Table 1). In addition, we anticipate a poor precision of the estimates of the effects of categorical covariates since their design matrix, shown by Table 1, is unfortunately very strongly unbalanced.

C-Field C-Forest C-Grassland Forest Natural N-forest N-Grassland N-Savanna Total Andosol 0 2 1 0 1 4 0 0 8 Arenosol 0 2 0 0 1 0 0 1 4 Cambisol 2 0 0 1 0 4 2 0 9 Chernozem 2 0 0 0 0 0 11 0 13 Ferralsol 0 0 0 1 0 9 1 2 13 Fluvisol 2 0 0 2 0 0 0 0 4 Gleysol 2 1 0 0 1 0 0 0 4
To illustrate the composition of the dataset, the boxplots in Figure 2 show the average F 14 C variation versus the mean levels of non overlapping soil layers, for the most frequent types of profiles collected. This figure only shows average profiles for some specific combinations and prevents any strict interpretation as the number of observations differs from top to depth, and as soil horizon width differs from one profile to another (we do not expect the intensity of processes to be the same at the same depth between two profiles). Figure 2 shows as expected that the radiocarbon decreases with depth: with higher input, topsoil OM is more rapidly renewed (and thus shows a younger age) than deep soil OM.

A multivariate hierarchical non linear model

The statistical model structure that mimics (eqs 1 and 2) variations of F 14 C with depth along a profile within a given site is similar to the one considered in Mathieu et al. (2015). It differs only in the homogeneous variance for the measurement error and in the unit chosen to report radiocarbon concentration. Let S = 104 be the total number of carbon soil profiles under study. We note m s the number of measurements available for site s. Therefore, for each site s ∈ {1 : S} and each depth x ∈ {1 : m s }, the F 14 C content experimental record y(s, x) is modeled by:

y(s, x) = g φ (s), x + ε(s, x), ε(s, x) ∼ N(0, σ 2 )
(1)

g φ (s), x = φ 1 (s) + (φ 2 (s) -φ 1 (s))exp - x φ 3 (s) φ 4 (s) (2)
As indicated in Fig 3, the structure of the previous statistical model is interpreted:

φ 1 represents F 14 C in deep soil, φ 2 refers to the topsoil F 14 C, φ 3 is related to the depth at half maximum of the F 14 C peak, φ 4 describes the more or less rapid decrease of F 14 C.

The ε terms represent the within-site discrepancies between the observed and the adjusted F 14 C profiles.

To express the variability between the different sites, a linear link is considered between each of β i = (β i1 , . . . , β iP ) ′ ∈ R P , where i = 1, 2, 3, 4, represents the fixed covariate effect relative to each latent variable, and E i ∈ R P the corresponding centered random effect. φ i and E i are defined as the following vectors:

φ i = (φ i (1), φ i (2), . . . φ i (S)) ′ and E i = (E i (1), E i (2), . . . , E i (S)) ′ .
In this case study, the number of columns P in the design matrix X is equal to 23 (P = 1 + (11 -1) + (8 -1) + 6). In fact, "1 + (11-1) + (8-1)" is the dimension of the two-way explanatory subspace spanned by the categorical factors "Soil type " and "Ecosystem" that includes the constant. 6 is the number of quantitative regressors. The quantitative regressors in X are normalized to allow comparison of their effects in a rescaled unit. Due to the presence of dummy variables generated by the two categorical factors, the number of columns of the design matrix ( 23) is greater than the number of explanatory covariates (6+2).

Bayesian selection model: The variable selection procedure is expected to reveal the most influential explanatory variables for the assemblage of the four latent sub-models given with 2 categorical covariates and 6 quantitative ones by equation 3. The idea is to consider a "slab and spike" prior [START_REF] Dellaportas | Bayesian variable selection using the Gibbs sampler[END_REF] for each β i parameter, with a spike centered at 0, and a flat slab elsewhere. Each combination of included variables corresponds to a different model, so variable selection amounts to choosing among all possible 2 P sub-models if the model considered were a simple linear model with P regressors. For a large number of covariates P , it would be therefore not feasible to consider each possible model separately. In our case, it may seem at first glance that P = 8, leading to only 2 8 = 256 sub-models for each of latent model given by Eqs.3 and 4. Hence the idea of a Bayesian variable selection, where we consider a stochastic exploration of this immense combinatorial set of possible models [START_REF] O'hara | A review of bayesian variable selection methods: what, how and which[END_REF]. In this article, we concentrate on the Stochastic Search Variable Selection introduced by George and McCulloch (1993). This approach is applied to the latent layers φ 1 , φ 2 , φ 3 and φ 4 , in presence of categorical covariates. For the selection procedure, we need to define an indicator variable I i j where i = 1,2,3,4 and j = 1, . . . , P as follows:

I i j =
1 if variable X j has an effect on φ i 0 otherwise

The mixture prior for β i j depends on I i j :

P(β i j |I i j ) = (1 -I i j )N(0, τ 2 i j ) + I i j N(0, c 2 i j τ 2 i j ) (5)
where i = 1, 2, 3, 4 and j = 1, . . . , P. Based on this Gaussian mixture, τ i j must be small, in order to sample β i j around 0 in situations when variable X j is not influential, but not strictly restricted to zero, though, otherwise the Gibbs sampler will rarely be able to flip from I i j = 0 to visit I i j = 1. Furthermore, c i j must be large enough for β j to be given a flat prior when X j is needed in the model. A semi-automatic approach to selecting τ i j and c i j was proposed by George and McCulloch (1993) considering the interaction point and relative heights at 0 of the marginal densities. They recommended "good" choices for the couple (σ β i j /τ i j , c i j ), where σ β i j is the observed standard error associated with the least squares estimate βi j . However, a more appropriate prior for β suggested later is the hyper-g prior proposed by [START_REF] Liang | Mixtures of g priors for bayesian variable selection[END_REF] based on the g-prior introduced by [START_REF] Zellner | On assessing prior distributions and bayesian regression analysis with g-prior distributions. Bayesian inference and decision techniques[END_REF]. This extension of the g-prior has been widely studied and widely used in a regression context. The specification of g is mostly based on a model selection criterion such as the Akaike Information Criterion (AIC, see [START_REF] Burnham | Aic model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons[END_REF]), the Bayesian information criterion (BIC, see [START_REF] Bhat | On the derivation of the bayesian information criterion[END_REF], the Deviance Information Criterion (DIC, see [START_REF] Spiegelhalter | Bayesian measures of model complexity and fit[END_REF], etc. Here, the β prior can be understood as a mixture of spike and slab of g-priors. In order to specify g and to ensure a reasonable order of magnitude for β , the hierarchical model without the selection step is first adjusted with a hyper-g prior (with a vague uniform prior at the upper level of the hierarchy). The value of g will be fixed as the posterior mean of this preliminary estimation and used afterwards for the Bayesian selection approach. In that respect, when I i, j is equal to 1, β i, j will be generated from the following g-prior N(0, g i σ 2 i (X ′ X) -1 j, j ), to be considered as the slab prior. In contrast, according to the concept of the spike prior, which should be more centered at 0, the β i, j corresponding to I i, j = 0, will be generated from a g-prior, where the variance is much smaller N(0, (1/c) * g i σ 2

i (X ′ X) -1 j, j ). The hyperparameter c is specified by the user based on a model comparison with different values of c according to the previously cited selection model criteria or to a cross validation study. A hyper prior can also be proposed for c (uniform prior). The model for Bayesian selection of variables can be finally summed up as follows:

• Likelihood: for each site s ∈ {1 : S} and each depth x ∈ {1 : m s }:

y(s, x) ∼ N(g(φ (s), x), σ 2 ) with φ (s) = (φ 1 (s), φ 2 (s), φ 3 (s), φ 4 (s))
• Latent variables:

φ i ∼ N S (Xβ i , σ 2 i I) i = 1, 2 log(φ i ) ∼ N S (Xβ i , σ 2 i I) i = 3, 4 with φ i = (φ 1,i , . . . , φ s,i , . . . , φ S,i ), φ i ∈ R P .
• Priors:

• 1/σ 2 ∼ G(0.001, 0.001) • 1/σ 2 i ∼ G(0.001, 0.001) for i = 1, 2, 3 and 4 G( , ) refers to the gamma distribution.

• An intercept is always included and common across all sub-models, for j = 1,2,3,4

β j1 ∼ N(0, 10000) • for quantitative covariates j = 2, . . . , K • β 1 j |I 1 j ∼ (1 -I 1 j ) * N(0, g 1 σ 2 1 (X ′ X) -1 j, j c 1 ) + I 1 j * N(0, g 1 σ 2 1 (X ′ X) -1 j, j ) • β 2 j |I 2 j ∼ (1 -I 2 j ) * N(0, g 2 σ 2 2 (X ′ X) -1 j, j c 2 ) + I 2 j * N(0, g 2 σ 2 2 (X ′ X) -1 j, j ) • β 3 j |I 3 j ∼ (1 -I 3 j ) * N(0, g 3 σ 2 3 (X ′ X) -1 j, j c 3 ) + I 3 j * N(0, g 3 σ 2 3 (X ′ X) -1 j, j ) • β 4 j |I 4 j ∼ (1 -I 4 j ) * N(0, g 4 σ 2 4 (X ′ X) -1 j, j c 4 ) + I 4 j * N(0, g 4 σ 2 4 (X ′ X) -1 j, j )
For j = 2, . . . , K and i = 1, 2, 3, 4:

I i j ∼ B(p i j = p) with B(.) the Bernoulli distribution (6)
i.e. all models are a priori equiprobable.

• For the categorical covariates numbered j = K + 1, . . . , P, with covariate C j having n j levels, the algorithm ensures that the n j modalities are either taken or dropped all together during Monte Carlo Markov Chain (MCMC) iteration:

• for each level s = 1, . . . , n j :

• β 1s |I C j ,1 ∼ (1 -I C j ,1 ) * N(0, g 1 σ 2 1 (X ′ X) -1 j, j c 1 ) + I C j ,1 * N(0, g 1 σ 2 1 (X ′ X) -1 j, j ) • β 2s |I C j ,2 ∼ (1 -I C j ,2 ) * N(0, g 2 σ 2 2 (X ′ X) -1 j, j c 2 ) + I C j ,2 * N(0, g 2 σ 2 2 (X ′ X) -1 j, j ) • β 3s |I C j ,3 ∼ (1 -I C j ,3 ) * N(0, g 3 σ 2 3 (X ′ X) -1 j, j c 3 ) + I C j ,3 * N(0, g 3 σ 2 3 (X ′ X) -1 j, j ) • β 4s |I C j ,4 ∼ (1 -I C j ,4 ) * N(0, g 4 σ 2 4 (X ′ X) -1 j, j c 4 ) + I C j ,4 * N(0, g 4 σ 2 4 (X ′ X) -1 j, j )
For j = k + 1, . . . , P and i = 1, 2, 3, 4:

I C j ,i ∼ B(p C j ,i = p)
All levels of a categorical factor therefore receive the same prior selection probability, but more informative priors can be designed, if prior expertise is available to tune the respective importance of the explanatory variables.

The calculation of the posterior distributions of the parameters is based on MCMC algorithms such as the Metropolis-Hastings and Gibbs Sampler [START_REF] Dellaportas | Bayesian variable selection using the Gibbs sampler[END_REF]. The SSVS is easily implemented in JAGS (Just Another Gibbs Sampler), as exemplified in Ntzoufras et al. (2002, pp.13-17).

Results and Discussion

Performing SSVS on artificial data

In this section, we illustrate the performance of SSVS on latent layers for artificial data generated according to the non linear multivariate statistical structure model ( 1)+( 2)+( 5)+(6) when:

1. all independent covariates are quantitative;

2. all covariates are quantitative, and some of them are correlated;

3. the covariates are mixed: some are quantitative and the others are categorical.

The purpose of this artificial data generation is to understand and study the challenges in the application of SSVS when the selection aims at hidden sub-models and the model structure is more complex than a simple univariate regression.

SSVS on latent layer models with independent quantitative covariates:

• Example 1: The artificial dataset mimics the real one by taking the same number of sites (104 sites) and depth measurements (951 records). In this example, 6 quantitative (continuous) predictors are considered. The predictors are generated as independent standard normal vectors, X 1 , . . . , X 6 iid N 104 (0, 1), so that they are practically uncorrelated. The regression effects are set to β 1 = (0, 1, 0, 1, 0, 1), β 2 = (0, 0, 1, 1, 0, 0), β 3 = (1, 0.8, 0, 0.7, 0, 1) and β 4 = (1, 0, 0, 1, 0.8, 0.8) with standard deviations σ 1 = σ 2 = σ 3 = σ 4 = 0.1 and σ = 0.1. The intercept is equal to 1 and will always be kept in the proposals of the latent layer models.

SSVS on latent layers with correlated quantitative covariates:

• As shown in Fig 8, for the real case, covariates may be correlated. Example 2 is designed to illustrate how SSVS reacts in the presence of high collinearity. The only difference with example 1 is that the matrix design X contains 2 correlated explanatory variables. X 5 and X 6 are defined as follows:

X 5 = 2 × X 3 X 6 = X 2 + 1.5 × Z, Z ∼ N(0, 1)
SSVS on latent layers with mixed covariates:

• Example 3 introduces categorical variables: this time, the latent linear models φ 1 , φ 2 , φ 3 and φ 4 contain 6 quantitative (X 1 , . . . , X 6 ) covariates and 2 qualitative factors (F 1 and F 2 ) with respectively 8 and 11 levels. Contrast-sum coding was considered to remain coherent with the presence of quantitative covariates. Regression effects were set to β 1 = (1, 0, 0, 1, 0, 1, 0, 1), β 2 = (0, 1, 0, 0, 1, 1, 0, 0), β 3 = (1, 1, 1, 0.8, 0, 0.7, 0, 1) and β 4 = (0, 0, 1, 0, 0, 1, 0.8, 0.8). 0 and 1 are the index vectors of length 7 or 10 with 0 and 1's corresponding to categorical covariates (position 2 and 3 of the regression coefficients vector).

The first position in β 1 , β 2 , β 3 and β 4 is always equal to 1 and refers to the intercept. σ i , i =1,2,3,4 and σ are fixed as in Example 1. Similar to real data, the experimental design of artificial data is strongly unbalanced.

Sensitivity analysis of the prior for SSVS latent layers on independent quantitative covariates

In order to suggest reasonable values of g 1 , g 2 , g 3 and g 4 for the spike and slab g-priors on the regression effect parameters, the inference of the linear model with a vague uniform prior (g i ∼ U(10, 1000), i = 1,2,3,4) on g was run. The posterior means of g 1 , g 2 , g 3 and g 4 were plugged into the SSVS model. The prior inclusion probability was fixed to 0.5 in the paper of George and McCulloch (1993). This choice is common for Bayesian selection models since it ensures for all explanatory covariates the same probability of being included in the model. Yet, this prior is informative and favors sub-models with half of the covariates included. For the purpose of studying the impact of the inclusion probability p on the selection results, the SSVS was tested under three different prior specifications:

Bayesian selection model applied to soil carbon dynamics 141

Sensitivity analysis prior for SSVS latent layers on independent quantitative covariates

In this section, we test the "best" choice of the hyperparameter c for the β prior specification.

We consider the following values of c: 10, 100, 1000 and 5000. The MCMC is run for 30, 000 iterations after a burn-in of 10, 000 iterations. In addition, a Beta prior B(2, 2) is proposed for the inclusion probability p. The four panels in Table 4 show, for Example 1 of artificial data, the SSVS performance under different priors on β 1 , β 2 , β 3 and β 4 . These tables show the three most frequent models with the false inclusion (False positive) or exclusion (False negative) rates of predictors.

For the different spike and slab priors, SSVS performs extremely well for c i = 10, 100 and 1000 ( i = 1, 2, 3, 4) since the best sub-models identified for each of the four latent layers contain TABLE 5. Comparison between the three SSVS models with different values of c according to the DIC criterion. The best model is identified by the lowest DIC estimation.

c DIC no false detections (see the first line of the panels a), b), c) and d)). The best sub-models do not contain any false choice. As expected, as the value of c increases, the posterior distribution becomes more peaked, which can be explained by the increase in probability appearance along these settings. In fact, the probability of the most visited model increases with higher values of c (see the probability values in the first row of the previous four tables). For example, in , the best sub-model under c = 10 is visited 870 times throughout 30, 000 iterations, while the best sub-model under c = 5000 is visited 1770 times. The SSVS with c = 1000 is identified as the best according to the DIC estimations. Moreover, a vague uniform prior can be proposed on parameter c in order to have a better estimation. Generally speaking, SSVS performs well on latent layer models with independent quantitative covariates.

3.1.3. The presence of collinearity increases false detection on SSVS in the latent layer George and McCulloch (1993) showed that collinearity may reduce the efficiency of SSVS by increasing the number of promising models in a linear model framework. Collinearity between some covariates in a latent layer model can also increase the rate of false positives/negatives especially when one of the correlated covariates is influential but the other is not. The SSVS model is now considered with a Beta prior on the probability selection p (p ∼ B(2, 2)) and a vague uniform prior on c (U(5, 1000)).

Figure 5 illustrates how correlated covariates restrict SSVS performances. The SSVS model provides a probability judgment about the most frequent explanatory covariates combination. In addition to that, the SSVS also provides a probability judgment about the inclusion of each of the explanatory covariates on the different sub-models identified throughout MCMC iterations. Here, the Posterior Inclusion Probabilities (PIP) for each covariate separately are illustrated in Fig. 5. In the first and third panels, the selected covariates correctly specify the influential covariates taken a priori into account to generate artificial data. Outputting, both X 3 and X 5 as non influential, and X 2 and X 6 as influential for φ 1 was expected since the correlated covariates were a priori both influential/not influential at the same time. With regard to the second panel, φ 2 was generated taking into account X 3 , while X 5 is omitted a priori. Therefore as X 5 is correlated with X 3 , SSVS misleads and selects X 5 . Likewise, X 2 and X 3 were not taken into account when generating φ 4 . As a result, two false choices are reported, the exclusion of X 5 and the inclusion of X 2 .

SSVS performance within latent layer mixed covariates (quantitative and qualitative)

The algorithm for mixed covariates was developed to give the same inclusion probability to all levels of the same categorical covariate. The results obtained in Example 3 highlight some limitations of SSVS with regards to the presence of categorical covariates in the latent layer. It can be clearly seen that SSVS may fail to detect some influential explanatory categorical covariates. However, SSVS does not seem to induce false choice inclusion. In our case study, it considers a categorical covariate as influential only if it is actually influential: it can miss some of them but does not induce false positives. The new dummy covariates needed to handle the presence of categorical covariates F 1 (8 levels) and F 2 (11 levels) strongly increase the dimensions of the space of competing models to be as false negative detections respectively for φ 3 and φ 4 linear models. These results highlight a limitation of SSVS related to the presence of categorical covariates in latent layers. It is clear that SSVS fails to detect some influential explanatory categorical covariates. However, SSVS does not induce false choice inclusion in this case study. In other words, it considers a categorical covariate as influential only if it is actually influential. Such avoidance of false choice inclusion might stem from the fact that SSVS with even prior weights tends to dampen the selection probability of a categorical covariate with a big number of modalities. In fact, the prior distribution of β k ∈ R M when covariate k is selected (i.e.

I k = 1) is proportional to 1 (gσ k (X ′ k X k ) -1 ) M .
Consequently, when M becomes large, the prior distribution P(β k |I k = 1) will vanish to 0. For that reason, SSVS may seem to be reluctant to select a categorical covariate with a high number of levels.

Variance sensitivity analysis for SSVS

As mentioned above, George and McCulloch (1993) designed and applied SSVS to detect explanatory covariates directly linked to the observed response whereas we applied it to covariates buried in latent layers in the framework of a hierarchical Bayesian model. To complete the assessment in our specific case, we evaluated the sensitivity of SSVS to the variance within the latent layer. Overall, sensitivity variance analyses highlight that an increase in variability between sites (expressed by the σ 1 , σ 2 , σ 3 and σ 4 of the latent layer models) does decrease SSVS robustness to select the best subset of covariates.

In our specific case, two sources of variability are to be distinguished: variability between sites expressed by σ 1 , σ 2 , σ 3 and σ 4 and variability within the same site expressed by σ . In order to test SSVS sensitivity to intersite variability changes, we simplified the proposed statistical model by fixing φ 2 , φ 3 and φ 4 . SSVS was applied only on φ 1 , which has a linear effect on the F 14 C response. We tested SSVS for four different values of σ 1 = (0.01, 0.1, 2.5, 3). Figure 6 shows the posterior inclusion probability for one of the considered covariates "X 2 ", for different σ 1 settings. Figure 6 clearly illustrates the impact of σ 1 on the posterior inclusion probability (PIP): the more σ 1 increases, the more PIP decreases. It even reaches a PIP close to 0.5 for σ 1 = 3, leading to a potential false choice (exclusion) of an important variable.

SSVS on observed radiocarbon profiles

Application of SSVS on soil F 14 C profiles

The aim of this section is to highlight the contribution of SSVS to understanding which climatic and environmental factors are likely to control soil carbon dynamics. Based on the results obtained on artificial data, it can be claimed that the presence of categorical covariates in the model can produce false exclusions of some of the influential categorical covariates. In addition, the correlation between some covariates such as temperature and latitude, may yield false detection, especially if they do not have the same effect on latent layers as we showed in subsection 3.1.3. visited sub-models for φ 1 and φ 2 . φ 1 and φ 2 are respectively related to the deep and topsoil F 14 C. All explanatory covariates are selected for φ 3 in its latent model. All the categorical covariates (land use or soil type) selected with a probability higher than 0.5 are included in the best submodel. So, land use is very surely included in the best sub-models of φ 1 , φ 3 , φ 4 and soil type in the φ 3 best sub-models. In contrast, every categorical covariate not selected (PIP smaller than 0.5), may be significant for the model since the SSVS approach can yield negative false detection for categorical covariates. For example, soil type is a priori not included in the best sub-model of φ 1 but might still be significant to explain deep soil radiocarbon. Moreover, a posterior probabilistic beliefs on the association of explanatory covariates is provided by looking at the most frequent covariate combinations throughout the MCMC iterations (see Table 9). According to the Table 9, the frequency visits to the best sub-models are very small with respect to the total number of iterations (180,000) and maybe not all the sub-models are explored by the MCMC. Moreover, the full models are detected as the best sub-models for three of the latent layers φ 2 , φ 3 and φ 4 . However, the covariates Posterior Inclusion Probabilities (PIP) highlight that the best model chosen should contain the covariates with a PIP higher than 0.5. Furthermore, for more detailed investigations, the unknown parameters of the statistical radiocarbon model are re-estimated, taking into account all the covariates for which the PIP is higher than 0.5 (see Fig. 7). In addition, as the SSVS may miss the inclusion of some influential categorical covariates, one may wonder whether the soil type has really no effect on the φ 1 latent linear linear model or whether it is perhaps simply not detected by the SSVS model. The answer to this question is reported in the following table.

Results of 2 most frequent combinations of covariates identified by Stochastic Search Variable Selection

Comparison of DIC for 5 sub-models taking into account for some sub-models the drawback of SSVS when categorical covariates are present in the model TABLE 10. Model* contains the explanatory covariates with a PIP higher than 0.5. To investigate whether a non selected categorical covariate is significant, we add respectively to Model*, the non included categorical covariates (land use or ,soil type) identified with a PIP smaller than 0.5. The Table displays the DIC criteria comparisons between the different models.

Models DIC

Most frequent model (denoted Model1 for each of latent layers in Table 9) -1703 Model* = the model adjusted on the covariates where their PIP are higher than 0.5 (see Fig. 7) -1837

Model* + considering the soil type for φ 1 -1897 Model* + considering the soil type for φ 1 and φ 2 -1890 Model* + considering the land use for φ 2 and soil type for φ 1 -1968 Model* + considering the soil type for φ 1 and land use for φ 2 and φ 4 -1879

The DIC comparison in Table 10, shows that the best model is the one that includes both PIP> 0.5 detected explanatory covariates, i.e. "soil type" for deep soil radiocarbon (φ 1 ) and "land use" for topsoil radiocarbon (φ 2 ) (DIC = -1968). In addition, this result highlights that the SSVS is misleading in that it detects two significant categorical covariates (2 false negatives). The final selection of covariates for the radiocarbon model is summed up in Table 11.

Selection results for the best sub-model: the climatic and environmental factors that affect soil radiocarbon dynamics TABLE 11. The final selected covariates for each of the four latent layer models φ 1 , φ 2 , φ 3 and φ 4 . For the third latent layer φ 3 all explanatory covariates are selected. Furthermore, for φ 1 and φ 4 four covariates are identified among 8 as significant while 5 covariates are detected for φ 2 as influential towards the 8 potential climatic and environmental factors. A further point is the correlation among covariates. For example, temperature and seasonal shift are positively correlated (see Fig. 8). This could suggest that temperature may not be really influential for φ 1 as its inclusion may be the result of its correlation with the highly influential covariate "seasonal shift". However, if we take a look at the second panel of Fig 7, we can see that seasonal shift has an effect on φ 2 , which is not the case for temperature, indicating that the correlation between temperature and seasonal shift does not seem to affect SSVS performance that much.

Best

soil. An increase in microbial activity that leads to higher mineralization will result in a weaker weight of older components relative to newly input ones in the age distribution of the mixture of soil components within the same soil layer (panel c in Fig. 12). This will result in an increase of radiocarbon in deep soil but a decrease in topsoil radiocarbon where the weight of the peak-bomb derived components decreases due to a higher mineralization (panel c in Fig. 12). We face the opposite effect in the case of processes that will enhance the organic matter stabilization and will better preserve old material (panel d in Fig. 12). Keeping in mind that point, our results for the deep soil highlight a positive posterior effect of mean annual temperature and a negative posterior effect of seasonal shift. In practice, an increase of 1°C in the mean annual temperature will result in an increase of radiocarbon of 0.12 and an increase of 1°C between the highest and the lowest monthly temperature will result in a decrease of radiocarbon by 0.03. This increase of deep soil radiocarbon with temperature is in agreement with a higher mineralization associated to an enhancement of microbial activity under higher temperature. Likewise the decrease of radiocarbon with seasonality matches what is known about the impact of seasonality on soil dynamics with much younger soils, i.e. with a higher turnover under the tropics than in boreal, i.e. continental areas, where soil shows a much lower turnover and thus yields much lower radiocarbon. Topsoil is negatively impacted by atmospheric radiocarbon, seasonal shift and mean annual temperature and positively impacted by aridity. Most of the profiles included in the database were sampled posteriorly to the 1960s, i.e. for years during the bomb peak decrease with an overrepresentation of the 1990s. The bomb peak gradually penetrates into soil layers with a time lag that depends of the mean residence time of components in the different layers. With a mean residence time of 100 yrs, the maximum of F 14 C will be in the early 2000's. Thus, the negative impact of the atmospheric F 14 C reflects the fact that an increase in the atmospheric F 14 C means that sampling was made some years before, when the bomb peak had not yet reached its maximum in soil. The dilution effect of bomb-peak derived components is thus higher, yielding a lower (closer to 1) mean radiocarbon. However, this effect remains very low (-0.01 decrease of topsoil radiocarbon associated to an increase of atmospheric radiocarbon by 1) reflecting the dilution effect of the bomb-peak and the disequilibrium of the database in which sites sampled in the 1990s are overrepresented. Negative impacts of seasonal shift and mean annual temperature by -0.02 and -0.01 respectively are the counterpart for topsoil of what is observed for deep soil. An higher mineralization for the mean annual temperature, leading for an higher loss of bomb-derived organic matter and a better preservation for seasonal shift yielding for a relative gain of the oldest components. It is noteworthy that impacts for topsoil appear much smaller than for deep soil. This result is counter-intuitive and no reason for that can be advanced. The positive impact of aridity is in agreement with a well-known low microbial efficiency in arid environments compared to humid ones. An increase in aridity results in a better preservation of the bomb-peak derived components and thus to an increase in the topsoil radiocarbon. The effect of aridity remains very low at +0.01. A large difference exists between the magnitudes of the posterior estimation of the influential covariates of the latent variable for topsoil and deep soil. While an explanation stemming from the database disequilibrium can be put forward to explain the low magnitude of atmospheric radiocarbon, no clear evidence can be provided for the other covariates.

Extensions and challenges

Database: To better predict the evolution of soil carbon dynamics with climate change and land use change practices, there is a need to collect more data for the type of soil (arenosol, fluvisol and gleysol) and ecosystem (natural/savanna, cultivated/grassland and forest) about which we do not have much information. In this study, the experimental design was strongly unbalanced, which affects the precision when estimating the quantities of interest: φ 1 , φ 2 , φ 3 and φ 4 . Furthermore, optimization of the experimental design should take into account the type of contrast used to solve the redundancy of the model caused by the presence of categorical explanatory covariates. An interesting new track will be to know where to take new samples and for which climatic and environmental conditions in order to improve the overall estimation. Another issue associated with data is correlation. Some of the explanatory covariates are naturally correlated (see Figure 8). For example, the aridity index (AI) is proportional to the mean annual precipitation (MAP) by definition (see eq. 4) since:

AI = MAP ET p ET p : potential evapotranspiration rates
SSVS is sensitive to the presence of correlated covariates as already seen in Section 3.1.3 (see Fig 8). More investigation can be done considering other Bayesian predictive criteria for model selection according to the paper by [START_REF] Piironen | Comparison of bayesian predictive methods for model selection[END_REF].

Improving the Bayesian selection model. The test carried out on artificial data shows that SSVS does not always detect influential categorical explanatory covariates. This issue could be thoroughly explored using the Bayesian effect fusion approach introduced by Pauger and Wagner (2017). They proposed a Bayesian approach for a sparse representation of the effect of a categorical predictor in linear models. The originality of their work is that it not only allows selection of categorical covariates but also induces fusion among the categorical covariate levels which have essentially the same effect on the response. Besides this approach, Bayesian variable selection for group Lasso presented in the paper by [START_REF] Xu | Bayesian variable selection and estimation for group Lasso[END_REF] selects variables both at the group level and also within a group. Revisiting the traditional Bayesian approach to the group Lasso problem, they developed a Bayesian group Lasso model with spike and slab priors for problems that also require selection of categorical explanatory variables.

Conclusion

In this paper, we have discussed the performance and limitations of SSVS on latent layers in the framework of a hierarchical Bayesian model applied to soil radiocarbon. The results on artificial data show that collinearity may lead to false inclusion or exclusion in the best sub-model selected.

Besides collinearity, if variability on the latent model response is high, the posterior inclusion probability may blur the effect of influential explanatory covariates as exemplified in Section 3.1.5. Furthermore, SSVS is not always able to select the influential categorical covariates, but at least does not seem to consider a covariate as influential unless it is indeed the case. Despite the complexity of SSVS compared to the full model, we show that the Bayesian selection approach has a better adjustment and prediction level in our case study. Finally, the application of SSVS to soil F 14 C profiles highlighted the influence of soil types on soil carbon dynamics by impacting deep soil F 14 C, topsoil F 14 C and F 14 C incorporation. Our results also indicate that temperature affects deep soil F 14 C more than topsoil.

Note: In figure 5 of the article, the vector β 3 at the top of the third panel should be β 3 = (1, 0.8, 0, 0.7, 0, 1) instead of β 3 = (0, 0.8, 0, 0.7, 0, 1)

Rectified references:

• Burnham, K. P., Anderson, D. R., and Huyvaert, K. P. ( 2011). Aic model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behavioral Ecology and Sociobiology, 65(1):23-35.

• Martin, W., Smith3, M., Jolivet, B., and Arrouays (2011). Spatial distribution of soil organic carbons stocks in France.

Chapter 5

In Chapter 3, the performance of the Stochastic Search Variable Selection (SSVS), originally built to select the significant numerical predictors, was adapted to the hierarchical non linear model with latent variables proposed for modeling soil radiocarbon dynamics. The SSVS had to be modified to handle the presence of categorical predictors such as the soil type and the ecosystem type of the sampled profile. In addition, the results obtained on artificial data generated according to the proposed statistical model highlighted that the detection of some significant categorical predictors can be misleading with SSVS. Further investigation done on the real data indicated that the prediction and model fitting were better after the inclusion of some categorical predictors detected as non significant by the SSVS. The results of the SSVS on the soil radiocarbon dynamics, published in the SFDS journal (see Chapter 3), gave rise to new questions: How can Bayesian selection methods handle the presence of categorical predictors?

What are the Bayesian selection methods appropriate for categorical predictors that already exist in the literature?

Which Bayesian selection methods to choose?

These questions led to the submission of the second article "Bayesian selection approaches for categorical predictors using JAGS" from which Chapter 4 derived. Three Bayesian Selection approaches appropriate for categorical predictors: Bayesian Group Lasso with Spike and Slab, Bayesian Sparse Group Selection and Bayesian Effect Fusion were tested on a simple linear model with categorical predictors in order to highlight the prior specifications for each Bayesian Selection method rather than the complexity of the statistical model.

In this chapter, we now test the performance of the Stochastic Search Variable Selection, introduced in Chapter 3, for the numerical predictors (mean annual temperature, aridity index, etc.) and the three Bayesian selection approaches explored in Chapter 4 for the categorical predictors (soil type and ecosystem type) on the latent layers within the framework of a non linear hierarchical model. Here, it should be pointed out that the framework of the application becomes more complex. This chapter will be organized as follows: Section 5.1 recaps the structure of the statistical model proposed for the soil radiocarbon dynamics and the climatic and environmental factors considered for this statistical analysis. It also introduces the structure of the statistical model for soil carbon content modeling and the potential climatic and environmental factors. Section 5.2 will explore the Bayesian Full model, in which all the explanatory predictors are included, and the implementation of Bayesian Selection approaches previously introduced in the manuscript. This section will be divided into two parts: the first part applies these methods to the soil radiocarbon dynamics and the second one will address the application of these methods to the soil carbon content dynamics. Following this outline, all the Bayesian methods are then compared with respect to both soil radiocarbon and soil carbon content dynamics based on the Bayesian selection criteria for model comparisons.

After choosing the best subset of climatic and environmental factors for soil radiocarbon and soil carbon content dynamics, Section 5.3 will touch on the physical interpretation of the selected climatic and environmental factors. This section will also comment on the expectations of soil scientists versus the selected and the non selected climatic and environmental predictors. For the Bayesian inference, with no prior information about the precision of the parameters, the scale and the hyperparameter shape of the Gamma distribution must be chosen so as to give a very disperse prior. The most widely used parametrization of the Gamma distribution is to have the same number in both hyperparameter positions. Thus, by choosing 0.001 for both hyperparameters, the precision has a mean 1 and a large variance of 1000.

For the regression effects for latent linear models, we assumed Zellner's g priors. In fact, Zellner's g prior [START_REF] Zellner | On assessing prior distributions and bayesian regression analysis with g-prior distributions. Bayesian inference and decision techniques[END_REF]) is based on the idea that the regression effect estimation should be invariant to changes in the scale of the regressors. Some linear algebra shows that this condition is satisfied if the mean and the variance of the Normal distribution prior on the regression effect β are equal to 0 and k(X ′ X) -1 respectively. A popular specification is to set k = gσ 2 for positive values of g. The choice of g can be based on many popular model selection criteria, such as the Akaike information criterion (AIC), the Bayesian information criterion (BIC) and others. However, assuming a prior on g has the advantage of avoiding paradoxes such as "Bartlett's Paradox" and the "Information Paradox".

Briefly, "Bartlett's Paradox" states that the null model would always be preferred to any other model when g→ ∞. On the other hand, when the coefficient of determination R 2 → 1, the Bayes factor converges to a constant instead of going to ∞ as we consider that the datasets fit the model better. This problem is called the "Information Paradox".

For the Full Bayesian model, we consider a vague Uniform prior on g, assigning the same weight to all possible values of g.

FULL BAYESIAN MODEL FOR SOIL RADIOCARBON MODELING

• Likelihood: for each site s ∈ {1 : S} (S = 131) and each depth x ∈ {1 : m s }, the likelihood is written as:

-F 14 C(s, x) ∼ N + g(φ (s), x), σ 2 (truncated Normal distribution) or log F 14 C(s, x) ∼ N log(g(φ (s), x)), σ 2 (logarithm transformation)

where, g(φ (s),

x) = φ 1 (s) + (φ 2 (s) -φ 1 (s)) * exp -x φ 3 (s) φ 4 (s)
• Latent variables: for each site s, the linear latent models are defined as:

φ i (s) ∼ N + X[s, ] * β i , σ 2 i i = 1, 2, 3 log(φ 4 (s)) ∼ N X[s, ] * β 4 , σ 2 4
• Priors:

-1/σ 2 ∼ Gamma(0.001, 0.001) for i = 1, 2, 3, 4, P = 20 and S = 131:

-1/σ 2 i ∼ Gamma(0.001, 0.001) -β i ∼ N(0, g i σ 2

i (X ′ X) -1 ) where 0 ∈ R P and X ∈ M S,P (R) g i ∼ Uni f orm (5,5000) Chapter 5

The test quantity T (y, θ ) used to compute the p-value of the Posterior Predictive Check (P.P.C) is equal to the mean of the difference between the real radiocarbon observation and the non linear mean obtained according to the latent variables φ1 , φ2 , φ3 and φ4 sampled from their posterior distributions. A better agreement between the model and the dataset was achieved under the additive model (0.58 p-value closer to 0.5 and farther from 1).

Based on the Deviance Information Criterion and the p-value of the Posterior Predictive Check criterion, an additive model with a truncated Normal distribution on F 14 C was chosen for soil radiocarbon modeling.

Bayesian

Group Lasso with Spike and Slab prior

BGL-SS model specification and choice of hyperparameters

As stated in Chapter 4, the Bayesian Group Lasso with Spike and Slab prior (BGL-SS) is the simplest Bayesian selection approach used for both categorical and numerical predictors. Furthermore, this Bayesian selection method requires few hyperparameters to tune. Within the framework of the non linear hierarchical model, the BGL-SS was been applied to the latent linear models for φ 1 , φ 2 , φ 3 and φ 4 . Thus, the BGL-SS model is specified as follows:

BAYESIAN GROUP LASSO WITH SPIKE AND SLAB PRIOR • Likelihood: for each site s ∈ {1 : S} (S = 131) and each depth x ∈ {1 : m s }, the likelihood is written as:

-F 14 C(s, x) ∼ N + g(φ (s), x), σ 2 where, g(φ (s), x) = φ 1 (s) + (φ 2 (s) -φ 1 (s)) * exp -x φ 3 (s) φ 4 (s)
• Latent variables: for each site s, the linear models are defined us:

φ i (s) ∼ N + X[s, ] * β i , σ 2 i i = 1, 2, 3 log(φ 4 (s)) ∼ N X[s, ] * β 4 , σ 2 4 • Priors:
-1/σ 2 ∼ Gamma(0.001, 0.001) -for the fourth latent variable φ 4 , we propose a vague Normal prior:

β 4 ∼ N P (0, 100 * I) where 0 ∈ R P , I ∈ M P,P (R) and P = 20 -Intercepts:

β 0,i ∼ N(0, 1) for i = 1, 2 β 0,3 ∼ N(0, 1000)
for i = 1, 2, 3, 4: 1/σ 2 i ∼ Gamma(0.001, 0.001)

Chapter 5

For BGL-SS, there are two feasible criteria for selecting the best subset of explanatory predictors:

1. The posterior median estimation of regression effects: [START_REF] Xu | Bayesian variable selection and estimation for group Lasso[END_REF] showed that the sub-model selected according to the posterior median estimation of regression effects has an excellent performance for both variable selection and estimation (Table 5.4).

Latent variables Physical interpretation

Best subset of predictors

φ 1 deep radiocarbon MAT φ 2
topsoil radiocarbon Land, Soil, Dif_T, MAP φ 3 radiocarbon incorporation Land, soil, Dif_T, MAP, AI, MAT Table 5.4: The best subsets of climatic and environmental predictors for latent linear models of φ 1 , φ 2 and φ 3 chosen according to the posterior median estimation of regression effects.

2. The Posterior Inclusion Probability (P.I.P): [START_REF] Barbieri | Optimal predictive model selection[END_REF] showed that, for a linear model, the optimal predictive model is often the median probability model, which is defined as the model consisting of predictors which have overall posterior probabilities greater than or equal to 1/2 of being in a model (Table 5.5). Table 5.5: The best subsets of climatic and environmental predictors for latent layers φ 1 , φ 2 and φ 3 chosen according to the Posterior inclusion Probability (PIP). The significant predictors are detected with a PIP at least equal to 0.5. The predictors highlighted in blue were the ones detected in addition to those identified by the posterior median estimation of regression effects. The predictors are ordered according to the PIP.

Latent variables Physical interpretation

According to Table 5.5, the Aridity Index (AI) and the minimum precipitation (min_P) were included, in the φ 2 and φ 3 latent linear models respectively, in addition to the predictors identified for φ 2 and φ 3 in Table 5.4. The inclusion of the Aridity Index with a rather small PIP (PIP = 53) can be explained by the correlation between this predictor and the Mean Annual Precipitation (MAP) detected as significant (Figure 2.3 in Chapter 2) with both selection criteria (Pearson correlation, P.C(AI,MAT) = 0.66). The inclusion of the minimum precipitation (min_P) with a PIP equal to 55 can also be explained by the positive relationship existing between this predictor and the Aridity Index (P.C(min_P,AI) = 0.67).

The model selected according to the posterior median estimation of regression effects will be more robust to the collinearity problem than the sub-model selected according to the posterior inclusion probabilities of predictors.

Chapter 5

Comparison of "Best BGL-SS sub-models" Two sub-models are in competition, the one obtained using median probability criteria (Table 5.4) and the one using PIP criteria (Table 5.5). Which one to choose? In light of the Bayesian model checking and the comparison of the model criteria presented in the previous section, the most parsimonious model was chosen. The DIC was computed by testing the hierarchical non linear model (without the selection step) and considering only the predictors detected as significant with each of the BGL-SS selection criteria. The aim of this step is to readjust the estimation of regression effects by removing the predictors detected as non significant.

According to the DIC criterion, the model that best fits the dataset is the posterior median model including predictors with Posterior inclusion Probabilities (PIP) at least equal to 0.5 (the lowest recorded DIC is -2356). The Cross Validation (C.V), obtained by splitting the data into (k = 5) groups, showed no important difference between the Relative Error (R.E) computed for the sub-model selected according to the posterior median estimation of the regression effects and the Relative Error (R.E) for the model selected based on the posterior inclusion probabilities for predictors. We chose the PIP model since it has a better adjustment of data (Table 5 For Bayesian Group Lasso with Spike and Slab prior, we selected the sub-model according to the Posterior Inclusion probability in Table 5.5, as it shows a better prediction power than the model selected based on the posterior median estimation of regression effects.

Bayesian Sparse Group Selection

BSGS model specification and choice of hyperparameters

After identifying that the soil type influences the topsoil radiocarbon dynamics by applying the Bayesian Group Lasso with Spike and Slab prior, soil scientists may wonder:

Which are the levels with significant contributions (non null effects) within the soil type that affect the topsoil radiocarbon?

The Bayesian Sparse Group Selection is constructed to handle the presence of categorical predictors in the linear model and to answer this type of question. This method performs better than the Sparse Group Lasso with Spike and Slab prior in terms of selecting the active predictors as well as identifying the active levels within the selected predictors (Chapter 4). The Bayesian Sparse Group Selection model will be applied on the latent linear models within the framework of the hierarchical non linear model proposed for soil radiocarbon dynamics. For the prior specifications, the hyperparameters are chosen as follows:

• Prior inclusion probability: we propose a Beta distribution, with both hyperparameters equal to 2, on prior probability for predictors and levels inclusion. This induces a distribution for the number of included variables which has a heavier tail than the binomial distribution, allowing the model to learn about the degree of sparsity.

Chapter 5

• The variance hyperparameter of Normal distribution for regression effects: Zellner's g prior is proposed on the regression effect when the relative predictor is selected. The g * i , for i = 1, 2, 3, 4 represents the posterior mean estimation of g i obtained from the Bayesian Full model (including all predictors). This choice was made to allow plausible values for regression effects.

BAYESIAN SPARSE GROUP SELECTION

• Likelihood: for each site s ∈ {1 : S} (S = 131) and each depth x ∈ {1 : m s }, the likelihood is written as:

-F 14 C(s, x) ∼ N + g(φ (s), x), σ 2 where, g(φ (s), x) = φ 1 (s) + (φ 2 (s) -φ 1 (s)) * exp -x φ 3 (s) φ 4 (s)
• Latent variables: for each site s, the linear models are defined us:

φ i (s) ∼ N + X[s, ] * β i , σ 2 i i = 1, 2, 3 log(φ 4 (s)) ∼ N X[s, ] * β 4 , σ 2 4 • Priors:
-1/σ 2 ∼ Gamma(0.001, 0.001) for i = 1, 2, 3, 4:

-1/σ 2 i ∼ Gamma(0.001, 0.001)

for the fourth latent variable φ 4 : β 4, j ∼ N(0, g * 4 σ 2 4 (X ′ j X j ) -1 ) for j = 1,. . . , P and P = 20 -Intercepts, for i= 1, 2, 3:

β 0,i ∼ N(0, g * i σ 2 i (X ′ 0 X 0 ) -1 )
for the two categorical predictors g = 1, 2 and for the latent layer i = 1, 2, 3: * binary indicator for categorical predictor inclusion: ν g (i) ∼ Ber(p predictor (i))

* binary indicator for level inclusion:

λ lg (i)|ν g (i) ∼ (1 -ν g (i))δ 0 + ν g (i)Ber(p level (i))
* predictor prior inclusion probability:

p predictor (i) ∼ Beta(2, 2)
* level prior inclusion probability:

p level (i) ∼ Beta(2, 2)
Chapter 5

• Priors:

-regression effect of level l within the categorical predictor g:

β lg (i)|λ lg (i)ν g (i) ∼ (1 -λ lg (i))δ 0 + λ lg (i)ν g (i)N(0, g * i σ 2 i (X ′ g X g ) -1 )
• for the six numerical covariates, we used the Bayesian selection method of Kuo and Mallick (1998) . This approach is based on the Stochastic Search Variable Selection introduced by George and McCulloch (1993). The spike and slab prior proposed on regression effects is replaced by a mixture model between a mass point at 0 (Dirac distribution) and a Normal distribution.

-regression effect for numerical predictor n = 1. . . ,6:

β n (i) ∼ (1 -ν n (i)) * δ 0 + ν n (i)N(0, g * i σ 2 i (X ′ n X n ) -1
) -binary indicator for numerical predictor n:

ν n (i) ∼ Ber(p predictor (i))

Selection results of Bayesian Sparse Group Selection

Best subset of predictors selected

The best sub-model chosen according to the Posterior Inclusion Probability is summarized in Table 5.7: The best subsets of climatic and environmental predictors for latent linear models for φ 1 , φ 2 and φ 3 chosen according to the Posterior inclusion Probability (PIP). The significant predictors are detected with a PIP at least equal to 0.5.

Identification of the active levels within the selected categorical predictors

The active levels are detected with a posterior inclusion probability at least equal to 0.5. In Figure 5.7, the results of the Posterior Inclusion probability (PIP) obtained within the influential categorical predictors are presented as bars versus the real variation for deep and topsoil radiocarbon according to ecosystem and soil type shown as boxplots.

Chapter 5

The BEF model for the soil radiocarbon dynamics is written as follows:

BAYESIAN EFFECT FUSION MODEL-BASED CLUSTERING

• Likelihood:

for each site s ∈ {1 : S} (S = 131) and each depth x ∈ {1 : m s }, the likelihood is written as:

-F 14 C(s, x) ∼ N + g(φ (s), x), σ 2
where, g(φ (s),

x) = φ 1 (s) + (φ 2 (s) -φ 1 (s)) * exp -x φ 3 (s) φ 4 (s)
• Latent variables:

for each site s, the linear models are defined us:

φ i (s) ∼ N + X[s, ] * β i , σ 2 i i = 1, 2, 3 log(φ 4 (s)) ∼ N X[s, ] * β 4 , σ 2 4
• Priors:

• 1/σ 2 ∼ Gamma(0.001, 0.001) for i = 1, 2, 3, 4:

• 1/σ 2 i ∼ Gamma(0.001, 0.001)

• for the fourth latent variable φ 4 :

β 4, j ∼ N(0, g * 4 σ 2 4 (X ′ j X j ) -1
) for j = 1,. . . , P and P = 20

• Intercepts, for i= 1, 2, 3:

β 0,i ∼ N(0, g * i σ 2 i (X ′ 0 X 0 ) -1 )
• for the latent variable φ i i = 1, 2, 3 and the categorical covariate X g g = 1, 2 with C 1 and C 2 levels respectively:

-

β gl (i) = ∑ C g l=1 ν l (i)N(µ l (i), ψ g (i)) + ν 0g (i)δ 0 -ν l (i) ∼ Dirichlet C g +1 (e 0 ) where e 0 ∈ R C g +1 * ν 0g (i) = 0 * ν l (i) ∼ N(m g (i), M g (i)) * ψ g i = V g (i)/k where V g (i) = 1 C g -1 ∑ C g l=1 ( βgl (i) -βg (i)) 2 and βg (i) = 1 C g ∑ C g l=1 βgl (i)
Chapter 5

• k and e 0 highlighted in blue need to be tuned by the user. m g (i) and M g (i) are specified according to the suggestions of [START_REF] Malsiner-Walli | Effect fusion using model-based clustering[END_REF].

• Like the BSGS, the selection of numerical predictors is based on Kuo and Mallick (1998) approach.

For i = 1, 2 , 3 and n = 1, . . . , 6:

-regression effect for the numerical predictor n:

β n (i) ∼ (1 -ν n (i)) * δ 0 + ν n (i)N(0, g * i σ 2 i (X ′ n X n ) -1
) -binary indicator for the numerical predictor n:

ν n (i) ∼ Ber(p predictor (i))

Selection results of Bayesian Effect Fusion model-based clustering

Results of the sensitivity analysis of the variance parameters

The Bayesian Effect Fusion was tested with three different values of k for the Gaussian mixture distribution variances: k = 10, 50 and 100. Two levels within the same categorical predictor are fused if their Posterior Fusion Probability (PFP) is at least equal to 0.5. Furthermore, a given level is fused to the baseline if the estimation of its

Posterior Median Regression Effect (PMRE) estimation is exactly equal to 0. The best sub-models are identified based on both Posterior Fusion effect and Posterior Median Regression Effect.

For these three values, one sub-model is identified (Table 5.8). The lowest DIC was recorded for the BEF with k equal to 50. However, the difference in the DIC between the BEF with k = 10 (DIC = -2354) and k = 50 (DIC = -2363) is slight. A difficulty of convergence is detected for the BEF with k = 100. Even when increasing the number of iterations, there are still some parameters that do not converge according to the potential scale reduction factor defined by [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF]. Results of fusion of levels within the significant categorical predictors

Results of fusion of levels within the ecosystem type

The clustering of the levels of the ecosystem type identified according to the PFP and the PMRE are represented as pie charts versus the observed variation of the radiocarbon F 14 C according to the ecosystem type illustrated by boxplots (Figure 5.8). As for soil radiocarbon modeling, there are also constraints for soil carbon content:

1. The soil carbon content reported in g/kg is always positive. Thus, using a Normal truncated distribution or a logarithm transformation of the soil carbon content is recommended.

2. The latent variables ω 1 , ω 2 and ω 3 should have positive values in accordance with their physical interpretations.

For the Bayesian inference, we consider a vague Gamma prior on the precision parameters 1/σ 2 c and 1/sd 2 i for i = 1, 2, 3, where both hyperparameters are equal to 0.001. With regard to regression effects, we propose a Zellner's g prior under a vague Uniform prior on g. The Full Bayesian model for soil carbon content dynamics is written as follows:

FULL BAYESIAN MODEL FOR SOIL CARBON CONTENT DYNAMICS

• Likelihood: for each site s ∈ {1 : S} (S = 125) and each depth x ∈ {1 : m s }, the likelihood is written as:

-C(s, x) ∼ N + f (ω(s), x), σ 2 c (truncated normal distribution) or -log (C(s, x)) ∼ N log( f (ω(s), x)), σ 2 c (logarithm transformation) where, f (φ (s), x) = ω 1 (s) + (ω 2 (s) -ω 1 (s)) * exp -x ω 3 (s)
• Latent variables:

for each site s, the linear models are defined as:

ω i (s) ∼ N + X * [s, ] * β i , sd 2 i i = 1, 2, 3
• Priors:

-1/σ 2 ∼ Gamma(0.001, 0.001) for i = 1, 2, 3, 4, P' = 19 and S = 125:

-1/sd 2 i ∼ Gamma(0.001, 0.001)

-β i ∼ N(0, g i sd 2 i (X * ′ X * ) -1
) where 0 ∈ R P ′ and X ∈ M S,P ′ (R)

g i ∼ Uni f orm(5, 10000) (vague prior)

Results of the Full Bayesian carbon content dynamics model

In order to decide how to model the response of soil carbon contents, the additive model with truncated Normal distribution was tested against the multiplicative model with the log transformation. The Deviance Information Criterion (DIC) and the Posterior Predictive Check (P.P.C) for both models are given in Table 5.10:

Chapter 5 5.2.2.2 Bayesian Group Lasso with Spike and Slab prior

BGL-SS model specification and choice of hyperparameters

For the Bayesian Group Lasso with Spike and Slab (BGL-SS), we followed the same concept of prior specification of the BGL-SS used for the soil radiocarbon dynamics. Thus, a Beta prior is proposed for the prior inclusion probability with both hyperparameters equal to 2.

BGL-SS selection results

Best sub-models selection criteria

Two Bayesian selection criteria can be used to choose the best sub-predictors for soil carbon content dynamics.

We recall that the first criterion is based on choosing the sub-model according to the posterior median regression effect estimation while the second criterion consists in including in the model all the predictors with a Posterior Inclusion probability (PIP) at least equal to 0. Table 5.11: The best sub-predictors detected according to the two Bayesian selection criteria used for the Bayesian Group Lasso with Spike and Slab prior. The sub-predictors selected according to the posterior median estimation of regression effects are indicated in the column called "Median model". The sub predictors detected according to the Posterior inclusion Probability (PIP) are represented in the column called "PIP model". This column also contains the PIP for the significant predictors. The symbol C in the second column refers to the carbon content.

Readjustment and best sub-models comparison

In order to better estimate the regression effects, the best sub-models based on the posterior median estimation of regression effects and the posterior inclusion probability were readjusted by removing the redundant predictors.

The Deviance Information Criterion, for both sub-models, is given in Table 5.12.

Best sub-models DIC Median model 2072 PIP model 2078 Table 5.12: The Deviance Information Criterion (DIC) for the two sub-models detected by the Bayesian Group Lasso with Spike and Slab. The Median model is based on the posterior median estimation of regression effects while the PIP model is based on the Posterior Inclusion Probability (PIP) for predictors.

According to Table 5.12, the best model chosen with respect to the lowest DIC is the Median model (DIC = 2072). Furthermore, the two sub-models can be considered as a step of Stepwise regression. In fact, the only

Chapter 5

The non detection of atmospheric radiocarbon as an influential predictor, for the topsoil radiocarbon, is due to a poor representation of the potential sampling years in the database, resulting in an overrepresentation of profiles collected in the 1990s.

Mean Annual Temperature (MAT)

Reminder Mean Annual Temperature (MAT) was selected as an influential predictor for deep F 14 C (φ 1 ) and F 14 C incorporation depth (φ 3 ) of the F 14 C profile and for the topsoil (ω 2 ) of the carbon content profile

It is globally assumed that a cooler temperature is associated with slower decomposition and increases the mean residence time of soil carbon. The smallest change in the soil carbon content may have a large impact on the concentration of CO 2 in the atmosphere [START_REF] Trumbore | Rapid exchange between soil carbon and atmospheric carbon dioxide driven by temperature change[END_REF]. The modalities of temperature impact on the decomposition rate of soil organic matter remain an interesting topic of discussion. A study done by [START_REF] Giardina | Evidence that decomposition rates of organic carbon in mineral soil do not vary with temperature[END_REF] suggested that the recalcitrant carbon is not sensitive to temperature variation. In contrast, [START_REF] Fierer | Litter quality and the temperature sensitivity of decomposition[END_REF] suggested that the non-labile organic matter is more sensitive to temperature than the labile pool. [START_REF] Fang | Similar response of labile and resistant soil organic matter pools to changes in temperature[END_REF], on incubated soils under changing temperature, found similar results. Likewise, [START_REF] Lefevre | Higher temperature sensitivity for stable than for labile soil organic carbon-evidence from incubations of long-term bare fallow soils[END_REF] highlighted, on long-term (to 79 years) bare fallow experiments, a strong relationship between the residence time of carbon organic matter and the temperature sensitivity of its mineralization: the more stable the organic matter, the more sensitive to temperature the organic matter will be. And finally, [START_REF] Conen | Warming mobilises young and old soil carbon equally[END_REF] pointed out that recalcitrant and labile pools have a similar temperature sensitivity.

A recent study done by [START_REF] Yan | The temperature sensitivity of soil organic carbon decomposition is greater in subsoil than in topsoil during laboratory incubation[END_REF] also took a position in this debate, pointing out different fates between top and deep soil. This study was based on a sequential temperature (8 • C to 28 • C) changing method applied on cultivated fields in China. Results showed that the average SOC decomposition rate was 59% to 282% higher in the topsoil than in the subsoil layer. In contrast, the temperature sensitivity values in the topsoil layer were significantly lower than those in the subsoil layer.

The increasing temperature in humid climates increases both plant growth and decomposition of soil organic matter. However, the relative increase in the decomposition rate of organic matter remains greater than the net primary production [START_REF] Oades | The retention of organic matter in soils[END_REF]. In a laboratory study, [START_REF] Hagerty | Accelerated microbial turnover but constant growth efficiency with warming in soil[END_REF] showed that the microbial turnover accelerates with temperature while the growth efficiency is not sensitive to temperature changes.

It is certain that in order to extract the sole effect of temperature on soil carbon content and dynamics, it would be necessary to work with several temperatures on the same environmental conditions (soil type, ecosystem, vegetation cover, etc.) which was not always applicable for the studies carried out. Several parameters influencing soil carbon are thus mixed and it is difficult to deconvolve the signal. Furthermore, the sensitivity of soil carbon content to temperature is often studied by confronting soil carbon results with Mean Annual temperature (MAT) or the seasonal shift of temperature not immediate temperature [START_REF] Smith | Impact of global warming on soil organic carbon[END_REF].

In this debate, the statistical results obtained on our meta-analysis of radiocarbon and soil carbon content can provide elements for decision making Chapter 5

Aridity Index (AI)

Reminder Aridity Index (AI) was selected as an influential predictor for topsoil (φ 2 ) and incorporation depth (φ 3 ) of the F 14 C profile and for none of the latent variables of the carbon content profile.

Aridity influences soil carbon inputs by affecting the production of above (loss of plant cover) and below ground plant biomass, water infiltration, microbial biomass and mineralization processes and thus the biogeochemical cycle of nutrients [START_REF] Ren | Responses of soil total microbial biomass and community compositions to rainfall reductions[END_REF]. A study done by [START_REF] Maestre | Increasing aridity reduces soil microbial diversity and abundance in global drylands[END_REF] showed that the diversity and abundance of soil bacteria and fungi are reduced as aridity increases in global drylands. A recent study by [START_REF] Jones | Moisture activation and carbon use efficiency of soil microbial communities along an aridity gradient in the atacama desert[END_REF], to investigate the response of soil microbial communities to water and carbon availability across an aridity gradient (semi-arid, arid and hyper arid) within the Atacama, showed that even under hyper arid conditions, very low levels of microbial activity and carbon turnover occur. This result is in line with the one obtained by [START_REF] Rabbi | Climate and soil properties limit the positive effects of land use reversion on carbon storage in eastern australia[END_REF] who showed that aridity has a strong negative influence on the soil organic carbon stock, based on a Structural Equation modeling applied to 1482 sites surveyed across the major agricultural region in Eastern Australia (AI varies from less than 0.07 (Eastern Australia coast) to more than 0.65 (moving away from the Eastern Australia coast)).

Impact of Aridity Index (AI) on topsoil F 14 C

The Aridity Index (AI) is detected as influential for topsoil F 14 C with a posterior inclusion probability of 53%. The impact of the AI obtained by the posterior mean is estimated at -0.029±0.022 with a 91% probability of being negative. In that respect, the increase of 0.1 of AI leads to a F 14 C decrease of 0.0029. In other words, the increase of 0.1 in AI will lead to a higher mean residence time of soil organic matter. Therefore, this result does not make sense since we expected that topsoil radiocarbon in the humid regions to be characterized by a higher F 14 C value (fresh carbon input) than in the arid regions.

To investigate this point further, we plotted the distribution of topsoil F 14 C (depth higher than 10 cm) according to the aridity index for all 131 profiles (Figure 5.20). This figure underlines a visually distinguishable increasing trend of the topsoil F 14 C with the increment in the AI up to a value of 3. In contrast, the green circle highlights a weird behavior of some points with a higher Aridity Index (AI = 4.223). A closer look at the database profiles shows that these points belong to the same 7 sites from Moor House Nature Reserve, that we already pointed out as not real mineral soils. They are further characterized by a very wet climate (MAP of 1665 mm) without efficient evaporation. They thus appear to undergo an equatorial monsoon climate whereas they are in a temperate region. This explains why the effect of the Aridity Index has a negative impact on the topsoil F 14 C instead of being positive. Chapter 5

Minimum monthly precipitation (min_P)

Reminder

The minimum monthly precipitation (min_P) was selected as an influential predictor only for incorporation depth (φ 3 ) of the F 14 C profile and for none of the latent variables of the carbon content profile.

5.3.6.1 Impact of the Minimum monthly precipitation (min_P) on F 14 C incorporation depth

The minimum monthly temperature was detected as influential for the F 14 C incorporation depth with a posterior inclusion probability of 55%. The posterior mean effect is estimated at -4cm with 7cm of standard deviation (70% chance of having a negative effect). Thus an increment of 1 mm in the minimum monthly precipitation decreases the F 14 C incorporation depth by 4 cm. At a constant level, F 14 C thus decreases with the temperature increases, reflecting the loss of 14 C enriched components, maybe some decades old compounds whose carbon derives from the nuclear bomb peak.

Soil radiocarbon and soil carbon content incorporation (φ 3 and ω 3 )

The estimation of regression effects corresponding to φ 3 (the F 14 C incorporation depth) and the ω 3 (the soil carbon content incorporation depth) indicates a poor estimation of what can happen in reality. This bad estimation can be explained by the non linear link between these latent variables, φ 3 and ω 3 , and the responses of F 14 C and soil carbon content respectively.

Soil type (Soil)

Recall Soil type was identified as influential for all latent variables for both F 14 C profile and carbon content profiles Many soils are marked by climate and type of vegetation [START_REF] Legros | Les grands sols du monde[END_REF]. For example, Gleysols and Podzols are characteristic of cold regions. Cambisols, Luvisols and Podzols are conditioned by temperate climates [START_REF] Spaargaren | Major soils of the world[END_REF]. Ferralsols and Plinthosols are tropical soils with forest cover. In addition, Chernozem, Kastanozems and Phaeozems are associated to steppe and grassland vegetation cover (http://www.isric.org), under climatic regimes that range from cool temperate to warm Mediterranean (https://www.britannica.com). Various soil physical and chemical properties such as the clay content are reported to control the organic matter decomposition rates [START_REF] Balesdent | Relationship of soil organic matter dynamics to physical protection and tillage[END_REF]. The clay and silt content is assumed to be positively correlated with the soil organic carbon [START_REF] Paul | Stabilization of recent soil carbon in the humid tropics following land use changes: evidence from aggregate fractionation and stable isotope analyses[END_REF].

Impact of soil type on mean residence time (F 14 C)

• Impact of soil type on deep F 14 C Jobbágy and [START_REF] Jobbágy | The vertical distribution of soil organic carbon and its relation to climate and vegetation[END_REF] underline that 56% of soil carbon globally can be found below 1 meter. The soil type was detected as influential for the deep radiocarbon response with posterior inclusion probabilities of 51%. The fusion of the levels by soil type obtained by the BEF (with k = 50) discriminates between the levels (see colors in Figure 5.22 and major lines in Table 5.16) and proposes a clustering that outputs 4 different groups.

Recap of the result of levels fusion for soil type for the deep F 14 C Chapter 5 Figure 5.31: Synthetic view of the dependence of soil F 14 C and carbon content on soil-climate-biome. Ten sites from the database were selected as representative of 10 major biomes, taking into account only the explanatory covariates detected as significant for soil F 14 C and soil carbon dynamics respectively. The orange (the green) band corresponds to the confidence in the local (within site) estimate of F 14 C (soil carbon content), and the gray band corresponds to the between-site variability of soil with similar environmental variables. Dark lines represent the sampled horizon of observed data and the blue points the real F 14 C (soil carbon content) measurements.

Briefly, it is interesting to note that simulated profiles are very close to measured data, for both carbon content and F 14 C profiles. Simulation reflects the general shape of the original profile but misses some specific features, as is expected from a model. The mean estimated topsoil F 14 C ranges between 0.97 and 1.22 and the mean estimated deep radiocarbon varies between 0.10 and 0.45 in all ten soils. Generally speaking and excluding the Andosol-type soils, the surface F 14 C was the highest in the warm tropical climates and was the lowest in the coldest climates.

The mean estimated carbon content ranges between 3.88 g/kg and 294 g/kg for topsoil and between 0.02 g/kg and 5.56g/kg for the deep soil. Soil carbon content estimation at the top of the profiles is associated to a very large variability that reflects the natural worldwide variability of carbon content [START_REF] Tifafi | Large differences in global and regional total soil carbon stock estimates based on SoilGrids, HWSD, and NCSCD: Intercomparison and evaluation based on field data from Usa, England, Wales, and France[END_REF]). It appears thus that even by forcing the system by specifying soil type, climate and land-use, the variability remains very high. The variability of the simulated profile decreases with depth.

It is worth noting the specific fate of some profiles. For example, the "Vertic" fate of the Vertisol that shows a specific F 14 C profile with a very deep F 14 C incorporation depth, associated, as shown by the carbon content profile, with deep incorporation of fresh plant-derived organic material. Chernozem and boreal Podzol are deep soils with a long history, close to the fate of permafrost. This is also reflected by the very low F 14 C in depth. In contrary Ferralsol and Cambisol are young soils with rapid turnover and thus associated to a higher F 14 C in depth. Luvisol, also, shows a specific structure. It exhibits good drainage and this is clearly expressed in the carbon content profile of the Mediterranean Luvisol where a drastic decrease in carbon content occurs at the upper levels of the profile. A study done by [START_REF] Noojipady | Forest carbon emissions from cropland expansion in the brazilian cerrado biome[END_REF] based on the satellite data on cropland expansion, forest cover and vegetation carbon stocks in the Cerrado biome showed that 29% of the carbon emission which is equivalent to 16.28 tg C yr -1 (tg = teragram = 10 12 grams), between 2003 and 2013, is due to the conversion of forest to cropland.

Chapter 5

A slight decrease from 1.20 to 1.18 is observed for the topsoil radiocarbon when converting the forest to cropland associated to a clear decrease in the topsoil carbon content from 26 g/kg for the tropical forest to 16 g/kg when converting forest to croplands (Figure 5.32). Tropical forest brings more fresh carbon into the topsoil (2200 tons per km 2 per year, Table 5.18) than the cultivated land (650 tons per km 2 per year) as it is characterized by a higher Net Primary Productivity (NPP) (Table 5.18) and shallow rhizosphere.

A significant increase in the deep F 14 C, from 0.45 to 0.58 is observed at deep layers. In other words, the residence time of the carbon at deep layers is greater for tropical rainforest than for cultivated land. But no real change occurs for carbon content in depth.

These results are in line with [START_REF] Balesdent | Atmosphere-soil carbon transfer as a function of soil depth[END_REF] study that shows that land use for crops reduces the incorporation of carbon into the soil surface layer but not into deeper layers. Our results add further elements to the discussion by [START_REF] Balesdent | Atmosphere-soil carbon transfer as a function of soil depth[END_REF] by showing that organic matter in our current soils is the legacy of its management by several generations of farmers. The much higher F 14 C in depth will impact the global carbon cycle over a long period of time with a more rapid return of stored (and thus not stocked) carbon to the atmosphere. Table 5.18: The variation of the Net Primary Productivity (NPP) per unit area according to the land type [START_REF] Jackson | A global budget for fine root biomass, surface area, and nutrient contents[END_REF].

Land type

5.5.1.2 Reforestation of temperate cropland and pasture leads to a higher carbon stock on short and longterm duration

Afforestation is commonly regarded as a mitigation solution to address climate warming thanks to an a priori high sequestration potential. However many studies have reported contradictory findings: afforestation results either in a decrease [START_REF] Garcia-Franco | Carbon dynamics after afforestation of semiarid shrublands: Implications of site preparation techniques[END_REF][START_REF] Li | Converting native shrub forests to chinese chestnut plantations and subsequent intensive management affected soil C and N pools[END_REF][START_REF] Wiesmeier | Depletion of soil organic carbon and nitrogen under Pinus taeda plantations in southern Brazilian grasslands (Campos)[END_REF], an increase [START_REF] Nave | Afforestation effects on soil carbon storage in the united states: a synthesis[END_REF] or no clear effect [START_REF] Eclesia | Shifts in soil organic carbon for plantation and pasture establishment in native forests and grasslands of South America[END_REF] in soil organic carbon stocks. Many factors have been highlighted to explain the extent of change in the soil carbon stock: site preparation, site management, fire, time elapsed since the conversion, previous land use, climate, type of soil, etc. but no univocal relation appears to link these factors to organic carbon content. Recently, [START_REF] Song | Global land change from 1982 to 2016[END_REF] reported the increase in land surface covered by forest in temperate regions due to loss of agricultural land and loss of pasture. What will be the consequence of this return of forest in temperate regions?

To address this issue, we selected the ten sites available in the database that met both temperate climate and ecosystem type equal to field (cropland -three profiles) or cultivated grassland (pasture -seven profiles ). We are positioned here in the case of reforestation following an agricultural decline rather than the afforestation of a region of weak vegetation (mostly savannah).

Chapter 5

Annual Temperature can be lost.

The impact of global warming as shown in Figure 5.35 differs in extent for both F 14 C and carbon content profiles depending on the climate subgroup. According to the Bayesian selection, an increase in MAT results in i-a decrease in topsoil carbon content that might be related to an increase in microbial activity and higher soil carbon mineralization, ii-no impact on deep soil carbon content, iii-no real change in top soilF 14 C, iv-variable (from negative to positive, including null) effect on deep soil F 14 C. Some of these mathematical results are well reflected in the distribution shown in Figure5.35. We do indeed find a lower carbon content under equatorial than under snow climate, and no real change in depth. However, the effect on topsoilF 14 C is noticeable and not as small as the Bayesian selection returned and the effect on deep soil F 14 C, while not univocal, is less variable, and a trend towards a decrease in F 14 C with increasing MAT is observed.

To address the global warming issue, an increment of 1 • C, 1.5 • C and 2 • C was applied to the 11 mean profiles that correspond to the 11 climate subgroups present in the database. The same procedure as previously described was applied: projecting each individual profile in a climate with MAT +1 • C, and evaluating the mean profile climate sub-group by climate sub-group. All the results are shown in Figure 23 and Table 4 provides key numbers. The following observations can be made:

• It is important to highlight that if an impact is recorded for a change of +1 • Within the impacted profiles, the profiles under a Steppe climate are those that change the most. This sub group is however represented by only one profile from the database. This soil shows a large change in deep soil F 14 C with an increase of 55% in the current value, from 0.11 to 0.17 resulting in a massive destocking of deep old carbon associated to warming. This result is accompanied by an 17% increase in deep carbon content and a shift towards deeper depth of the incorporation depth for carbon content. No conclusions can be drawn however from a single profile.

• Profiles associated to Equatorial savannah climate with dry summer also evolve with the MAT increase. It results in a 22% decrease in the topsoil carbon content and a shift towards a lower F 14 C of deep soil (change of 29%, from 0.48 to 0.34). This decrease in F 14 C is associated to a slight increase in the carbon content in depth. So a MAT increase seems to imply a higher plant input in topsoil, likely associated to higher vegetation NPP and surprisingly to a higher residence time in depth. Does this reflect a shift towards a profile with less bomb-peak derived carbon?

• The carbon content profile under snow climate with dry summer shows an increase in deep soil F 14 C, leading to a lower residence time in depth. This is associated to a small decrease in carbon content. No change is recorded in topsoil.

• The carbon content profile under Warm temperate climate with dry summer shows odd changes with a huge increase in the incorporation change (from 8.8 to 15.4cm). In parallel, an increase in deep soil F 14 C is recorded, i.e., a 17% lower residence time.

In summary, the impact of an increase in MAT is visible on soil carbon profiles and noticeable for half of the climate sub-groups. However, there is no univocal signal associated to a MAT increase in our data analysis. Impact can result either in an increase, a decrease or no change for both the F 14 C profile and carbon content profile and Chapter 5 for all of the latent variables. There is no global signal and special attention should be paid at the regional scale of equatorial climate with a dry season that shows a better storage capacity.

Chapter 5 Table 5.19: The latent values for both F 14 C and soil carbon profile corresponding to current temperature and an increase of the MAT by +1 • C, +1..5 • C and +2.5 • C. φ 1 and ω 1 refer to deep F 14 C and carbon content respectively, φ 2 and ω 2 represent the topsoil F 14 C and soil carbon content and finally φ 3 and ω 3 underline the F 14 C and carbon content incorporation depth.

C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C φ 1 0.
C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C φ 1 0.
C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C Current +1 • C +1.5 • C +2 • C φ 1 0.
Chapter 6

Summary and Conclusions

Soil carbon is a solution to mitigate global warming

Understanding the dynamics of soil carbon is a major challenge, especially as the IPCC pointed out the large uncertainty on the soil carbon stock and its potential impact on future climate change. The large capacity for carbon exchanges with the atmosphere, the huge uncertainties about the response in soil carbon to global changes in climate and land use practices and lastly the fact that soil carbon is the only pool that humans can manage, show the crucial global interest of better understanding the fate of soil carbon. Increasing the soil carbon stock is foreseen as a solution to mitigate global warming but this is relevant only if the storage is sustainable.

Incomplete physical representation of soil carbon dynamics by mechanistic models

Several mechanistic models have been developed to better understand the behavior of soil carbon based on expressing the physical processes by differential equations. Among these mechanistic models, one can mention the Roth-C Model, and the CENTURY model. These soil models are also integrated into Land Surface Models, such as ORCHIDEE, which is the land component of the Institut Pierre Simon Laplace (IPSL) Earth System Model. However, the incomplete view of the physical protection processes (spatial inaccessibility, encapsulation in organic macromolecules, etc.), the disregard of deep soil carbon layers and the parametrization of mechanistic models under specific management and climatic conditions highlight that there are still missing unknown terms in the soil carbon dynamics representation. Furthermore, improving the mechanistic models by integrating the missing physical concepts requires years and years of research work.

Faced with all these challenges, how can statistical meta-analysis help soil scientists?

Statistical meta-analysis helps decision making

Statistical modeling can provide faster responses to the scientific questions of today and the near future about soil carbon and it can consider uncertainties.

Using Statistical modeling to explore a soil carbon database

Before speaking about the statistical model, the database used for this study required long-term collecting and updating by Christine Hatté from 85 articles in the soil science, archaeology and paleoclimatology fields. The database contains 343 soil carbon profiles distributed roughly over the entire globe. Each profile displays radiocarbon and carbon content measurements at different depths. Furthermore, climatic and environmental conditions were reported (temperature, precipitation, soil type, ecosystem type, etc.). Unfortunately, the experimental design usually employed in a meta-analysis is not balanced. In our study, the only factors are the soil and ecosystem types. In order to improve the experimental design, we merged some categories that share the same features within the same factor. Whatever the efforts that can be made, we can improve but not avoid this issue since the soil type and ecosystem are naturally associated. The database provides a homogeneous representativeness of intermediate climates, i.e. tropical, warm temperate and snow climates, leaving out extreme climates (arid and polar).

The database has the advantage of providing soil radiocarbon and carbon content for each profile. A parallel comparison of these two profiles afforded a twofold vision of soil carbon content and its mean residence time since the radiocarbon is characterized by its radioactive decay. In addition, among the profiles of radiocarbon ( 131Chapter 6 profiles) and carbon content (125 profiles) remaining after database cleaning (removing paleosol, profiles with unknown information, etc.), 58 profiles have measurements for depth levels exceeding 100 cm. Thus, the deep soil carbon is well represented in the database.

Bayesian inference better represents model uncertainties

The variation in both soil radiocarbon and soil carbon content with depth are separately represented by a hierarchical non linear model with latent variables which in turn are linked to the climatic and environmental potential predictors. The estimation of unknown model parameters was done using Bayesian inference. Bayesian inference gives us the possibility to integrate expert advice about the model parameters. Unfortunately, sometimes, it is not easy to provide an informative prior, either because no prior beliefs on model parameters are available or because the soil scientist prefers to be objective and see what we can learn from data. Thus, for each soil carbon profile, Bayesian inference makes it possible to compute the credible intervals which highlight uncertainties within and between soil carbon profiles. The statistical model has the advantage of predicting the soil carbon profile for both radiocarbon and carbon content for sites where no data have been collected. We only need to know the climatic and environmental site conditions.

Bayesian selection methods provide a probabilistic judgment about the contribution of climatic and environmental factors to soil carbon dynamics

As a matter of fact, soil scientists have an idea about the potential climatic and environmental factors that can impact the soil carbon dynamic but they do not know, in the first place how to prioritize these factors by their effects on soil carbon dynamics and in the second place there is still debate on some issues. For example, the soil science community is divided on the question of the soil carbon sensitivity to temperature changes. Some soil scientists think that the topsoil carbon is more sensitive to temperature changes than the deep soil carbon, others consider that the top and deep soil have the same sensitivity to temperature changes and the remaining ones think that the deep carbon is more influenced by temperature changes than the topsoil.

The first goal achieved in my thesis was to provide a probabilistic judgment on inclusion of the climatic and environmental factors which have a physical interpretation in the latent layer models. I first explored the Stochastic Search Variable Selection (SSVS). This approach introduced by George and McCulloch (1993), was designed to select numerical predictors in the framework of a linear model. However, even after adapting the SSVS to support the inclusion of the categorical predictors such as the soil and ecosystem types used in our study, the results obtained on artificial data, generated according to the proposed statistical model for soil carbon dynamics, proved that the detection of some significant categorical predictors with SSVS can be misleading. Thus, the second challenge achieved was to investigate other Bayesian selection approaches appropriate for categorical predictors. To ensure the best selection of the categorical and numerical covariates, three Bayesian Selection approaches were explored: Bayesian Group Lasso with Spike and Slab priors (BGL-SS), Bayesian Sparse Group Selection (BSGS) and Bayesian Effect Fusion model-based clustering (BEF). In addition to selecting categorical predictors, the BSGS also provides a selection by level within the same predictor and the BEF makes it possible to merge the levels within the same predictor having the same effect on the variable on interest.

The best sub-sets of climate and environmental factors selected were obtained by comparing the previous selection methods based on the Bayesian selection Criteria such as the Deviance Information Criterion (DIC), the 5-fold cross validation, etc.

The best sub-sets of selected climate and environmental factors show very interesting findings:

Chapter 6

1. The non detection of the F 14 C of the sampling year for topsoil F 14 C is explained by a non representative distribution of database profiles by sampling year.

2. The deep soil layer radiocarbon is identified as more sensitive to temperature changes than the topsoil radiocarbon.

3. The Mean Annual Precipitation affects both topsoil and incorporation depth in radiocarbon and carbon content profiles.

4. The Aridity Index is only an influential predictor for topsoil and for incorporation depth of the carbon content profile.

5. The result of merging soil types for both topsoil F 14 C and deep F 14 C underlines that the F 14 C is mainly dominated by the climate/vegetation and soil texture at the topsoil and by clay content for deep layers.

6. Soil type and land use affect both topsoil and deep layers of radiocarbon and carbon content profiles.

Besides these very encouraging results, it appears that the impact of a climatic or categorical factor on the latent variable representing separately the incorporation depth for F 14 C and soil carbon content (φ 3 and ω 3 latent variables) is highly overestimated.

The statistical model provides prediction in a context of global changes

The third achieved goal was to illustrate how the statistical model can be used to predict the changes occurring in the soil carbon profile when the land use or climatic conditions change. Here we show the impact of land use change with two examples: 1-when equatorial forest is replaced by cropland and 2-when cultivated grassland and field are replaced by cultivated forest in a temperate region.

The first example shows that the conversion of equatorial forest to cropland impacts both topsoil soil carbon content and deep radiocarbon dynamics. After conversion, the topsoil carbon content decreases from 26 g/kg to 16 g/kg but shows a decrease in the mean residence time of deep soil carbon (deep soil radiocarbon increases from 0.45 to 0.58). For the second example, the reforestation of temperate cropland and pasture yields a higher carbon stock and long-term duration. A third example highlights the impact of global warming on soil carbon dynamics.

Radiocarbon and carbon content profiles plotted according to the Köppen-Geiger climate classification showed that the deep carbon is more stable under a snow climate with dry summer and has a lower mean residence time under an equatorial savannah with dry winter climate. The temperature increment is more visible on the carbon profile on going from the current temperature to an increase of +1 • C than on going from an increase of +1.5 • C to +2 • C. The impact of temperature increase may be negligible for soil in a temperate climate but considerable for both deep and topsoil in a snow climate.

Finally, there are still some ideas and objectives that I have not had time to look at in detail and that I will discuss in the "Perspectives" section.

Improvements and Perspectives

Several propositions and possibilities can be investigated in order to improve the structure of the statistical model and to a further outlook for the use of the statistical model. I distinguish here, the following points:

• A bivariate Gaussian Process for both radiocarbon and carbon content dynamics can better express model uncertainties Chapter 6

• The database can be better managed -Due to lack of time, the analysis of soil bulk density variability required to convert carbon content into soil carbon stock could not be carried out. This step will have to be done to project the results at the global scale. It might consist of studying measured bulk density within the database or of using pedotransfer functions that differ with the type of soil rather than the universal Alexander equation.

-At the interpretation steps, some profiles were highlighted outside the general trend. These include, in particular, some profiles in the UK under very wet climate. These profiles likely biased the variable selection and highly likely the extent of impact value. A second set of evaluations can be conducted without these "outlier" profiles to get a statistical model that would better mimic the general trend.

• It is necessary to improve the experimental design of the database Another challenge will be to improve the current database. As we noted in Chapter 2 that describes the database and later on for the physical interpretation of the variable selection and fusion, the database is far from being perfect. For example, we don't have the same number of profiles by soil and ecosystem type. 37% of the database profiles are forest while only 8% are defined as cultivated grassland. Even for soil type, 9% of the database profiles are defined as Andosols and 7% as Regosols/Arenosol/Leptosol. However, only 1% of the total continental land area on earth is occupied by Andosols and 22% by Regosols/Arenosol/Leptosol.

Furthermore, as soil type and ecosystem are associated, it might also be of interest to divide Gleysol into two categories: tropical and boreal Gleysols. Some of their features are similar as they have the same name but some others, especially interaction with vegetation, are different.

Increasing the number of profiles by soil and ecosystem (land use + vegetation) category is not the only important point, however; there are several crucial features that also need to be addressed such as: the total continental land area occupation by soil and vegetation cover, the association between: soil, vegetation cover and the climate conditions. In addition, the current database does not contain profiles from the arid and hyper arid climate classes. This is a crucial lack, especially for the use of the model in a prediction mode in these regions that are particularly vulnerable to climatic changes. As stated in the database description, the inhomogeneous distribution of the sampling years in the database hampers a good representation of the F 14 C profile, in particular for the latent variable that mimics the incorporation depth (φ 3 ).

A further point to improve in the database is thus the distribution of the profiles by sampling year. In the current database, the profiles sampled between 1990 and 2000 represent 53% of the database profiles. As a consequence, the atmospheric radiocarbon related to the bomb peak in the early 1960s was not detected as influential for the topsoil F 14 C nor for the incorporation depth.Thus, why not consider archives samples, such Finally, the question remains as to how to consider permafrost in the soil carbon statistical model. In fact, permafrost accounts for about 20% of the Earth's land surface. This type of soil is particular since its temperature subsurface never rises above 0 • C for at least one year.

A deeper analysis of impact of the level width on the profile modeling might also help in defining the sampling pattern by pointing out the best compromise between fieldwork and analytical work and the best representation of the profile specificities. It might thus result in advice on the maximum depth to be reached according to the type of soil, for the level thickness according to the depth (likely a finer resolution for top and mid depth soil

Chapter 6

and a lower resolution for deep soil).

• The extrapolation of the Bayesian statistical model developed for the soil carbon content is useful to have a global (or regional) estimation of the soil carbon stock

The Bayesian statistical models for radiocarbon and soil carbon content dynamics allow us to predict the soil radiocarbon and carbon content profiles of a new site, knowing the corresponding climatic and environmental information. As no measurements is provided for this site, wider credible intervals will be observed for the unknown model parameters.

In the first place, the predicted profiles, when changing the land use or climatic conditions, are obtained without taking into account the observed measurements (Chapter 5.5.2, Section 5.5). Which means that the site is considered as a new site. Thus, new latent variables are generated for the soil radiocarbon and soil carbon content models. These considerations raise the following question: how can we predict the profile with the changes in climatic or environmental conditions, taking into account the measurements observed for the corresponding site? We can imagine that we have to add to the current and estimated latent variables under the ambient climatic and environmental conditions for site s, the change of effect due to replacing a forest by a cropland or to increasing the temperature by 1 • C. In contrast, following this proposition does not guarantee the positivity constraint of model latent variables.

In the second place, one can wonder how we can extrapolate the Bayesian statistical model to have a predicted profile of soil carbon content or radiocarbon at regional or global level. With a more complete database and a more balanced experimental design, other doors open. It then becomes possible to apply the statistical model designed on the whole database to models for each type of soil and each climate zone. This would increase the power of the projection capacity of the study. It would make it possible to better decipher the impact of land-use change according to the soil type, and to better project the impact of present global warming according to the world regions. GIS is thus possible.

Here, we are talking about an extrapolation to 3 dimensions: longitude, latitude and depth. Digital Soil Mapping (DSM) using contextual information spatial models in deep learning is very popular and can be used to generate maps [START_REF] Mcbratney | On digital soil mapping[END_REF]. In fact, there are deep learning methods such as convolutional neural networks which expand the classical DSM approach by including information about the vicinity of a site. Each site is characterized by n climate and environmental covariates with a 3-dimensional array for width, length in pixels of a window centered at point p (site coordinates) and the covariates. Multi-task learning can handle the notion of depth by providing prediction layer by layer. More research and investigation should be done to explore how to apply these deep learning approaches can be challenging competitions of our hierarchical non linear structure model with latent variables. The possibility of extrapolating the statistical model is very useful since optimizing the sampling design takes a lot of time and it is also expensive (acquisition of data and processing samples in the laboratory).

• Bayesian selection approaches can help to better understand the outputs of the mechanistic model for soil carbon dynamics

The coefficient of diffusion, which underlines soil bioturbation, and the advection coefficient, which is related to lessivation, are treated as constants in the mechanistic models developed for soil carbon dynamics. However, in reality, these coefficients are not constant and vary with depth. A challenge will be to transform these constant coefficients to functions that decrease exponentially with depth. The explored Bayesian selection methods explored can be used in order to define the influential climatic and environmental factors.

Conclusion

Moreover, one can distinguish the effect of isotopic fractionation. In fact, plants fix 12 C more easily than 13 C and 14 C during photosynthesis. In other words, the vegetation has not the same radiocarbon concentration as the atmosphere that it grew in. To compute the isotopic fraction of 14 C

12 C we must calculate the isotopic fraction of 13 C

12 C since it is easy to measure and is not driven by time, unlike 14 C. The 14 C fractionation factor is approximately given by the square of the 13 C fractionation factor:

14 C 12 C = 13 C 12 C 2 .
How to report radiocarbon activity: Several units have been proposed in the literature to express radiocarbon activity. The choice differs from one community to another. Whereas paleoclimatologists and archeologists will prefer yr BP or cal. yr BP, geoscientists dealing with recent samples will prefer pMC, F 14 C and ∆ 14 C.

To better understand the physical and the philosophical differences between the proposed units, it is necessary to come back to the physical measurement and briefly to the history of 14 C.

Let us assume A = 14 C/ 12 C (or 14 C/ 13 C depending on the type of measurement). The rate normalized to δ 13 C of -25 A N is defined as A N = A(1 -2 * (25 + δ 13 C))/1000). We will note A SN as the sample activity normalized to -25 and A ON as the atmospheric ratio in 1950. It is retrieved thanks to two international standards: HOxI (only a few labs still use this one) and HOxII [START_REF] Donahue | Isotope-ratio and background corrections for accelerator mass spectrometry radiocarbon measurements[END_REF]. It is equal to A ON = 0.95 × A HOxI or A ON = 0.7459 × R HOxII . Furthermore, the activity of the standards decreases with time and thus depends on the sampling year. This can be corrected by considering A abs defined as:

A abs = A ON * exp[λ (y -1950)] with λ = 1/(ln(2) * 5730)yr -1 (7.2)
• year BP (Before Present)

The first unit used is the conventional age expressed in units of years Before Present (BP). Standard practice is to use 1 January 1950 as the commencement date of the age scale. There is no particular reason for choosing the year 1950 as the reference except that it is before the bomb peak and it refers to the publication of the first radiocarbon dates. The conventional age (BP) considers the radiocarbon decay equation given in equation 7. • Percent Modern Carbon, pMC Because year zero is 1950, the conventional age calculated for a sample that is more recent than 1950 will be negative. This does not make sense for an age and can even cause misunderstanding. For that reason, this unit cannot be used for modern data analysis.

In 1972, the community agreed to use pMC to replace yr BP for modern samples. This unit is considered as a percentage of the ratio of the sample's normalized activity ASN versus the absolute normalized activity A abs , which corresponds to the specific activity of atmospheric carbon in the year 1950. This unit is specially used for post-bomb samples.

pMC = A SN /A abs

Problems arise when the dataset includes both old and modern samples as old samples might be reported as pMC = A SN /A 0N , i.e. without considering the decay between the measurement year and 1950. This unit should thus be avoided for datasets such as the one we deal with here.

• Per mil depletion, ∆ 14 C ∆ 14 C is a very useful way of reporting 14 C measurements for geochemical studies (Reimer, 2004). Under Stuiver and Polach's (1977) definition, the ∆ 14 C of a sample depends on the year in which it is measured whatever the age of the sample. For example, a sample formed in 1977 will give a different ∆ 14 C if measured today versus if it had been measured in 1977.

∆ 14 C = A SN e λ c (y-x) A abs -1 .1000% (with age correction) y = year of measurement;

x = year of formation or growth;

λ c = (1/ln2*5730) year -1 .
A SN : the normalized specific radiocarbon activity for the sample;

A abs : specific activity of atmospheric carbon of year 1950, measured in 1950 A problem might arise since ∆ 14 C is also used by the oceanography community to refer to the shift between the 14 C content in the ocean and the 14 C content in the atmosphere on the same day. As the same symbol has several meanings it might lead to confusion and one prefer to avoid it for this study.

• Fraction Modern, F 14 C Reimer et al. ( 2004) highlight an alternative unit that does not depend on the year of measurement and is corrected for isotope fractionation. They propose as an alternative solution to use the F 14 C unit : Here is a modern rewriting of Bayes' formula:

F 14 C = A SN /A
Pr(θ |Y ) = Pr(Y |θ ) × Pr(θ ) Pr(Y ) (7.3)
A few years after the death of Thomas Bayes, and independently, the French mathematician Pierre Simon de Laplace, rediscovered Bayes' formula which makes it possible to evaluate the relevance of what one believes one knows (θ ) in the light of the information provided by observations (Y ). Bayes, a clergyman like his father and an amateur mathematician, did not publish the one scientific paper for which he is nowadays famous, whereas Laplace resisted his father's desire for him to become a clergyman, embraced the career of -what could now be called-a professional researcher and published an impressive quantity of scientific papers in several domains. He made many advanced applications, as diverse as celestial mechanics by predicting the motion of planets and their satellites through the first statistical analysis of astronomical data, or demography by working from samples of the number of births, marriages and deaths, or reliability by studying testimony in court, and so on.

In the centuries that followed, interest in the statistical learning equation 7.3 seemed somehow to fall dormant in the academic world, despite many remarkable applications in engineering and operations research. Does history remember that Poincaré's plea for Dreyfus's innocence was based on Bayes' formula? That it was the key that allowed Alan Turing and his colleagues to break through the encrypted messages of the Enigma machine used by Nazi Germany in its military communications, giving the allies a decisive advantage? The book by McGrayne (2011) traces back through centuries the stunning story of Bayes'rule and its important, although almost secret, operational applications.

Bayesian reasoning

As a very striking example, let's skip two hundred and fifty years from Bayes and Laplace into our present world and consider, as [START_REF] Raftery | Less than 2 • c warming by 2100 unlikely[END_REF] did, the following hypothesis θ : The global temperature increase of our planet will be less than two degrees by the end of the century. Bayes' formula shows that, if you wish to infer the probability of this hypothesis θ , knowing that it is some set of observations Y that has occurred, you only need to be able to quantify the probability of this collection of observations Y under hypothesis θ (and, in this sense, in general, the reasoning of physicists goes very easily from cause to observable consequences), to weight it by the chance that you agree with cause θ before observations Y are revealed, then to re-normalize this product by the probability Pr(Y ).

We have here all the terms of eq7.3. The future emission Y of greenhouse gas CO 2 can be viewed as the product of the Gross Domestic Product per capita (GDP), the carbon intensity (CO 2 emission per unit of GDP) and the world population : these three driving components can be calibrated on the last 30 years of data and assessed for UN population predictions up to 2100 . The authors also modeled the links from θ to Y to calculate Pr(Y |θ )

with mild hypotheses: a decreasing trend for carbon intensity was considered and the GDP could not overcome 2003; [START_REF] Wikle | Hierarchical models in environmental science[END_REF][START_REF] Clark | Applications of Computational Statistics in the Environmental Sciences: Hierarchical Bayes and MCMC Methods[END_REF][START_REF] Cressie | Accounting for uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical modelling[END_REF].

HBM works through the conditional decomposition of high-dimension problems into a series of probabilistically linked simpler substructures. HBM makes it possible to exploit diverse sources of information to derive inferences from large numbers of latent variables and parameters that describe complex relationships while keeping as close as possible to the basic phenomena.

The three basic layers of hierarchical statistical models are as follows:

1. A parameter (often called θ ) level identifying the fixed quantities that would be sufficient, were they known, to mimic the behavior of the system and to produce new data statistically similar to the ones already collected.

2. A latent process level Z depicting the various hidden mechanisms (given the parameters θ ) that make sense of the data;

3. A data level that specifies the probability distribution of the observables at hand (Y ) given the parameters (θ ) and the underlying processes (Z) ; HBM stands out as an approach that can accommodate complex systems in a fully consistent framework and can represent a much broader class of models than the classical statistical methods from ready-to-use toolboxes that can be found in the frequentist literature. Eq 7.3 still applies by changing the notations and considering the block of unknowns (θ , Z) as the term that was formerly named θ .

Computational Bayes

Since the turn of the century, Bayes' formula has also revolutionized [START_REF] Brooks | Bayesian computation: a statistical revolution[END_REF]) the scientist's toolbox: simulation methods, also known as Monte Carlo techniques in statistical jargon, made possible by the successive waves of (personal) computer advances , are effectively helping to disentangle the complex networks of causes that appear in the scientific challenges at the beginning of our third millennium.

Obtaining the posterior distribution of the unknowns (the term Pr(θ |Y ) in equation 7.3) is a more difficult task than it appears at first sight : the denominator of the formula ((the term Pr(Y ) = θ Pr(θ |Y ) × Pr(θ )dθ ) is an integral which needs to be performed on the support of the unknowns but θ may vary over a high dimensional space!

Conclusion

Brute force numerical integration does not work beyond dimension 3 , and for long, only reduced dimension models -not to say baby models -were under study in the Bayesian paradigm. An alternative approach was models for which the association between the prior Pr(θ ) and the likelihood Pr(Y |θ ) appears in closed form; these are the so-called conjugate models. Working in such restrictive ad-hoc families of models curbed the impetus of the lay Bayesian modeler from the twentieth century. However algorithms to tackle the problem of generating random draws from equation 7.3 without computing its denominator were already known by the end of World War 2 because a lot of effort had been devoted to their development in the Manhattan project at Los Alamos where there were computers at the time, but computer power in ordinary research laboratories was not yet sufficient to allow for their routine implementation. It is only since the turn of the century (around 1990) that Markov Chain Monte Carlo (MCMC) algorithms, with the advent of personal computers, became easily available worldwide for scientists. For details on the modern theory of Bayesian inference, see [START_REF] Brooks | Handbook of markov chain monte carlo[END_REF] and [START_REF] Neal | Mcmc using hamiltonian dynamics. Chapter 5 of Handbook of Markov Chain Monte Carlo[END_REF].

Monte Carlo integration

Bayesian learning requires evaluating the expectation of functions of the unknowns θ with respect to the posterior distribution. Nowadays, Monte Carlo algorithms avoid the explicit computation of the Bayes formula denominator since they only need to know the distribution from which to sample, up to a constant. To focus on the technical problem , we will now write π(θ ) for the posterior probability density (because from now on we consider that the data are given once and for all and need not be recalled in Pr(θ |Y ) and review some techniques to evaluate the essential characteristics of π(θ ) such as

E f = θ π(θ ) f (θ )dθ
for any measurable function f when π(θ ) is known up to a constant. For instance, if the test function f (θ ) = θ one would obtain E f as the posterior mean, and more generally the moments of order p (when they exist) by taking f (θ ) = θ p . Such an expectation of f can be approximated by the Monte Carlo method, using a G-sample (θ (1) , . . . , θ (G) ) of independent values generated from π

E f ≈ 1 G G ∑ g=1 f (θ (g) )dθ
Unfortunately, a method to directly draw an independent sample of θ (g) from π is generally not available for π which is often a multivariate sophisticated distribution ( it is the unnormalized product of the likelihood and the prior)

. But the approximation will still work when the θ (g) are dependent, as long as the dependence is not too stringent, such as the one given by a Markov chain (with specific properties easy to enforce in the case of homogeneous Markov chains). An ergodic Markov chain, with π as a (unique) invariant distribution, will converge to π from any initial distribution of states given to the chain. MCMC are ergodic Markov chains designed to stochastically visit the support of θ , quickly disregarding the initial values of the chain (the so called burn-in period) and moving rapidlyi.e without too many correlations between the successive steps of the algorithm -so as to reach the equilibrium distribution π and make it possible to have access to the posterior distribution by means of the random sample of replicates {θ (g) } g=1:G [START_REF] Brooks | Handbook of markov chain monte carlo[END_REF]Gelman et al. (2013a); [START_REF] Roberts | Markov Chain Concepts related to sampling algorithms[END_REF]; [START_REF] Robert | Monte carlo methods[END_REF]).

The free software WinBUGS [START_REF] Lunn | Winbugs -a Bayesian modelling framework: Concepts, structure, and extensibility[END_REF][START_REF] Congdon | Bayesian Statistical Modelling[END_REF][START_REF] Spiegelhalter | Winbugs version 14 user manual[END_REF][START_REF] Ntzoufras | Bayesian modeling using WinBUGS[END_REF] is a tool of choice for Bayesian inference when the posterior distributions contain high-dimensional latent system states and parameters, which is the case of HBMs. OpenBUGS is the open source variant of WinBUGS which offers greater flexibility and extensibility. JAGS is a program developed by Martyn Plummer that relies on the same modeling grammar and language as WinBUGS, but its main advantages are its speed (since it is written in C++ instead of Component Pascal) and its platform independence (with regards to Windows systems) [START_REF] Plummer | Jags: A program for analysis of Bayesian graphical models using Gibbs sampling[END_REF]. Many step-by-step Bayesian initiation course tutorials have been developed using OpenBUGS/JAGS and their extensions (GeoBUGS, PkBUGS, etc.). ) is designed to deal with hierarchical models and the many problems they raise. The NIMBLE creators try to fill the gap between an abundant literature on these questions but a limited software offer that does not allow scientists to write their own blocks of the inference routine. The idea underpinning NIMBLE is to allow both a flexible model specification and a programming system that adapts to the model structure. Some frequentist authors even resort to Bayesian computation as an ultimate tool in the hand of the analyst to find solutions to problems not solved by other conventional analyses, see for example

Chapter 17 of [START_REF] Zuur | Mixed Effect Models and Extension in Ecology with R[END_REF] or the many papers on data cloning by [START_REF] Lele | Data cloning: easy maximum likelihood estimation for complex ecological models using bayesian markov chain monte carlo methods[END_REF].

Gibbs sampling

Gibbs sampling is the simplest MCMC method, widely used for Bayesian inference following its introduction by [START_REF] Geman | Stochastic relaxation, gibbs distributions and the bayesian restoration of images[END_REF] and [START_REF] Gelfand | Sampling-based approaches to calculating marginal densities[END_REF]. It is applicable when the multi-dimensional parameter θ = (θ 1 , . . . , θ j , . . . , θ p ) is such that one can easily sample from each of the full conditionals of the posterior distribution. The full conditional for one component θ j of θ is the conditional distribution (with respect to π) of this j th component given values for all the other components of θ . The Gibbs sampler generates a MCMC algorithm1 by stacking stochastic transitions from θ (g) to θ (g+1) relying on the full conditionals as follows:

• Pick θ

(g+1) 1 from the distribution of θ 1 given θ (g) 2 , . . . , θ (g) p • Pick θ (g+1) 2 from the distribution of θ 2 given θ (g+1) 1 , θ (g) 3 , . . . , θ (g) p
• . . .

• Pick θ (g+1) j from the distribution of θ j given θ (g+1) 1 , θ (g+1) 2
, . . . , θ

j-1 , θ (g) j+1 . . . , θ (g+1)
• . . . The Gibbs sampler has the favor of many practitioners since no tuning is necessary to run the MCMC algorithm. It grounds the WinBUGS expert system (BUGS meaning Bayes Using Gibbs Sampler) and the JAGS (Just Another

Gibbs Sampler) package. The Gibbs sampler takes advantage of the conditional independence structures encountered in the direct acyclic graph nodes of HBM to implement more efficient sampling in the parameter space [START_REF] Lauritzen | Local computations with probabilities on graphical structures and their application to expert systems (with discussion)[END_REF]: many conditional distributions are of standard forms for which efficient sampling procedures are readily available. However Gibbs sampling may be rather slow since drawing in conditionals generally yields only small moves in the parameter space and this drawback worsens with the parameters' dimension.

Metropolis

The Metropolis algorithm was first presented in the seminal paper of [START_REF] Geman | Stochastic relaxation, gibbs distributions and the bayesian restoration of images[END_REF] but statisticians had to await the development of personal computers for a generalized diffusion of this sampling routine, able to generate random draws from almost any distribution known up to a constant. Conversely to the Gibbs sampler, the Metropolis algorithm needs an auxiliary distribution (or jump function J(θ ′ → θ ) ) to make proposals that will be accepted or rejected. In the latter case, the new value is a replication of the previous state of the algorithm. The basic Metropolis generates an MCMC algorithm 1 by iterating singular probabilistic transitions from θ (g) to θ (g+1)

as follows:

• Generate a candidate θ ⋆ from J(θ (g) → θ ⋆ )

• Accept the candidate θ ⋆ with probability min(1, π(θ ⋆ ) π(θ (g) ) )2

Going Bayesian

The application of Bayesian concepts and methods was for long poorly developed because of the complexity of posterior calculations and the absence of closed form solutions. Advanced numerical computation for the implementation of Bayesian methods (Monte Carlo Markov Chain simulation methods: Gibbs Sampler, Metropolis

Hastings, Hamiltonian MCMC, Population Monte Carlo, etc. ) nowadays allow us to overcome this technical obstacle . This is reflected in a marked upsurge of Bayesian methods in all statistical applications [START_REF] Hoff | A first course in Bayesian statistical methods[END_REF].

When searching for the term " Bayesian" in PubMed, one can observe an exponential growth in the number of publications, which has doubled about every five years since 1963! Bayesian models and algorithms are now widely applied in all domains of engineering and machine learning [START_REF] Barber | Bayesian reasoning and machine learning[END_REF], as well as neural and psychological models. There even is a theory of the Bayesian brain [START_REF] Doya | Bayesian brain: Probabilistic approaches to neural coding[END_REF], following the seminal work of Dehaene (2011).

Of course, the merits of Bayesian reasoning are grounded on philosophical arguments regarding uncertainty such as the ones detailed in [START_REF] Kadane | Principles of uncertainty[END_REF] and [START_REF] Lindley | Understanding uncertainty[END_REF], but people tend to attach more importance to the demonstration of Bayesian effectiveness in practical contexts. Among the many advantages of Bayesian statistical inference that explains its seemingly inexorable success, the following ones are often pointed out:

• Bayesian learning is simply performed by applying the rules of probability and data assimilation is nothing but updating the scientist's probability distribution over all the unknown quantities. Bayesian inference offers a direct probabilistic interpretation for credible intervals and tests, as well as fruitful insight into frequentist estimators and statistical decision rules.

• The uncertainties on all unknowns are taken into account in the scientific reporting of the study and the predictive applications of a model, often delivering a less optimistic point of view than their frequentist counterparts.

• The missing values can be easily imputed using Bayesian MCMC simulation methods. The analysis is often more robust to outliers, by using more flexible distributions.

• For small sample surveys, Bayesian analysis can help to increase precision by bringing neighboring information or prior expertise into the analysis. For large sample surveys, central limit theorems often make it possible to obtain the frequentist classical results as asymptotic limits of Bayesian ones .

Bayesian model checking and Bayesian model comparison A-Bayesian model checking

After specifying the prior distributions based on the expert advice, constructing a reasonable model to explain the data and computing the posterior distributions of the unknowns based on the MCMC simulation methods mentioned previously, a last remaining step is to check the adequacy of the model fit.

A crucial point to check is the robustness of the results to the choice of prior distributions by a sensitivity analysis . It is important to examine how the posterior distributions change when another prior distribution is used instead of the present one (Gelman et al., 2013a). In fact, the prior distribution should represent our probabilistic judgment about the values that the parameter to be estimate may have. However, it is sometimes very difficult to come up with precise distributions to be used as prior. Thus, the analyst would like to test all possible combinations of prior distributions and likelihood functions selected from classes of priors and likelihoods considered empirically plausible. The Bayesian model is considered to be robust if the posterior distributions seem invariable with the different priors. If the posterior results differ substantially, the analyst must either make the prior informative using expert advice or build several prior + likelihood assemblies and specify the set of modeling assumptions considered Conclusion for each case. For example, in the case of Bayesian selection methods, it is useful to study the effect of fixing the prior inclusion probability "p" to 0.5 or proposing a prior distribution as a Beta [START_REF] Katsuno | A study of the carbon dynamics of Japanese grassland and forest using 14C and 13C[END_REF][START_REF] Katsuno | A study of the carbon dynamics of Japanese grassland and forest using 14C and 13C[END_REF] or a Uniform(0,1) distribution.

A common choice among Bayesian selection approaches is to set the prior inclusion probability to 0.5 as this choice makes all models equiprobable but on the other hand, it favors models where about half of the variables are selected. This choice, for example, is problematic when a large number of explanatory covariates are included in the model.

Another point is to investigate if the inferences from the model make sense. Sometimes, the expert does not intervene to specify prior distributions for reasons of convenience or objectivity or time, thus the posterior results may not reveal what they expect. For example, for the soil scientist, atmospheric radiocarbon is an influential predictor of topsoil radiocarbon. However, the application of a Bayesian selection model does not detect atmospheric radiocarbon as a significant predictor. This result may not make sense for soil scientists but there may be some explanations for this lack of significance. For example, on the one hand, the natural variation in atmospheric radiocarbon is not constant with time and on the other hand, it is marked by an artificial increase from above-ground nuclear testing from about 1950 until 1963. Consequently, it is clear that the database is not representative of the actual variation in atmospheric 14 C (the majority of data were sampled after 1963). Furthermore, the prior inclusion probability is fixed at 0.5 which means that the same probability is given to select or remove the predictor. It is not a good idea to set non informative priors and rely on Bayesian inference as a data analysis approach. In such a case, to achieve a more coherent statistical selection, the expert should first express his own beliefs on the prior inclusion probability for atmospheric 14 C , setting the prior inclusion of this predictor to a value (higher than 0.5) that corresponds to his own belief, and see if the updated posterior distribution after data assimilation goes against his prior judgment .

Furthermore, investigations should be carried out to test whether the model is consistent with the data. In fact, a model is consistent if replicated data generated from the model look similar to the observed data. An important discrepancy can highlight misfitting of the model. The discrepancy between model and data can be calculated based on the tail-area probability or the so-called Bayesian p-value applied to the test quantities T (y, θ ) and T (y rep , θ ) , where y represents the observed data, y rep represents the replicated data and θ the unknown parameters that we seek to estimate (Gelman et al., 2013a). The Bayesian p-value is defined as the probability that the replicated data will exceed and be more extreme than the observed data:

p B = Pr(T (y rep , θ ) ≥ T (y, θ )|y)
If the model is true or close to true, the posterior predictive p-value will almost certainly be very close to 0.5 (Gelman et al., 2013b). However, a clear guide-line on how best to interpret the p-value is not available. In fact, the posterior predictive posterior check remains one of the most misunderstood and confusing concepts in statistics.

Gelman stressed the fact that a p-value -say equal to 0.2 -is not, and should not be interpreted as a claim that the model is "true"; rather, it should be interpreted as a statement that the model (probabilistically speaking) fits one particular aspect of the data (summed by the T statistics).

B-Bayesian model comparison

The usefulness of the model is evaluated by its ability to provide a good data fitting and prediction about the future. There exist several Bayesian criteria to compare different models in competition. We distinguished among them:

• The Deviance Information Criterion (DIC) is a measure of predictive accuracy. This criterion is defined as follows:

DIC = 2 Davg (y) -D θ (y) Davg (y) = 1 L L ∑ l=1 D(y, θ l ) and D θ (y) = D(y, θ (y))
where D(y, θ ) = -2log(P(y|θ )) and L is the number of iterations. θ is fixed to the posterior mean of θ under the posterior distribution. The estimated average Davg is a better summary of discrepancy than D θ of the point estimate since it takes into account the model uncertainties. The DIC is a popular criterion among the community of Bayesian practitioners, used to compare models through its easy implementation in the graphical modeling package BUGS/JAGS/Stan. The model with the smallest DIC is considered to be the model that would best predict a replicated dataset with the same structure as the currently observed one (for more details see the paper by [START_REF] Spiegelhalter | Bayesian measures of model complexity and fit[END_REF]). However for most Bayesian theorists, this criterion is not theoretically grounded since it is defined using a pointwise estimate (the posterior mean estimation θ (y)) without any probabilistic justification.

• Widely Applicable Information Criterion (WAIC) is seen as an improvement of the DIC. As mentioned previously, the DIC is considered as non fully Bayesian and does not work for singular models. The WAIC claims to be fully Bayesian since it is based on the entire posterior distribution. WAIC can be calculated as:

WAIC = 1 n n ∑ i=1 log(P(y i |y)) - V n
where, V = ∑ n i=1 Var(log(P(y i |y)). The package "loo" developed in R allows fast computation of the WAIC. The function compares 2 models by estimating the difference between prediction errors. The difference will be positive if the expected predictive accuracy for the second model is higher (for more details, see Watanabe (2010); [START_REF] Vehtari | Efficient implementation of leave-one-out cross-validation and waic for evaluating fitted bayesian models[END_REF]).

• Bayes Factor (BF) can be interpreted as how much the data favor model M 1 over M 2 . It is equal to the ratio of the evidence (prior mean of the likelihood) of one particular hypothesis to the evidence of another:

BF = P(y|M 1 ) P(y|M 2 ) = θ P(y|θ , M 1 )P(θ |M 1 ) θ P(y|θ , M 2 )P(θ |M 2 )
Its most common interpretation is the one --first proposed by [START_REF] Jeffreys | The theory of probability[END_REF] • k-fold Cross Validation (C.V) estimates the predictive power of the model: for example, to be able to use a Bayesian model to predict the soil carbon dynamics for sites where no measurements were taken. This procedure considers splitting the data into k groups of sub-sets. The model is adjusted on k -1 sub-sets "y (training) " and a cost function will be evaluated on the remaining sub-set "y (validation) ". This procedure of training and validation is repeated k times and the k performances are averaged in order to assess the predictive

Conclusion

performance of the model, avoiding overfitting [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF]. k is usually taken as 5 or 10, but there is no general formal rule. As k gets larger, the difference in size between the training set and the re-sampling subsets diminishes. As this difference decreases, the bias of the technique becomes smaller. The relative error is defined as follows:

RE(%) = |y (training) -y (validation) | |y (test,real) | * 100
When multiple models are compared, the model with smaller Relative Error (RE) is preferred to models with larger RE.

• Coverage of Bayesian credible intervals percentage on validation sets

The credible interval covers a pre-specified credibility range of a given unknown. Since the validation sets are not taken into account to adjust the model, one can check that the probability that the actual observation will lie within the credible interval corresponds to its theoretical value. The percentage of coverage of the Bayesian credible interval is defined as follows:

p_cov(%) = 1 L Y (validation) > born_in f & Y (validation) < born_sup * 100
where L represents the number of MCMC iterations, P(Y (validation) < born_sup) = 0.975 and P(Y (validation) < born_in f ) = 0.025. Y (validation) indicates the sub sets of data which are not used to adjust the model. born_sup and born_in f correspond to the 2.5% and 97.5% quantiles.

C-Convergence of Bayesian computations

In addition, an essential point to check is the Monte Carlo convergence of the estimated unknown model parameters. Indeed, in theory, the posterior distribution will converge to the target distribution after an infinite number of MCMC iterations. Practically, the number of MCMC iterations is of course finite and does matter for computational cost and efficiency. Assessing MCMC convergence is crucial since Bayesian statistical reporting is entirely based on the posterior distribution : controlling the numerical accuracy of parameter estimation is therefore necessary to quantify uncertainty. The best way to check the convergence of the Bayesian model is to run multiple Markov chains initialized from different initial conditions. Theoretically, these Markov chains should forget their starting points and converge to the same target distribution. There exist several warning signals which indicate that convergence is not established. One can distinguish the following visual inspections:

• A low or high acceptance rate for the Metropolis Hastings simulation algorithm,

• Poor mixing when observing the trace-plots of the estimated model parameters,

• High autocorrelation between states of the Markov chain. The higher the autocorrelation in the chain, the larger the MCMC variance and the worse the approximation. The lag-t autocorrelation function of the sequence that are t steps apart is defined as [START_REF] Hoff | A first course in Bayesian statistical methods[END_REF], chapter 6):

ac f t (θ ) = 1 S-t ∑ S s=1 (θ s -θ )(θ s+t -θ ) 1 S-1 ∑ S s=1 (θ s -θ ) 2
A high value of the acf indicates that the MCMC is only slowly moving around the parameter space and may take a long time to explore the parameter space adequately.

• Suspicious tails or shapes when examining the posterior distributions for unknown parameters.

There also exist several formal tests to judge MCMC convergence. The most popular are : the Geweke diagnostic [START_REF] Geweke | Evaluating the accuracy of sampling-based approaches to the calculations of posterior moments[END_REF], the Raftery and Lewis diagnostic [START_REF] Raftery | practical markov chain monte carlo]: comment: one long run with diagnostics: implementation strategies for markov chain monte carlo[END_REF], the Gelman and Rubin diagnostic [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF], etc. The Geweke diagnostic takes two non overlapping parts (usually the first 0.1 and last 0.5 proportions) of the Markov chain and compares the means of the two parts, using a difference of means test to see if the two parts of the chain may stem from the same distribution (as a null hypothesis). The disadvantage of this diagnostic is that it is sensitive to the specification of the spectral window. In addition, Geweke does not suggest a quantitative rule to conclude to convergence. On the other hand, the Raftery and Lewis diagnostic estimates the minimum chain length needed to estimate a percentile to some precision. One needs to select a posterior quantile of interest q, an acceptable tolerance r for this quantile and a probability s, which is the desired probability of being within (qr, q + r). Thus the minimum length is given by applying the following recipe:

n min = φ -1 ( s + 1 2 ) q(1 -q) r 2
where φ -1 (.) is the inverse of the normal cumulative distribution function.

The convergence is based on the dependence factor I obtained by the raftery.diag() function in R. A high dependence factor (for example > 5) may be explained by bad starting values, poor mixing or a high correlation between parameters. A review paper by [START_REF] Cowles | Markov chain monte carlo convergence diagnostics: a comparative review[END_REF] pointed out some weaknesses of the Raftery and Lewis convergence diagnostic. For example, variable estimates can be produced given different initial chains starting points. Added to that, it is not realistic to impose that the convergence should be tested for every quantile of interest. where V ar(θ ) = (1 -1 n )W + 1 n B (n is the number of discarded iterations), W is the mean of the variances of each chain, defined as W = 1 m ∑ m j=1 s 2 j with s 2 j the variance of the jth chain given by: s 2 j = 1 n-1 ∑ n i=1 (θ i j -θ j ) 2 . The Variance between chain is given by B = n m-1 ∑ m j=1 ( θ j -θ ) 2 where θ = 1 m ∑ m j=1 θ j and m is the number of total iterations.

The convergence is satisfied when the potential scale reduction factor R for each of the unknown parameters goes to 1. A high value of R (greater than 1.2, according to the recommendations of the authors) underlines the need to run the Markov chain longer so as to reach convergence to the stationary distribution.

One of the challenges of the Gelman and Rubin convergence test is to propose initial values from a starting distribution that is sufficiently dispersed with respect to the target distribution to adequately explore the parameter space ...but not too far from the target because we wish the proposals of the MCMC auxiliary distribution to be finally accepted after a reasonable time so that all the chains start evolving after initialization.

Conclusion 7.3 Effect of the unbalanced experimental design on the quality of the estimators

To carry out an experimental design, we have to clarify the following points:

1. The definition of the purpose and conditions of the experiment.

2. The definition of the factor or factors to study, and its or their levels.

3. The definition of individuals or, more generally, the experimental units that we propose to observe.

4. The definition of the number of observations. 5. How to assign the different objects to the different experimental units.

We consider a two-way analysis of variance (ANOVA) with two additive factors A and B with I and J levels respectively and without considering the interaction effect between the factor A and B. Thus, the m-th observed response when factor A is of level i and factor B is of level j is written as : y i, j,m = µ i, j + ε i, j,m ε i, j,m ∼ N(0, σ 2 )

The mean µ i, j is expressed as:

µ i, j = µ + a i + b j
• µ i, j : is the mean of the observation for the ith level of factor A and the jth level of B.

• a i : is the ith effect of level i within the factor A.

• b i : is the jth effect j within the factor B.

• µ: is the grand mean of observations.

• n: is the total number of observations. The 2-factor experimental design (A and B) is orthogonal if it verifies: n i, j = n i+ × n + j n ∀i = 1, . . . , I and ∀ j = 1, . . . , J

where n i+ is the total number of observations of level "i" within the factor A (respectively n + j is the total number of observations of level "j" within the factor B). This definition means that if an experimental unit (observation) is selected at random, the events "observation of i level of factor A" and "observation of j level of factor B" are independent.

In multiple regression, the experimental design is orthogonal if the explanatory variables are not correlated. The orthogonality of experimental design prevents any confusion of regression effects in the model. However, this case is infrequent since the explanatory variables are not controlled (except in the case where the variables are set by the experimenter).

The presence of categorical predictors in the linear model requires the definition of new variables called: dummy variables. To handle the presence of categorical predictors and the redundancy of the model, two solutions can be considered: 1 -elimination of the intercept µ from the ANOVA model, 2-use of a contrast matrix (the constant term is part of the model).

The most common contrast matrices are: treatment contrast (the most widely used by statisticians), sum contrast and Helmert contrast. The interpretation of the regression effects depends on the type of contrast considered. For example, the treatment contrast is recommended when comparing the level effects to a reference level.

What is the effect of an unbalanced experimental design on the quality of the estimation?

In the following part of this appendix, we will provide an answer for these questions on an illustrated example.

To achieve our purpose and for simplicity, we consider a one-way analysis variance of response Y∈ R n with factor "F" characterized by four levels of size n i for i = 1, 2, 3, 4 respectively. A general expression of the model is:

Y = Xβ + E E ∼ N(0, σ 2 )
Here, the design matrix X can be written in different ways:

• If we decide to remove the intercept, the design matrix X is written by using the binary coding X * . In this case, dummy variables are created and the m-th observation of level j is associated X * [m, ] = (0, . . . , 1 j , . . . , 0), where the only value of 1 corresponds to the jth column of matrix X*.

• If we decide to use a contrast matrix, the design matrix X is written as X = (1 X * C), where 1 is the unit vector and C is the contrast matrix. For example, the matrices corresponding to treatment C t and sum C s contrasts for a factor of 4 levels are written respectively as:

C t =       0 0 0 1 0 0 0 1 0 0 0 1       , C s =       1 0 0 0 1 0 0 0 1 -1 -1 -1      
After constructing the design matrix X, there are many criteria that can be used to optimize the experimental design. Here, we choose to define n 1 , n 2 , n 3 and n 4 , the number of observations that should be collected for each level of factor F, by minimizing the sum of the variance for regression effects under the following constraint:

∑ 4 i=1 n i = n.
The variance of the regression effect β is given by:

V (β ) = (X ′ X) -1 σ 2
Thus, minimizing the sum of the variance of the regression effects consists in minimizing the trace of the matrix (X ′ X) -1 . Here, the Lagrange multiplier is used as a strategy to find the minima of the trace of the matrix (X ′ X) -1 . Thus, under the so-called A optimal design:

• For an ANOVA model without intercept, for a total number of observations n, the number of observations must be equiprobable between the 4 levels:

n 1 = n 2 = n 3 = n 4 = n 4
• For an ANOVA model using contrasts:

-For a treatment contrast, where the first level is taken as the default choice in R to be the baseline reference, the design is optimal if: n 2 = n 3 = n 4 = n 1 2 , 2 5 of the total number of observations n is collected for the baseline level and each of the remaining three levels should be given 1 5 of n.

Conclusion

-For a sum contrast, the design is optimal if the number of observations collected for the baseline level is equal to 1 1+3 √

3 n and each of the remaining three levels should have

√ 3 1+3 √

n observations

How is the regression effect variance affected by the experimental design?

In order to illustrate the impact of an unbalanced experimental design on the quality of the estimation, we summarize in Tables 7.2, 7.3 and7 According to Tables 7.2, 7.3 and 7.4, an unbalanced experimental design with regards to the optimal one decreases the precision of the regression effects. For uncontrolled experiments (observational data), the design can be far away from the optimal design and the scientist may suffer from the poor precision of the regression effects. -120.37 37.46 (Baisden et al., 2002)

1

  Le cycle global du carbone. Les flèches rouges représentent les flux anthropiques et celles représentent les flux naturels (Ciais et al., 2014). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2Stock de carbone total (kg C m -2 ) sur la couche supérieure (0,1 m) par latitude et pour les différentes bases de données[START_REF] Tifafi | Large differences in global and regional total soil carbon stock estimates based on SoilGrids, HWSD, and NCSCD: Intercomparison and evaluation based on field data from Usa, England, Wales, and France[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 La quantité du 14 C de l'année x mesurée en 2018, pour tout matériau atmosphérique. Le sol reproduit la variation du radiocarbone atmosphérique dans le cadre du réservoir de radiocarbone mais, en mélangeant des matières organiques de plusieurs périodes, le signal atmosphérique est fortement atténué dans le sol. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Localisation des sites étudiés (points bleus). La distribution spatiale des sites est plutôt hétérogène,avec un manque évident de données dans les zones extra-tropicales, en Afrique et en Russie. . . . 5 Superposition des sites sur la carte de classification climatique Köppen-Geiger. La "classification climatique de Köppen" est fondée sur l'analyse climatique de la période allant de 1951 à 2000. Elle divise les climats en cinq groupes climatiques principaux, chaque groupe étant divisé en fonction des précipitations saisonnières et des régimes de température. Les cinq groupes principaux sont A (équatorial), B (aride), C (tempéré chaud), D (neige) et E (polaire). Tous les climats, à l'exception de ceux du groupe E, sont ensuite déclinés en sous-groupés associés à des régimes saisonniers spécifiques de précipitations (représenté par la deuxième lettre). Par exemple, Af indique un climat de forêt tropicale humide. Le système assigne un sous-groupe de température pour tous les groupes, sauf ceux du groupe A, indiqué par une troisième lettre pour les climats B, C et D et une seconde lettre pour les climats E. Par exemple, Cfb indique un climat océanique avec des étés chauds comme indiqué par la dernière lettre b. Les climats sont classés en fonction de critères spécifiques propres à chaque type de climat. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Profil statistique du F 14 C en fonction de la profondeur. . . . . . . . . . . . . . . . . . . . . . . . 7 Probabilités d'inclusion a posteriori pour toutes les variables explicatives obtenues en appliquant la SSVS à aux profils F 14 C de la base de données. La taille des points dépend de l'importance de la probabilité d'inclusion a posteriori. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 The global carbon cycle. Red arrows represent anthropogenic fluxes and black arrows represent natural fluxes in PgC yr -1 (Ciais et al., 2014). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The total carbon stock (kg C m -2 ) on the (0, 1 m) upper layer per latitude and for the different databases (Tifafi et al., 2018). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Amount of F 14 C of the year x measured in 2018, for any atmospheric material. The soil reproduces the variation in atmospheric radiocarbon as part of the radiocarbon reservoir but, by mixing organic material from several periods, the atmospheric signal is greatly smoothed in soil. . . . . . . . . . 2.1 Geographical locations of studied sites are represented in blue dots. The spatial distribution of sites is very heterogeneous, with a lack of data in extra-tropical zones, Africa and Russia. . . . . . . . . 10 LIST OF FIGURES 2.2 Correlation matrix corresponding to the original climatic numerical covariates for soil F 14 C profiles.

2. 3

 3 Correlation matrix after removing the strongly correlated covariates based on the F 14 C profiles. . 48 2.4 Graphical distribution of the 131 soil radiocarbon profiles before (left panels) and after (right panels) the merging of original ecosystem and soil type groups. . . . . . . . . . . . . . . . . . . . . . . . 50 2.5 Superimposition of site locations on the Köppen-Geiger climate classification. The "Köppen climate classification" is based on the 50-year period 1951-2000. It divides climates into five main climate groups, with each group being divided based on seasonal precipitation and temperature patterns. The five main groups are A (equatorial), B (arid), C (warm temperate), D (snow), and E (polar). All climates except for those in the E group are assigned to a seasonal precipitation subgroup (represented by a second letter). For example, Af indicates a tropical rainforest climate. The system assigns a temperature subgroup for all groups other than those in the A group, indicated by a third letter for climates in B, C, and D, and a second letter for climates in E. For example, Cfb indicates an oceanic climate with warm summers as indicated by the ending b. Climates are classified based on specific criteria unique to each climate type. . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.6 The cumulative number of F 14 C profiles for the following six numerical climatic covariates: mean annual temperature (panel a), seasonal shift of temperature (panel b), minimal temperature (panel c), mean annual precipitation (panel d), minimal precipitation (panel e) and aridity index (panel f). 53 2.7 Distribution of profiles according to the sampling year (grouped by decade). The majority of the F 14 C profiles were sampled between 1990 and 2000. Only two profiles were collected before the "bomb peak". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.8 Distribution of F 14 C profiles according to the number of measurements. Most sites have between 3 and 10 measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

2. 11

 11 The variation in topsoil carbon content (g/kg) recorded among the different soil types (FAO classification) according to the SoilGrids database (http://www.isric.org/explore/soilgrids). . . . . 59 5.1 Statistical profile of soil F 14 C versus depth obtained from Equation 5.1. . . . . . . . . . . . . . . 94 5.2 The climatic and environmental factors that potentially impact soil radiocarbon dynamics. . . . . 95 5.3 The left panel shows the real variation of the carbon content with depth for all studied sites. The right panel underlines the structure of the statistical model proposed according to the left panel. ω 1 defines the deep soil carbon content, ω 2 represents the topsoil soil carbon content and ω 3 is related to the point from which the curve changes decay rate. . . . . . . . . . . . . . . . . . . . . . . . . 96 LIST OF FIGURES 5.4 Comparison of the same F 14 C profile with different values of layer thickness (2cm-step F 14 C measurement to 36 cm depth, then a 4 cm-step (first panel), 8 and 12 cm (second panel) and 48 cm-step sampling(third panel)) where the sample was collected in 1986 by Becker-Heidmann et al. (2002) under a cultivated field (Vertisol). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Residuals versus fitted values for all F 14 C profiles and all depth measurements. Positive values for the residual mean that the prediction is too low, and negative values mean that the prediction is too high; 0 means that the estimation was exactly correct. . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Beta distribution density for various choices of the two hyperparameters. Beta(2,2) used as a prior distribution for inclusion probability is highlighted in orange. . . . . . . . . . . . . . . . . . . . .

F

  14 C incorporation (φ 3 ) are represented by pie charts. The actual variation of the topsoil and deep radiocarbon according to the ecosystem type are illustrated by boxplots. The boxplot or piechart parts of ecosystem categories having the same color belong to the same cluster. . . . . . . . . . .5.9The fusion of levels for soil type for the deep radiocarbon (φ 1 ), topsoil F 14 C (φ 2 ) and the F 14 C incorporation (φ 3 ) are represented by pie charts. The real variation of the topsoil and deep radiocarbon according to the soil type are illustrated by boxplots. The boxplot or piechart parts of soil categories having the same color belong to the same cluster. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Residuals versus fitted values for soil carbon content model. Positive values for the residual mean the prediction was too low, and negative values mean the prediction was too high; 0 means the estimation was exactly correct. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  and after 2000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Topsoil (depth ≤ 10 cm, panel a) and deep (depth > 100 cm, panel b) F 14 C from the database versus the Mean Annual Temperature (MAT) for the database radiocarbon profiles. . . . . . . . . LIST OF FIGURES 5.16 The impact of the increment of 1 • C in the Mean Annual Temperature (MAT) on the deep F 14 C (φ 1 ) and the F 14 C incorporation (φ 3 ). The profile is plotted for φ 1 = 0.3, φ 2 = 1.1, φ 3 = 20 and φ 4 = 2 . The modified profile is given in the blue curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 5.17 Histogram of the posterior distribution of the Mean Annual Precipitation (MAP) effect on the topsoil radiocarbon (φ 2 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.18 Distribution of the topsoil radiocarbon (less than 10 cm depth) F 14 C from the database versus the Mean Annual Precipitation (MAP). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.19 The impact of an increment of 1 mm in the Mean Annual Precipitation (MAP) on the topsoil F 14 C (φ 2 ) and the F 14 C incorporation (φ 3 ). The profile is plotted for φ 1 (deep F 14 C) = 0.3, φ 2 (topsoil F 14 C) = 1.1, φ 3 (the distance corresponds to half of the topsoil F 14 C amount) = 20 and φ 4 = 2. The modified profile is shown in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

5. 22

 22 Distribution of deep F 14 C (at least 100 cm of depth) according to the soil type. The boxplots of soil categories having the same color belong to the same cluster. . . . . . . . . . . . . . . . . . . . . 134 5.23 Distribution of the topsoil F 14 C (less than 10 cm of depth) according to the soil type. The boxplots of soil categories having the same color belong to the same cluster. . . . . . . . . . . . . . . . . . 136 5.24 The fusion of levels for soil type for the F 14 C incorporation depth (φ 3 ) according to the Bayesian Effect Method with k = 10. Pie slices of soil categories having the same color belong to the same cluster. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 5.25 Distribution of both topsoil (panel a) and deep soil carbon content (panel b) (depth ≤ 10 cm and depth > 100 cm) according to the soil type is illustrated by boxplots. The boxplots of soil categories having the same color belong to the same cluster. . . . . . . . . . . . . . . . . . . . . . . . . . . 137 5.26 The results of the soil type levels fusion for the carbon content incorporation depth (ω 3 )) obtained by the BEF. The same fusion of levels was identified under k = 10, 50 and 100. Pie slices of soil categories having the same color belong to the same cluster. . . . . . . . . . . . . . . . . . . . . 138 5.27 Distribution of both topsoil F 14 C (panel a) and deep F 14 C (panel b) (depth lower than 10 cm and depth higher than 100 cm, respectively) according to the ecosystem type. The boxplots of soil categories having the same color belong to the same cluster. . . . . . . . . . . . . . . . . . . . . 139 5.28 The fusion of levels of ecosystem type for the F 14 C incorporation depth (φ 3 ). The slices of ecosystem categories having the same color belong to the same cluster. The fusion of ecosystem levels is the same for k = 10 and 50. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 5.29 Distribution of both topsoil (panel a) and deep soil carbon content (panel b) (depth ≤ 10 cm and depth > 100 cm) according to the ecosystem type is illustrated by boxplots. The boxplots of soil categories having the same color belong to the same cluster . . . . . . . . . . . . . . . . . . . . . 140 5.30 The fusion of levels of ecosystem type for the carbon content incorporation depth (ω 3 ). The slices of ecosystem categories having the same color belong to the same cluster. The same levels fusion of ecosystem type is obtained under different values of k : k = 10, 50 and 100. . . . . . . . . . . . . 141

5. 32

 32 Soil radiocarbon and carbon content profiles of each of the nine Brazilian natural forest profiles under rather equatorial climate before (continuous green line) and after the conversion into cropland (dashed brown line) (panel a). The average profile characterized by 25.6 • C of Mean Annual Temperature, 2057 mm of Mean Annual Precipitation and 1.56 of Aridity Index is given in panel b. 5.33 The soil radiocarbon and soil carbon content dynamics after (green dashed curve) and before (yellow curve) the conversion of field and cultivated grassland to natural forest for the average temperate profile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.34 Influence of climatic variables on estimated average F 14 C illustrated for Köppen-Geiger climate classes summarized in the Table below the panels. The dark line represents the average profile relative to the Köppen-Geiger climate class. The colored band corresponds to the between-site variability of soil profiles with similar Köppen-Geiger climate classes. The gray band corresponds to the confidence within the same site. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.35 Superposition of the average profiles obtained for all the Köppen-Geiger climate sub-groups from the database for F 14 C (panel a) and carbon content (b). . . . . . . . . . . . . . . . . . . . . . . . 5.36 Average profiles of F 14 C and soil carbon content under the current Mean Annual Temperature (MAT, black curve) and predicted for an increment of l • C (red curve), 1.5 • C (green curve) and 2 • C (blue curve) of MAT. Each panel is relative to a sub-group of the Köppen-Geiger climate classification . 6.1 All F 14 C and soil carbon content database profiles. . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Prediction of two F 14 C profiles chosen from the database according to the Gaussian Process statistical modeling of soil radiocarbon dynamics. The orange band presents the uncertainty inter profile, while the gray band represents the uncertainty intra profile. The red points highlight the real observed radiocarbon measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 The soil radiocarbon activity variation against the variation in soil carbon content based on 125 database profiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Percentages show the fraction of the total carbon reservoir of each type. Numbers after the slash

  14 C cycle with time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Development the atmospheric radiocarbon in the Northern Hemisphere in the last 50 years. Data before 1959 were derived from tree rings (Stuiver and Quay, 1981). From 1959 to 1983 measurements were performed at the Alpine site Vermunt; data from 1984 onwards are from the Schauinsland in the Black Forest (http://slideplayer.fr/slide/12319388/) . . . . . . . . . . . . . . . . . . . . . . . 7.5 The New Zealand curve is representative of the Southern Hemisphere; the Austrian curve is representative of the Northern Hemisphere. Atmospheric nuclear weapon tests almost doubled the concentration of 14 C in the Northern Hemisphere (Currie, 2004) . . . . . . . . . . . . . . . . . . LIST OF FIGURES 7.6 Posterior predictive distribution of the global mean temperature increase relative to 1861-1880 in Celsius degrees. This histogram is extracted from Fig. 3 of Raftery et al. (2017) . . . . . . . . . . . . . . . . . . . 7.7 General statement of Bayes theorem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.8 Illustration of the fundamental steps of Bayesian analysis (Healy and Proctor, 2003). . . . . . . . 7.9 Illustration of the Bayesian analysis for the wreckage search of Air France Flight 447 (High likelihood areas given in Fig 13 from Stone et al. (2011)) . . . . . . . . . . . . . . . . . . . . . . 7.10 Hierarchical modeling strategy. Factorization of the complexity and Bayesian inferences. This representation is extracted from Fig. 1.12 of Parent and Rivot (2012) with bracket notations for probability distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2. 8

 8 Number of F 14 C profiles according to the number of measurements for each aggregated ecosystem.5.1 Abbreviations of the climatic and environmental predictors considered for soil radiocarbon modeling. 5.2 Abbreviation of the eight climatic and environmental predictors considered for soil carbon content modeling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Deviance Information Criterion (DIC) comparison between the Full additive model and the multiplicative Full model. The model with the lowest DIC is preferred to models with larger DIC. . . . 16 LIST OF TABLES 5.4 The best subsets of climatic and environmental predictors for latent linear models of φ 1 , φ 2 and φ 3 chosen according to the posterior median estimation of regression effects. . . . . . . . . . . . . . 102 5.5 The best subsets of climatic and environmental predictors for latent layers φ 1 , φ 2 and φ 3 chosen according to the Posterior inclusion Probability (PIP). The significant predictors are detected with a PIP at least equal to 0.5. The predictors highlighted in blue were the ones detected in addition to those identified by the posterior median estimation of regression effects. The predictors are ordered according to the PIP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.6 Bayesian criteria comparison for the posterior median model and the model selected according to the posterior inclusion probabilities of predictors. The model with the lowest Deviance Information Criterion (DIC) is preferred to models with higher DIC. The model with a p-value of the Posterior Predictive Check (P.P.C) close to 0.5 is preferred to models with p-values close to 0 or 1. R.E refers to the Relative Error computed for both learning and validation sets. . . . . . . . . . . . . . . . . 103 5.7 The best subsets of climatic and environmental predictors for latent linear models for φ 1 , φ 2 and φ 3 chosen according to the Posterior inclusion Probability (PIP). The significant predictors are detected with a PIP at least equal to 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.8 The same sub-set of predictors is identified under the three choices of k values. The sub-model identification is based on the Posterior Fusion Probability (PFP) and the Posterior Median Regression Effect (PMRE). The Deviance Information Criterion for model fitting is given in the DIC column. The column named "R" indicates Gelman & Rubin's potential scale reduction factor for model convergence. The check-mark underlines that the convergence has been achieved while the Xmark indicates a poor model convergence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.9 Comparison of the Full Bayesian model and the sub-models identified by the Bayesian selection approaches for soil radiocarbon dynamics using the Bayesian selection criteria. The model with the lowest Deviance Information Criterion (DIC) is preferred to models with a higher DIC. A p-value of the Posterior Predictive Check (P.P.C) close to 0.5 indicates a good model fitting. The model having the smallest Relative Error (R.E) on validation sets has the best predictive power. The posterior coverage of the credible intervals on the validation sets should be around 95%. . . . . . . . . . . . 113 5.10 DIC comparison between the full additive model and the multiplicative full model. The model with the lowest DIC is preferred to models with a larger DIC. The p-value of the Posterior Predictive Check is calculated according to the statistical quantity T (C, ω) = E(Cf (ω, depth)) where f is the non linear deterministic form used to model soil carbon content and C is the soil carbon content response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.11 The best sub-predictors detected according to the two Bayesian selection criteria used for the Bayesian Group Lasso with Spike and Slab prior. The sub-predictors selected according to the posterior median estimation of regression effects are indicated in the column called "Median model". The sub predictors detected according to the Posterior inclusion Probability (PIP) are represented in the column called "PIP model". This column also contains the PIP for the significant predictors. The symbol C in the second column refers to the carbon content. . . . . . . . . . . . . . . . . . . 116 5.12 The Deviance Information Criterion (DIC) for the two sub-models detected by the Bayesian Group Lasso with Spike and Slab. The Median model is based on the posterior median estimation of regression effects while the PIP model is based on the Posterior Inclusion Probability (PIP) for predictors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.13 The best sub-predictors for each latent linear model identified according to the Posterior Inclusion Probability (PIP). A predictor is considered as significant if its PIP is at least equal to 0.5. . . . . . 117 5.14 The best sub-set of predictors is identified for each choice of k. The Deviance Information Criterion for model fitting is given in the column DIC. The column named "R" indicates Gelman & Rubin's convergence. The check-mark underlines that convergence has been achieved while the Xmark indicates a difficulty in achieving model convergence. The blue Xmark indicates a poor convergence. 5.15 Comparison of the Full Bayesian model and the sub-models identified by the Bayesian selection approaches for soil carbon content dynamics using the Bayesian selection criteria. The model with the lowest Deviance Information Criterion (DIC) is preferred to models with a higher DIC. A p-value of the Posterior Predictive Check (P.P.C) close to 0.5 indicates a good model fitting. The model having the smallest Relative Error (R.E) on validation sets has the best predictive power. The posterior coverage of the credible intervals on the validation sets should be around 95%. . . . . . 5.16 Amount and type of clay generally observed in the soil types from the database. Types of soil are ranked according to the result of Bayesian Effect Fusion of soil type levels for deep soil radiocarbon activity (φ 1 ). Cluster number and color are the ones used in Figure 5.22. The number of profiles associated to each type of soil and each group of type of soils are provided in brackets. Column six refers to the median value of deep (higher than 100 cm deep) soil F 14 C from the database (line inside rectangle in Figure 5.22). The last column refers to 25% and 75% quantiles (q= quantile) (the upper and lower rectangle bounds in Figure 5.22) * The result of the clay amount and the type of clay corresponds to Arenosol soil type, only. . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.19 The latent values for both F 14 C and soil carbon profile corresponding to current temperature and an increase of the MAT by +1 • C, +1..5 • C and +2.5 • C. φ 1 and ω 1 refer to deep F 14 C and carbon content respectively, φ 2 and ω 2 represent the topsoil F 14 C and soil carbon content and finally φ 3 and ω 3 underline the F 14 C and carbon content incorporation depth. . . . . . . . . . . . . . . . . . 7.1 Interpretation of Bayes Factors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Sum of variances of the regression effects for the constant model under different experimental designs. n 1 refers to the number of observations for the baseline level. The total number of observations is set to 20. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Sum of variances of the regression effects under different experimental designs. A treatment contrast is used to design the matrix design. n 1 refers to the number of observations for the baseline level. The total number of observations is set to 20. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Sum of variances of the regression effects under different experimental designs. A sum contrast is used to design the matrix design. n 1 refers to the number of observations for the baseline level. The total number of observations is set to 20. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . background-corrected and δ 13 C-normalized (the 1950 atmospheric activity) sample activity, background-corrected and δ 13 Csample activity to the standard activity measured in the same year, both activities background-corrected and δ 13 C-normalized Atmospheric F 14 C of the sampling year Gram per La représentation actuelle de la dynamique du carbone du sol par les modèles mécanistes n'est pas entièrement satisfaisante

  La deuxième partie souligne les propositions et possibilités qui peuvent être examinées afin d'améliorer la structure du modèle statistique et d'optimiser les perspectives d'utilisation du modèle statistique. détecté comme significatif pour le F 14 C en surface ni pour l'incorporation du F 14 C en profondeur. Alors, pourquoi ne pas envisager l'analyse d'échantillons d'archives, comme ceux de Rothamsted? Cette collection d'échantillons a été créée par Lowes et Gilbert en 1843. Plusieurs milliers de sols recueillis dans les années 1920-1950 sont stockés dans le référentiel. Environ 1200 cultures et 200 échantillons de sol sont ajoutés chaque année aux archives.

  Carlo (MCMC) simulations and the fourth chapter applies the Bayesian Selection methods to soil radiocarbon and soil carbon content dynamics. Chapter 4 intends to propose possible physical interpretation of the selected climatic and environmental factors that explain the shape of radiocarbon and carbon content profiles. These variables are the ones selected by the best subset by the Bayesian selection. It also provides a synthetic view of profiles shape under different climate. The predictive capacity of the model is also tested under two scenarios of land use change (deforestation in equatorial region, agricultural decline and reforestation in temperate region) and under 3 scenarios of global warming. The manuscript ends by a conclusion of my thesis work and with some propositions and perspectives in order to improve the soil carbon statistical modeling and the database.

  carbon profiles database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Original database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Processed database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Potential explanatory covariates affecting soil carbon dynamics . . . . . . . . . . . . . . 2.1.3.1 Potential climatic numerical predictors . . . . . . . . . . . . . . . . . . . . . . 2.1.3.2 Atmospheric 14 C concentration of the sampling year . . . . . . . . . . . . . . . 2.1.3.3 Potential environmental categorical predictors : ecosystem and soil type . . . . 2.1.4 Final selected climatic and environmental potential predictors for the soil carbon dynamics 2.1.5 Database evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.5.1 World climatic zones and distribution of the profiles . . . . . . . . . . . . . . . 2.1.5.2 Soil type diversity and distribution of profiles . . . . . . . . . . . . . . . . . . 2.1.5.3 Ecosystem diversity and distribution of profiles . . . . . . . . . . . . . . . . . 2.1.5.4 Sampling year and distribution of profiles . . . . . . . . . . . . . . . . . . . . 2.1.5.5 Variation in the number of observations per profile . . . . . . . . . . . . . . . . 2.1.5.6 Large variation of the topsoil carbon content . . . . . . . . . . . . . . . . . . . 2.2 The take-home messages of the database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Finally, Section 5 .

 5 5 will underline the predictive model applications in a context of global climatic and land use changes. Chapter 5 (b) A logarithm transformation is used to model the fourth latent variable, φ 4 .

5. 5

 5 Predictive model applications in a context of global changes 5.5.1 Impact of a land use change 5.5.1.1 Conversion of equatorial forest to cropland impacts both topsoil carbon content and deep carbon dynamics

  as the ones of the Rothamsted archives? This collection of samples was established by Lowes and Gilbert in 1843. Several thousands of soils collected in the 1920s-1950s are stored in the repository. About 1200 crop and 200 soil samples are added annually to the archive.

Figure 7 .

 7 Figure 7.10: Hierarchical modeling strategy. Factorization of the complexity and Bayesian inferences. This representation is extracted from Fig. 1.12 of Parent and Rivot (2012) with bracket notations for probability distributions.

  Finally, the Gelman and Rubin diagnostic compares intra and inter variances of Markov chains. The implementation of the Gelman and Rubin test is available in the programs developed for MCMC simulations such as JAGS/ BUGS/ Stan, etc. The potential scale reduction factor R is defined as:

. 4 ,

 4 the sum of variances of the regression effects under different experimental designs: balanced, unbalanced and strongly unbalanced for n = 20. Experimental design n 1 n 2 n 3 n 4 sum of regression effect

1

  Facteurs climatiques et environnementaux explicatives potentiellement explicatifs de la dynamique du radiocarbone et de la teneur en carbone dans le sol. Le F 14 C atmosphérique de l'année d'échantillonage est mis en évidence en italique gras car il n'affecte que la dynamique du radiocarbone de sol. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Les facteurs climatiques et environnementaux sélectionnés pour les variables latentes du modèle hiérarchique φ 1 , φ 2 , φ 3 et φ 4 . Pour la troisième couche latente φ 3 , toutes les variables explicatives

ont été sélectionnées. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Comparaison de BGL-SS, BAGS et BEF, trois méthodes bayésiennes récentes appropriées aux variables catégorielles: quoi, comment et pour quoi faire? . . . . . . . . . . . . . . . . . . . . . . 2.1 Merging of WRB soil type groups for soil radiocarbon profiles according to expert advice. For ease of reference, we will hereafter use the soil type group "short name" (e.g. Chernozem) to refer to the concatenation of the merged groups (e.g. Chernozem, Kastanozem, Phaeozem). . . . . . . . . . . 2.2 Ecosystem aggregated types for soil radiocarbon profiles according to expert advice. . . . . . . . 2.3 The categories of soil and ecosystem types for the six profiles removed for the soil carbon content statistical modeling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The potential explanatory climatic and environmental covariates that may impact soil radiocarbon carbon dynamics as well as soil carbon content dynamics. Atmospheric radiocarbon is highlighted in bold italics since it only impacts the radiocarbon soil dynamics. . . . . . . . . . . . . . . . . . 2.5 Description of the Köppen-Geiger classification (1st and 2nd letter description only) and number of radiocarbon profiles selected from the database that correspond to the different subgroups (last column). P min (P max ) and T min (T max ) are for the minimum (maximum) monthly precipitation and temperature, P ann is for the MAP, S and W subscripts are for summer and winter respectively. P th = 2 * MAT + a, with a = 0 if at least 2/3 of MAP occurs in winter, a = 28 if at least 2/3 of MAP occurs in summer and a = 14 otherwise. The calculation key implies that the polar climates (E) have to be determined first, followed by the arid climates (B) and subsequent differentiations into the equatorial climates (A) and the warm temperate and snow climates (C) and (D), respectively. 2.6 Percentage of the total continental land area on Earth by soil types (first two columns) and by merged groups of soil type (columns 3 and 4). The last column gives the number of profiles by merged group. Note for Vertisol: the total land area differs according to the classification, values according to the FAO [*] and USDA [ + ] classifications are provided Source : https ://www.britannica.com . 2.7 Number of F 14 C profiles according to the number of measurements for each aggregated soil type.

  Ce chapitre présente donc trois approches bayésiennes de sélection appropriées aux variables catégorielles. Parmi ces approches, le Bayesian Sparse Group Selection (BSGS) et le Bayesian Effect Fusion (BEF) permettent d'aller au-delà de la simple sélection des variables catégorielles. Le BEF peut être utilisé afin de fusionner les modalités ayant le même effet au sein de chaque facteur et le BSGS nous permet d'établir un jugement probabiliste sur l'inclusion des effets des différentes modalités relatives au même groupe. Pour la derniére approche, le Bayesian Group Lasso with Spike and Slab (BGL-SS), l'estimateur a posteriori de la médiane présente une excellente performance, à la fois pour estimer et pour sélectionner les effets de régression. Notre objectif est d'appliquer

	ces techniques de sélection sur la couche latente du modèle non linéaire hiérarchique afin d'identifier les facteurs
	climatiques et environnementaux significatifs.
	Plusieurs méthodes de sélection bayésienne ont été proposées dans la littérature afin de choisir le meilleur
	sous-modèle. On peut distinguer par exemple le Stochastic Search Variable Selection (SSVS) proposé par George
	and McCulloch (1993), la méthode de Kuo and Mallick (1998) ainsi que le Gibbs Variable Selection suggéré par
	Dellaportas and Ntzoufras (1997). Cependant, ces techniques capturent seulement les effets de régression relatifs à
	chaque variable continue et non pas les effets d'un regroupement des modalités associées aux variables catégorielles.
	Ici, on s'intéresse aux méthodes de sélection appropriées aux variables catégorielles qui exigent l'introduction
	de variables fictives (dummy variables) dans le modèle. Considérons d'abord le modèle d'analyse de la variance
	suivant:
											curve shape
	Posterior inclusion probability	0.4 0.5 0.6 0.7		Land use						Soil type		• Aridity Seasonal shift Precipitation 14C atm Temperature Latitude	•
		Cultivated-Filed	Cultivated-Forest	Cultivated-Grassland	Forest	Natural	Natural-Grassland	Natural-Savanna	Andosol	Arenosol	Cambisol	Chernozem	Ferralsol	Fluvisol	Gleysol	Nitisol	Podzol	Vertisol	14Catm	MAT	MAP	Lat	AI	Dif_T
	Figure 7: Probabilités d'inclusion a posteriori pour toutes les variables explicatives obtenues en appliquant la SSVS à aux profils F 14 C de la base de données. La taille des points dépend de l'importance de la probabilité d'inclusion a
	posteriori.
	Exploration de trois méthodes bayésiennes de sélection de variables catégorielles et leurs codages sous
	JAGS
	Le troisième chapitre de la thèse est présenté sous la forme d'un article, soumis le 15 août au journal Bayesian
	Analysis.
	Le SSVS exploré dans le chapitre précédent peut conduire à faux négatifs pour certaines variables explicatives
	catégorielles. Afin de résoudre ce problème, un intérêt particulier a été porté à l'exploitation d'autres méthodes
	bayésiennes de sélection appropriées aux variables catégorielles.

  (Y 1 , . . . ,Y n ) est le vecteur réponse , G est le nombre de variables catégorielles et µ représente la constante associée au vecteur unitaire 1. Chaque variable catégorielle g renferme C g modalités telles que ∑ G l=1 C l = p. β g est le vecteur des effets de taille (C g × 1) relatif au facteur g. X g est la matrice de design de taille (n ×C g ) définie avec un contraste traitement et ε représente l'erreur.Dans cet article, on présente et discute la performance de ces trois techniques de sélection appropriées aux variables catégorielles: Bayesian Sparse Group Selection proposé par[START_REF] Chen | Bayesian sparse group selection[END_REF], Bayesian Group Lasso with Spike and Slab suggéré par[START_REF] Xu | Bayesian variable selection and estimation for group Lasso[END_REF] et le Bayesian Effect Fusion using model-based clustering défini par Malsiner-

	Walli et al. (2017).

Table 3 :

 3 Comparaison de BGL-SS, BAGS et BEF, trois méthodes bayésiennes récentes appropriées aux variables catégorielles: quoi, comment et pour quoi faire? En premier lieu, ce chapitre décrit les structures des modèles statistiques des profils du radiocarbone et des teneurs en carbone indépendamment. Une modification du modèle statistique des profils du F 14 C est mise en oeuvre par rapport au modèle statistique publié dans le journal SFDS. Par exemple, la latitude ne fait plus partie des variables explicatives et une loi normale tronquée est proposée pour modéliser la réponse F 14 C.Les meilleurs sous-modèles identifiés pour les profils de F 14 C et des teneurs en carbone Les méthodes bayésiennes de sélection ont été comparées entre elles en se fondant sur des critères bayésiens de comparaison de modèles comme le DIC, l'erreur relative sur les sites de validation (cross-validation), etc. Les deux meilleurs sous-modèles obtenus après comparaison, sont résumés dans les tableaux suivants:

	Méthodes bayéseinnes de sélection
	BGL-SS BSGS	BEF

Table 2 .

 2 .2).

	ecosystem database profiles nbr short name for the merged ecosystem profiles nbr
	forest natural-forest	7 42	natural forest	49
	cultivated-forest	10	cultivated forest	10
	natural natural-savanna	8 4	others	13
	natural-desert	1
	natural-grassland	33	natural-grassland	33
	cultivated-grassland	8	cultivated-grassland	8
	field	18	field	18

2: Ecosystem aggregated types for soil radiocarbon profiles according to expert advice.

Table 2 .

 2 6: Percentage of the total continental land area on Earth by soil types (first two columns) and by merged groups of soil type (columns 3 and 4). The last column gives the number of profiles by merged group. Note for Vertisol: the total land area differs according to the classification, values according to the FAO [*] and USDA [

	WRB soil type	% of the total merged WRB group % of the total nb of F 14 C profiles land area of soil type land area
	Arenosol	7	Areno-
	Regosol-Arenosol	2	Regosol	22	7
	Leptosol	13
	Fluvisol Cambisol	2.8 12	Cambisol	14.8	16
	Nitisol Ferralsol	1.6 6	Ferralsol	7.6	18
	Gleysol	5.7
	Planosol	1.0	Gleysol	7.2	9
	Plinthosol	0.5
	Chernozem	1.8
	Kastanozem	3.7	Chernozem	7	19
	Phaeozem	1.5
	Luvisol	5	Luvisol	5	27
	Podzol	4	Podzol	4	16
	Vertisol	2.7(*)-2( + )	Vertisol	2.7(*)-2( + )	7
	Andosol	1	Andosol	1	12

+ ] classifications are provided Source : https ://www.britannica.com

2.1.5.3 Ecosystem diversity and distribution of profiles

The database offers a correct representativeness of ecosystems (Table

2

.2). Each of the 6 selected categories is represented by 8 to 49 profiles, i.e. by 6 to 37% of the available profiles. However, it remains unbalanced with a high weight (45% of the selected database) of forests, mostly natural forests, whereas grasslands account for 31% leaving the remaining 24% for both fields and other types of ecosystems. This chapter introduces a published article on the 13 th of September 2018 in the "Journal de la Société Française de Statistique".

TABLE 4 .

 4 SSVS evaluation for artificial data including only independent quantitative covariates. Panels a, b, c, d are the results obtained for φ 1 , φ 2 , φ 3 and φ 4 latent layers, respectively. Rows give the three most visited sub-models. Columns correspond to the different tested priors. F.C. (False Choice) sums up both false inclusion and false exclusion. Prob. is the probability appearance of model subsets throughout iterations. The best sub-models detected by the SSVS with the three proposed values of c do not contain any false detection.

		c = 10	c = 100	c = 1000	c = 5000
		Prob F.C Prob F.C Prob F.C Prob F.C
	a)	0.31	-	0.66	-	0.81	-	0.84	-
		0.13 0.12	X 3 X 5	0.09 0.08	X 3 X 1	0.06 0.05	X 1 X 3	0.06 0.04	X 3 X 5
		c = 10	c = 100	c = 1000	c = 5000
		Prob F.C Prob F.C Prob F.C Prob F.C
	b)	0.30	-	0.46	-	0.40	-	0.39	-
		0.10 0.08	X 2 X 1	0.17 0.06	X 2 X 6	0.25 0.08	X 2 X 6	0.28 0.07	X 2 X 6
		c = 10	c = 100	c = 1000	c = 5000
		Prob	F.C	Prob F.C Prob F.C Prob F.C
	c)	0.33	-	0.61	-	0.72	-	0.74	-
		0.23 X 3 , X 5 0.21 X 3	0.16 0.14	X 5 X 3	0.13 0.10	X 5 X 3	0.12 0.11	X 5 X 3
		c = 10	c = 100	c = 1000	c = 5000
		Prob	F.C	Prob F.C Prob F.C Prob F.C
	d)	0.29	-	0.53	-	0.58	-	0.59	-
		0.25 X 2 , X 3 0.21 X 2	0.18 0.17	X 2 X 3	0.17 0.16	X 2 X 3	0.18 0.15	X 2 X 3

TABLE 9 .

 9 High 2 frequency models (Model1 and Model2) for each of the latent linear models. It represents the 2 most frequent combinations of explanatory covariates among all the MCMC iterations. The linear models with all explanatory covariates are identified for φ 2 , φ 3 and φ 4 .

	Latent linear model	High frequency model	frequency (n.iter = 180,000)
	φ 1	Model1: land use, temperature and seasonal shift	12,549
		Model2: land use, seasonal shift	10,822
	φ 2	Model1: all covariates	6,606
		Model2: seasonal shift	4,272
	φ 3	Model1: all covariates	36,819
		Model2: all covariates except land use	12,587
	φ 4	Model1: all covariates	14,782
		Model2: land use, F 14 C atmospheric, latitude	7,336

Table 5

 5 .6).

	Models	DIC	p-value of P.P.C	R.E on validation sets R.E on learning sets 5-fold C.V(%) 5-fold C.V	Posterior coverage on validation sets(%)
	Median model -2340	0.568	36.22	11.89	95.5
	PIP model	-2356	0.566	36.87	11.62	95.1

.6: Bayesian criteria comparison for the posterior median model and the model selected according to the posterior inclusion probabilities of predictors. The model with the lowest Deviance Information Criterion (DIC) is preferred to models with higher DIC. The model with a p-value of the Posterior Predictive Check (P.P.C) close to 0.5 is preferred to models with p-values close to 0 or 1. R.E refers to the Relative Error computed for both learning and validation sets.

Table 5

 5

	.7:

Table 5 .

 5 8: The same sub-set of predictors is identified under the three choices of k values. The sub-model identification is based on the Posterior Fusion Probability (PFP) and the Posterior Median Regression Effect (PMRE). The Deviance Information Criterion for model fitting is given in the DIC column. The column named "R" indicates Gelman & Rubin's potential scale reduction factor for model convergence. The check-mark underlines that the convergence has been achieved while the Xmark indicates a poor model convergence.

	latent variables	k = 10 selected predictors	DIC	k = 50 R selected predictors	DIC	k = 100 R selected predictors	DIC	R
		MAT		MAT		MAT
	φ 1	Soil Land		Soil Land		Soil Land
	φ 2	Land Soil	-2354	Land Soil	-2363	Land Soil	-2336 ✗
		MAP		MAP		MAP
	φ 3	AI Land		AI Land		AI Land
		Soil		Soil		Soil

  5. The BGL-SS selection results, in latent linear models, for both Bayesian selection criteria are summarized in Table 5.11:

	Latent variables Physical interpretation	Median model best sub-predictors best sub-predictors PIP PIP model
			Land	Land	71
	ω 1	Deep C	Soil MAP	Soil MAP	71 63
				AI	50
			Soil	Soil	100
	ω 2	Topsoil C	MAT Dif_T	MAT Dif_T	100 70
			Soil	Soil	100
	ω 3	C incorporation	Dif_T MAP	Dif_T map	62 56

  • C, the impact does not propagate in extent for +1.5 • C and +2 • C. The +1.5 • C and +2 • C profiles are superimposed on the +1 • C profile.

• The impact of warming differs greatly from one sub-group to another. The profiles that evolve the most with the global warming are under Equatorial savannah climate with dry winter, Steppe climate, Warm temperate dry climate with dry summer, Snow climate with dry winter or summer. This means that about half of the climate sub group does not show a noticeable change in carbon profiles associated to global warming.

  It was quickly recognized that yr BP were not equivalent to calendar years before 1950. It is necessary to use the calibration (e.g. IntCal13,[START_REF] Reimer | Intcal13 and marine13 radiocarbon age calibration curves 0-50,000 years cal bp[END_REF]) to transform a yr BP into a real age. The real age is output as cal.yr BP. As already stated, the calibration takes into account the real period of 5730 yrs and the variations in both sources and sinks of 14 C with time.

	2.
	Furthermore, it was calculated according to the Libby half-life of 5568 years.
	t = -5568/ln(2) * ln(A SN /A ON )

  ON A ON : normalized sample activity (corrected for radioactive decay to 1950)The Reverend Thomas Bayes was born in England in 1701 or 1702. The son of a Presbyterian minister from whom he took over, he was a nonconformist intellectual who seems to have lived a peaceful life as a studious bachelor, passionate about theology, natural sciences, mechanics and mathematics. His name is now associated with an elementary mathematical formula taught in any introductory course to Probability and Statistics, but Thomas Bayes would not have received the fame he now enjoys among all Bayesian statisticians if, two years after his death, a friend of his (Richard Price) had not searched his archives for Bayes' famous posthumous text -An Essay towards solving a Problem in the Doctrine of Chances -published after a reading before the Royal Society in 1763.

	7.2 Bayesian modeling and inference
	Bayes and Laplace
	This unit will be adopted for our statistical analysis since it best represents the value that is closest to raw
	radiocarbon activity.
	Conclusion

  (Numerical Inference for Statistical Models using Bayesian and Likelihood Estimation) is an interesting R package developed by Perry De Valpine, Christopher Paciorek, Duncan Temple, and Daniel Turek whose first version was published in June 2014 . The package (de Valpine et al., 2017

[START_REF] King | Bayesian Analysis for Population Ecology[END_REF]

is of special interest for the practitioner because, for each case study, they provide their inference routines both in R (as a stand-alone program) and in WinBUGS. NIMBLE

  and slightly modified by Lee et al. (2014)-given in Table7.1:

	Bayes factor	Evidence category
	≥ 100 30-100 10-30 3-10 1-3 1	Extreme evidence for M 1 Very strong evidence for M 1 Strong evidence for M 1 Moderate evidence for M 1 Anecdotal evidence for M 1 No evidence
	1/3-1 1/10-1/3 1/30-1/10 1/100-1/30 ≤ 1/100	Anecdotal evidence for M 2 Moderate evidence for M 2 Strong evidence for M 2 Very strong evidence for M 2 Extreme evidence for M 2
	Table 7.1: Interpretation of Bayes Factors.

Table 7 .

 7 2: Sum of variances of the regression effects for the constant model under different experimental designs. n 1 refers to the number of observations for the baseline level. The total number of observations is set to 20. Experimental design n 1 n 2 n 3 n 4 sum of regression effect variances Table 7.3: Sum of variances of the regression effects under different experimental designs. A treatment contrast is used to design the matrix design. n 1 refers to the number of observations for the baseline level. The total number of observations is set to 20. Experimental design n 1 n 2 n 3 n 4 sum of regression effect variances

	Optimal Unbalanced Strongly unbalanced 3 10 5 8 4 4 4 3 6	4 7 2	1.25 σ 2 1.64 σ 2 2.13 σ 2
	Optimal Unbalanced Strongly unbalanced 3 10 5 3 5 6 4 3 6	6 7 2	0.48 σ 2 0.54 σ 2 0.68 σ 2

Table 7 .

 7 4: Sum of variances of the regression effects under different experimental designs. A sum contrast is used to design the matrix design. n 1 refers to the number of observations for the baseline level. The total number of observations is set to 20.

Afin d'améliorer l'expression des incertitudes du modèle statistique, j'ai proposé un modèle de processus gaussien bivarié pour une modélisation conjointe du radiocarbone et de la teneur en carbone. En effet, notre modèle statistique ne prend pas en compte la dépendance entre les différentes mesures au sein d'un même site entre ces deux éléments. De plus, une visualisation de la variation du radiocarbone en fonction de la variation de la teneur en carbone met en évidence la présence d'une corrélation positive entre les deux réponses.En ce qui concerne l'amélioration de la base de données, lors de l'analyse statistique, certains profils sortent nettement hors de la tendance générale. Ces horsains incluent notamment certains profils au Royaume-Uni sous un climat très humide. Ces profils ont probablement biaisé la sélection des variables et très probablement l'ampleur des diverses estimations. Une deuxième série d'évaluations peut être réalisée sans ces profils aberrants pour obtenir un modèle statistique qui mettrait mieux en lumière la tendance générale.Un autre défi consistera à améliorer la base de données actuelle. Comme nous l'avions noté dans le chapitre 2, la base de données est loin d'être parfaite. Par exemple, nous n'avons pas le même nombre de profils par type de sol et d'écosystème. 37% des profils de la base de données sont des forêts, alors que seulement 8% sont définis comme des prairies cultivées. Même pour le type de sol, 9% des profils de la base de données sont définis en tant qu'Andosols et 7% en tant que Régosols / Arénosol / Leptosol. Cependant, seulement 1% de la superficie terrestre continentale de la planète est occupée par des Andosols et 22% par des Regosols / Arenosol / Leptosol. De plus, étant donné que le type de sol et l'écosystème sont associés, il pourrait également être intéressant de diviser Gleysol en deux catégories: les gleysols tropicaux et les gleysols boréaux. Certaines de leurs caractéristiques sont similaires car elles portent le même nom, mais d'autres, en particulier les interactions avec la végétation, sont différentes. L'augmentation du nombre de profils par catégorie de sol et d'écosystème (utilisation des sols + végétation) n'est toutefois pas le seul point important; il faut également tenir compte de plusieurs caractéristiques cruciales, telles que: l'occupation totale de la surface terrestre continentale par le sol et la couverture végétale, et l'association entre le sol, la couverture végétale et les conditions climatiques. De plus, la base de données actuelle ne contient pas de profils de classes climatiques arides et hyper-arides. C'est un manque crucial, en particulier pour l'utilisation du modèle en mode de prévision dans ces régions particulièrement vulnérables aux changements climatiques. Comme indiqué dans la description de la base de données, la répartition non homogène des années d'échantillonnage dans la base de données empêche une bonne représentation du profil F 14 C, notamment pour la variable latente qui donne la profondeur d'incorporation.. Un autre point à améliorer dans la base de données est donc la distribution des profils par année d'échantillonnage. Dans la base de données actuelle, 53% des profils sont échantillonnés 1990 et 2000. En conséquence, le radiocarbone atmosphérique lié aux essais nucléaires au début des années 1960 n' a pas été

site véritable. Ainsi, de nouvelles variables latentes sont générées pour les modèles statistiques des profils F 14 C et de la teneur en carbone. Ces considérations suscitent la question suivante: comment prédire le profil F 14 C et de la teneur en carbone en fonction de l'évolution des conditions climatiques ou environnementales, en tenant compte des mesures déjà observées pour le site correspondant? Nous pouvons ajouter aux variables latentes actuelles et estimées, dans les conditions climatiques et environnementales actuelles des sites, le changement d'effet résultant du remplacement d'une forêt par une terre cultivée ou de l'augmentation de la température de 1 • C. En revanche, suivre cette proposition ne garantit pas la contrainte de positivité des variables latentes du modèle.En second lieu, on peut se demander comment extrapoler le modèle statistique bayésien pour avoir un profil prédit de la teneur en carbone et du radiocarbone au niveau régional ou mondial. Avec une base de données plus complète et un dispositif expérimental plus équilibré, d'autres portes s'ouvrent. Il devient alors possible d'appliquer le modèle statistique conçu sur l'ensemble de la base de données à des modèles pour chaque type de sol et chaque zone climatique. Cela augmenterait la puissance de projection de l'étude. Cela permettrait de mieux déchiffrer l'impact du changement d'affectation des sols en fonction du type de sol et de mieux prédire l'impact du réchauffement climatique actuel selon les régions du monde. Se relier à un système d'information géographique (SIG) est également possible. On parle ici d'extrapolation à 3 dimensions: longitude, latitude et profondeur. La cartographie numérique des sols (DSM) utilisant des modèles spatiaux d'informations contextuelles de l'apprentissage profond (deep learning), est très populaire, et a déjà été utilisée pour générer des cartes[START_REF] Mcbratney | On digital soil mapping[END_REF]. En effet, il existe des méthodes d'apprentissage approfondi, telles que les réseaux de neurones à convolution, qui développent l'approche DSM classique en incluant des informations sur la proximité d'un site. Chaque site est caractérisé par des covariables climatiques et environnementales avec une matrice tridimensionnelle pour la largeur, la longueur en pixels d'une fenêtre centrée en un point (coordonnées du site) et en connaissant les covariables. L'apprentissage multitâche peut gérer la notion de profondeur en fournissant des prédictions, couche par couche. La possibilité d'extrapoler le modèle statistique serait très utile car l'optimisation de la conception de l'échantillonnage prend beaucoup de temps et est également coûteuse (acquisition de données et traitement des échantillons en laboratoire).En outre, les approches de sélection bayésiennes peuvent aider à mieux comprendre les résultats du modèle mécaniste pour la dynamique du carbone des sols. Le coefficient de diffusion, qui traduit la bioturbation du sol, et le coefficient d'advection, lié à la diminution de la motilation, sont traités comme des constantes dans les modèles mécanistes développés pour la dynamique du carbone. Cependant, en réalité, ces coefficients ne sont pas constants et varient avec la profondeur. Un défi consisterait à transformer ces coefficients constants en fonctions, par exemple de type exponentielles décroissantes, de la profondeur. Les méthodes bayésienne de sélection explorées peuvent alors être utilisées pour définir les facteurs climatiques et environnementaux significatifs au sein de modèles statistico-mécanistes.

Journal de la Société Française de Statistique, Vol. 159 No. 2 128-155 http://www.sfds.asso.fr/journal © Société Française de Statistique et Société Mathématique de France (2018) ISSN: 2102-6238 Chapter 3

Jreich, Hatté, Parent and Balesdent

Journal de la Société Française de Statistique, Vol. 159 No. 2 128-155 http://www.sfds.asso.fr/journal © Société Française de Statistique et Société Mathématique de France (2018) ISSN: 2102-6238

to be checked for ergodicity although obtained in most cases

if π(θ ⋆ ) > π(θ (g) ) the candidate is accepted anyway, else the ratio π(θ ⋆ ) π(θ (g) )is compared to a uniform random draw

Acknowledgments

We are indebted to two anonymous reviewers for their attentive and detailed reviews and for their numerous constructive comments. The data set was kindly provided by Christine Hatté. This work is supported by the research contract ANR 14-CE-01-0004 and is part of the first author's PhD study financed by CEA.

JAGS

This chapter introduces a submitted paper, it is for this reason that it will not be available in this manuscript. difference between the two sub-models is that the PIP model considers, in addition, the Aridity Index (AI) as a significant predictor for the ω 1 latent linear model with a Posterior Inclusion Probability equal to the selection threshold 0.5. An increase in DIC is observed for the PIP model after adding the Aridity Index (AI) (DIC = 2078).

Thus, the final sub-model kept is the Median model.

Bayesian Sparse Group Selection

BSGS model specification and choice of hyperparameters

The Bayesian Sparse Group Selection (BSGS) model has the same structure as the BSGS proposed for the soil radiocarbon modeling with the exception that for the soil carbon content three latent variables are defined instead of four. This approach provides a probabilistic judgment about the inclusion of categorical predictors as well as levels. The Bayesian method of Kuo and Mallick (1998) was used to select the numerical predictors. For the prior specification, a Beta distribution with both hyperparameters equal to 2 was considered for the prior inclusion probability for categorical predictors as well as for levels. If the categorical predictor is significant, the regression effect is generated from a Zellner's g prior where the g value for each latent linear model is replaced by the posterior mean estimation of g obtained from the Full Bayesian model.

BSGS selection results

The best sub-model identified by the BSGS In the best sub-model, all the predictors have a Posterior Inclusion Probabilities (PIP) at least equal to 0.5. The best sub-predictors for each latent linear model are summarized in 5.16: Amount and type of clay generally observed in the soil types from the database. Types of soil are ranked according to the result of Bayesian Effect Fusion of soil type levels for deep soil radiocarbon activity (φ 1 ). Cluster number and color are the ones used in Figure 5.22. The number of profiles associated to each type of soil and each group of type of soils are provided in brackets. Column six refers to the median value of deep (higher than 100 cm deep) soil F 14 C from the database (line inside rectangle in Figure 5.22). The last column refers to 25% and 75% quantiles (q= quantile) (the upper and lower rectangle bounds in Figure 5.22) * The result of the clay amount and the type of clay corresponds to Arenosol soil type, only. If the proposal distribution does not fulfill the symmetry condition J(θ ′ → θ ) = J(θ → θ ′ ) , the ratio of acceptance has to be adapted. Many choices can be made for the proposal distribution , the simplest one being the Gaussian distribution centered on θ (g) with a variance to be carefully tuned so as to monitor the ratio of acceptance of the candidate. Tuning can be rather challenging, however, following the Goldilocks principle, according to [START_REF] Rosenthal | Optimising and adapting the metropolis algorithm[END_REF]. If the candidate is too rarely accepted (because the average jump to make a proposal -controlled by the variance of the proposal distribution-is too big), the chain remains stuck for long periods , but on the contrary, if the candidate is too frequently accepted (because the variance is too small) the chain wobbles within a small area of the parameter space). In both cases, the tuning of the jump function is bad, resulting in a slow convergence of the MCMC targeting the posterior distribution π(.).

Beyond Gibbs and Metropolis

There is a huge continuous development of stochastic algorithms devoted to improving the tools needed to perform Bayesian inference. Many paths for future developments can be anticipated.

• There have been many attempts to develop adaptive versions of MCMC [START_REF] Rosenthal | Optimising and adapting the metropolis algorithm[END_REF][START_REF] Atchadé | Limit theorems for some adaptive mcmc algorithms with subgeometric kernels[END_REF] but automatically tuning the algorithm remains a challenging task, because convergence properties of non homogeneous Markov Chains are difficult to establish . Relying on weighted independent draws such as the one obtained through the family of importance sampling techniques seems a more promising avenue of research [START_REF] Liu | A theoretical framework for sequential importance sampling with resampling[END_REF]. To design efficient importance sampling algorithms, the auxiliary proposal distribution should be chosen as close as possible to the posterior distribution but as the posterior distribution is unknown, choosing such a proposal is a tricky crucial task [START_REF] Gelman | Efficient metropolis jumping rules[END_REF][START_REF] Roberts | Weak convergence and optimal scaling of random walk metropolis algorithms[END_REF]. The adaptive multiple importance sampling algorithm of [START_REF] Cornuet | Adaptive multiple importance sampling[END_REF] is a good example of the efficiency of Population Monte Carlo methods [START_REF] Cappé | Population monte carlo[END_REF]. It consists in iteratively generating parameters under an adaptive proposal distribution and assigning weights to the parameter replicates. The AMIS algorithm sequentially tunes the coefficients of the proposal distribution (selected from a parametric family of distributions, generally the Gaussian one) and recomputes the weights of the cumulated posterior sample at the end of each iteration.

• Hamiltonian Monte Carlo algorithm -originally known as hybrid Monte Carlo [START_REF] Duane | Hybrid monte carlo[END_REF]-, is a most remarkable Markov chain Monte Carlo method for reducing the correlation between successive sampled states by using properties from Hamiltonian dynamics [START_REF] Neal | Mcmc using hamiltonian dynamics. Chapter 5 of Handbook of Markov Chain Monte Carlo[END_REF]. It allows larger moves between states at the cost of doubling the dimension of the state space and being able to efficiently compute the gradient of the logposterior density. This causes the algorithm to converge more quickly to the targeted posterior probability distribution. [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF] developed the No-U-Turn sampler, an adaptive form of Hamiltonian Monte Carlo sampling that is encoded in the package Stan. Stan, now available in a stabilized version, provides useful modern tools for Bayesian inference for continuous-variable models that are used in a wide range of application domains, e.g. in ecology [START_REF] Monnahan | Faster estimation of bayesian models in ecology using hamiltonian monte carlo[END_REF].

• Variational Bayes methods [START_REF] Beal | Variational algorithms for approximate Bayesian inference[END_REF] drop the idea of targeting the exact posterior distribution but search for an approximate (with respect to a Kullback-Leibler divergence) solution belonging to some convenient multivariate parametric family, by alternating quick optimization and expectation steps. A particularly successful variational approximation in pattern recognition and machine learning is the factorized form (Jordan et al., 1999a,b). The idea for such a factorized approximation stems from theoretical physics where it is called mean field theory [START_REF] Parisi | Statistical field theory[END_REF]. If one wishes nevertheless to get the exact posterior distribution, [START_REF] Donnet | Using deterministic approximations to accelerate smc for posterior sampling[END_REF] proposed a rather straightforward implementation of a bridge sampling scheme [START_REF] Gelman | Simulating normalizing constants: From importance sampling to bridge sampling to path sampling[END_REF] starting from a variational approximation of the posterior distribution and targeting the true one. Title : Vertical dynamics of soil carbon -Combined use of isotopic tracers and statistical meta-analysis Keywords : carbon ; geochemistry ; probabilistic modeling ; bayesian approach ; meta-analysis Abstract : Although it is the largest land reservoir interacting with the atmosphere, the response of the soil carbon reservoir to climate change and land use evolution remains uncertain. To better understand the soil carbon dynamics and assess the impact of. climate and environmental factors on residence time and soil organic carbon content, a non-linear hierarchical random effects model was proposed to model the variation in the responses of soil carbon content and soil radiocarbon as a function of depth. The selection of climatic and environmental factors was based on three Bayesian Selection techniques (Bayesian Group Lasso, Bayesian Sparse Group Selection and Bayesian Effect fusion) appropriate for categorical predictors (soil type and ecosystem type) and on the Stochastic Search Variable Selection for the numerical predictors (temperature, precipitation, etc.). The statistical modeling also enabled the effect of temperature increase and land-use conversion on soil carbon dynamics to be investigated

Conclusion Conclusion

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France