crestr An R package to perform probabilistic climate reconstructions using fossil proxies

Statistical climate reconstruction techniques are practical tools to study past climate variability from fossil proxy data. In particular, the methods based on probability density functions (PDFs) are powerful at producing robust results from various environments and proxies. However, accessing and curating the necessary calibration data, as well as the complexity of interpreting probabilistic results, often limit their use in palaeoclimatological studies. To address these problems, I present a new R package (crestr) to apply the CREST method (Climate REconstruction SofTware) on diverse palaeoecological 5 datasets. crestr includes a globally curated calibration dataset for six common climate proxies (i.e. plants, beetles, chironomids, rodents, foraminifera, and dinoflagellate cysts) that enables its use in most terrestrial and marine regions. The package can also be used with private data collections instead of, or in combination with, the provided dataset. It also includes a suite of graphical diagnostic tools to represent the data at each step of the reconstruction process and provide insights into the effect of the different modelling assumptions and external factors that underlie a reconstruction. With this R package, the CREST 10 method can now be used in a scriptable environment, thus simplifying its use and integration in existing workflows. I t is hoped that crestr will contribute to producing the much-needed quantified records from the many regions where climate reconstructions are currently lacking, despite the existence of suitable fossil records.

. (a) Conceptual illustration of the differences between a modelling approach based on the estimation of the most likely climate where all the probabilities are concentrated around the mode (light grey), and a modelling approach focused on the full spread of the data with the probabilities more spread along the climate gradient (dark grey). In both cases, the area under the curve sums to one. The two types of approaches are illustrated in (b) and (c) with two theoretical fossil samples (in blue and purple) representing two independent reconstructions of the same climatic parameter for the same time period. (b) The two reconstructions are derived from a method focused on only estimating the most likely climate value, resulting in 'apparently' incompatible reconstructions. (c) The same samples are reconstructed using an approach that estimates their complete uncertainty distributions. In this case, the response of the blue sample is broader, and the response of the purple sample becomes bimodal. When the full spread of these uncertainties is considered, the two reconstructions are not incompatible anymore, and a joint climate estimate (gold) can be derived from their overlapping sections (hashed polygons).
illustrating the successive stages of a CREST analysis and how to use the diagnostic tools to reproduce a recently published pollen-based temperature reconstruction .

CREST method & Calibration data
2.1 The probability density function (PDF) approach 60 As is standard with statistical climate reconstruction techniques, the core process of CREST can be decomposed into two major stages: 1) estimating the modern climatic responses of the proxies observed in the fossil sequence ( Fig. 2a-c) and 2) combining these responses to reconstruct past climates (Fig. 2d). In the following sections, the main elements of these two stages are presented along with all the parameters and/or modelling assumptions that can be modified in crestr. For an in-depth description of the method and its assumptions, the reader is, however, referred to . presented in the following sections, along with all the parameters and/or modelling assumptions that can be modified in 'crestr,.
For an in-depth description of the method and its assumptions, the reader is, however, referred to Chevalier et al. (2014). In CREST, PDFs are used to transform the information contained in the modern observations of biological climate proxies into probabilistic climate responses. A PDF thus represents a weighted ensemble of all the conditions where the proxy is observed today. PDFs can be fitted in one or two steps depending on the nature and taxonomic resolution of the studied proxy. Climate responses are first fitted at the species level (hereafter PDF sp (c, s) with c representing the studied climate variable and s a 75 species), and when necessary, these PDF sp (c, s) are then combined together to meet the taxonomic resolution of the fossil taxon (hereafter PDF tx (t, c) with t representing the observed taxon).
For all the species, their empirical mean climate and associated variance are calculated from their modern distributions following Eq. 1 and 2. These two parameters can be intuitively interpreted as the climate preference and tolerance of the species. For a robust estimation, it is recommended to exclude species with too few observations. Different studies have shown 80 that a threshold of a minimum of N s ≥ 20-25 distinct occurrences usually leads to robust estimates (e.g. Chevalier et al. (2014), Chevalier et al. (2021)). However, this number can vary between regions and climates. Each observation can also be weighted to account for the uneven distribution of modern climate (Kühl et al. (2002), Bray et al. (2006)). Extreme values are usually under-represented (see, for instance, the inset histogram on Fig. 2c), and this bias can push the estimation of the PDF sp (c, s) towards the mean climate observed across the study area (i.e. towards the centre of the "climate space"). It can also artificially 85 4 https://doi.org/10.5194/cp-2021-153 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.
shrink the posterior range of the reconstructions. In CREST, this weighting can be accounted for by first sorting the N climate values that compose the climate space into bins of equal sizes (e.g. 2 • C or 50 mm). Then, each climate value is given a weight defined as the inverse of the abundance of the bin it belongs to (Eq. 3). With this correction, abundant climate values are down-weighted, while the rarer ones are up-weighted, and the distribution of modern climate is more 'balanced'.
These two parameters, m s,c and s s,c , are used to define a regular, unimodal distribution for the PDF sp (s, c). In CREST, the shape of these distributions is pre-defined to be either normal (Eq. 4) or log-normal (Eq. 5). The decision should be made based on the climate variable to reconstruct. In general, log-normal distributions are recommended for variables that are undefined 95 for negative values, such as precipitation variables. represents the age or depth of the sample to reconstruct, and Fig. 2d). Due to their exponential nature, the PDF tx (t, c) rapidly converge to zero beyond the range of the modern observations of its composing species. As such, the multiplication of these PDF tx (t, c) ensures that only the climate values where all the taxa can coexist have a non-null (or, more precisely, a noninfinitesimal) probability. In addition, the presence or absence of each taxon is considered independent from the others. As such, it is possible to select a subset of climatically sensitive taxa to reconstruct each climate variable and maximise the 115 reconstruction signal (Chevalier and Chase, 2015), even if it is not always mandatory . These definitions of sensitive taxa are always specific to a specific region and variable and should not be generalised (Chevalier et al., in press).
Finally, the selected PDF tx (t, c) can be weighted with the parameter ω(t, z), which can take any positive values. For presence/absence observations, the weights of taxon t will be either one (the taxon is observed) or zero (the taxon is not observed).

120
For compositional data, the weights can be the observed percentages (i.e. values between 0 and 100) or rescaled percentages when the observed percentages are not directly proportional to the taxa composition across the catchment. In most cases, using raw percentages to weight taxa can give considerable weight to abundant, ubiquitous taxa that do not have a well-defined climate response and strongly limit the influence of rarer taxa with more explicit climate preferences.
An empirical normalisation is proposed in CREST to account for the varying production rates, distribution and preservation 125 processes impacting the relative proportions of fossil taxa observed in the sediments (Chevalier et al., 2014). This scaling entails calculating the average percentage of each taxon when it is present and dividing all its percentages by this taxon-specific scaling factor (Eq. 8, where O(t, z) represents the observed percentage of taxon t at age or depth z). With this normalisation, all the taxa vary on a standardised scale. The average presence is given a weight of 1, and values below and above this average presence threshold are assumed to represent lower and higher abundance in the environment. While the empirical nature of

CREST calibration dataset
A multiproxy calibration dataset to estimate PDFs from a global collection of geolocalised presence-only data (hereafter proxy distributions) was first presented in Chevalier (2019). These data were obtained from the Global Biodiversity Information  Tables 1 and   2). The QDGC spatial resolution is an empirical trade-off between numerous factors, including the resolution of the presence data, the quality of the data or the spatial representativity of the studied proxy (see discussions in Chevalier et al. (2014) and Chevalier (2019)). However, this trade-off may be suboptimal in some situations, and for that reason, crestr can also be used with the raw GBIF data and even alternative calibration datasets.

155
In its current version (V2), the gbif4crest calibration dataset contains about 25.3 million unique presence data for the six proxies. Unfortunately, the density of data available varies between proxies and regions ( Fig. 3 This solution is the recommended option, as users without any SQL knowledge can benefit from the package's interface to automatically query the database by providing study-specific parameters (e.g. the name of the taxa or geographical boundaries for the study area) to import all the necessary data in the correct format to the R environment (see Section 4.2). Second, advanced users can also directly query the database to extract and curate data from the DISTRIB or DISTRIB_QDGC tables 180 Table 2. List of marine variables available in the gbif4crest database. Each one can be selected in crestr using its associated code. silicate Silicate Concentration (µmol/L) Garcia et al. (2018b) using the dbRequest() function and subsequently associate these data with climate variables. Finally, the full gbif4crest calibration dataset can be downloaded as an SQLite3 portable database file from Chevalier (2020).
3 The crestr R package

Philosophy of the package
The crestr package has been designed for two independent but complementary modelling purposes. The probabilistic proxy-185 climate responses can be used to quantitatively reconstruct climate from their statistical combination, such as in Chevalier et al. (2021) or Chevalier and Chase (2015), or, alternatively, they can be used in a more qualitative way to determine the (relative) climate sensitivities of different taxa in a given area to characterise past ecological changes, such as in Chevalier et al. (in press) or Quick et al. (2021). For a simpler access to the different functionalities, default values are provided for all parameters to enable a rapid generation of preliminary results that users can then use as a starting point to fine-tune the model to their data. To guide users in this task and avoid the common 'black box' criticism faced by many complex statistical tools, a large array of publication-ready, graphical diagnostic tools was designed to automatically represent CREST data in a standardised way. These tools include plots of the calibration data, the estimated climate responses, the reconstructions and more. They can be used to look at the data from different perspectives to help interpret the results, and possibly to help identify potential issues or biases in the selected data and/or parameters. Such diagnostic tools are available for every stage of the process and, 195 as exemplified in section 4, they can be generated with a single line of R code.

The central element: the crestObj object
In crestr, all the CREST-related data are stored within a single S3 object of the class crestObj. Most functions of the package will take a crestObj as their main input and will return an updated version of that object. In practice, a crestObj is a nested list that contains five sub-lists, each one grouping a specific type of information, such as the calibration data, the  fitted climate responses, and the reconstructions (Fig. 5). Wrapper functions have been implemented to manipulate and/or modify the information contained in a crestObj so that, users are never expected to manually modify their crestObjeven if it is possible. The five sub-lists contain the following information: inputs: contains the raw data (e.g. the counts/percentages, the ages of the samples or the names of the fossil taxa).
parameters: contains the parameters provided at the different stages of the analysis (e.g. the data extraction or the 205 fitting and combination of the PDFs) modelling: contains all the data related to the estimation of the PDFs (e.g. the distribution data for calibration, the climate space or the PDFs themselves).
reconstructions: contains all the results (e.g. best estimates, synthetic error measurements as well as the full posterior distribution of the uncertainties).

210
misc: contains some additional metadata relative to the reconstruction (e.g. the site location or, most importantly, information relative to the proxy-species equivalency described in section 3.3.2).

Input data for crestr
Five different input data files are compatible with crestr. However, most applications will only require two file (the df and PSE files, see below) to be created. More specific applications may require up to four of these files. All the files can be prepared 215 outside the R environment and imported using standard R functions.

The fossil data (df)
The df data frame is only required if crestr is used for reconstructing climate and can be omitted if the objective is limited to modelling the climate response(s) of different taxa. df should be a data frame with the different samples entered as rows, with either the age, depth, or sample ID as the first column and the fossil data in the subsequent columns. df can contain raw 220 counts, percentages, presence/absence (1s and 0s) or even relative weights to be used in the reconstruction (see examples in Section 4.5).

The proxy-species equivalency (PSE) table
The PSE data frame is required to use the gbif4crest calibration dataset. It is used to associate fossil taxon names to the species distributions of various species stored in the TAXA and DISTRIB_QDGC tables (Fig. 4). When all the fossil taxa are identified 225 at the species level, the PSE table should be a simple data frame with one row per taxon. In most cases, however, fossil taxa are identified at a higher taxonomic level (sub-genus, genus, sub-family, family). These varying levels of identification should be encoded in the PSE file to link one or more (groups of) species to their common fossil taxon name.
In practice, a PSE file is composed of five columns (Table 3). The first one (Level) contains an integer that indicates the level of taxonomic resolution of the row (1 for Family, 2 for Genus, 3 for Species and 4 for taxa that should be excluded from the  reconstruction, e.g. 'Triletes spores' in the example case study below). The fifth column ProxyName contains the name of the taxon and columns two to four contain the taxonomic classification of that taxon as Family, Genus and Species, respectively.
For simplicity, a pre-formatted version of the PSE table with the names of all the taxa to study can be generated by crestr using the createPSE(list_of_taxa) function that generates a spreadsheet with the correct structure and with the Level and ProxyName columns automatically filled in.

235
When performing the species to proxy classification with the crest.get_modern_data() function (see Section 4.2), crestr first classifies the taxa with the lowest taxonomical resolution (i.e. when Level is equal to one), and then increases the resolution Level by Level. In the example in Table 3, different taxonomic resolution levels are provided for different plant species belonging to the highly diverse Asteraceae family (the daisy family). To distribute all the Asteraceae species observed across the study area to their appropriate category, all the species are first classified as 'Asteraceae undiff.'. In a 240 second time, the classification of some of these Asteraceae species is refined when reaching the better resolved sub-groups (Stoebe-type and Artemisia at Level 2). At the end of the process, the 'Asteraceae undiff.' group only contains Asteraceae species that grow in the study area but are not part of the genuses Stoebe, Elytropappus or Artemisia. The latter are categorised separately as Stoebe-type or Artemisia. Additional taxa can also be included to the PSE file to exclude species that are known to not be part of a group. For instance, this could have been used to simplify the climate response of the 'Asteraceae undiff.' 245 group by excluding more species from it, even if the pollen grains corresponding to these species have not been observed.
This categorisation process can be time-consuming, as all the taxa must be classified in a unique PSE table, and will often require many iterations to be optimised. The results of the different assignments are stored in the crestObj returned by the crest.get_modern_data() function and can be evaluated by checking rcnstrctn$misc$taxa_notes.

250
Users that prefer fitting proxy-climate relationships from their personal calibration data instead of the proposed gbif4crest dataset should prepare a distributions dataset following the specific structure presented in should contain the species names (or unique identifiers) and the corresponding proxy name, respectively. If more than one species correspond to one taxon, the PDFs will be fitted in two steps, as explained in section 2.1. The following two columns contain the coordinates of the species occurrence data. Finally, the last columns contain the climate values to be reconstructed.

255
An optional column called weight can be added to distributions in fifth position (i.e. between the coordinates and the climate variables) if one wants to weight the different observations. For example, the (relative) abundance of taxa observed from modern proxy samples can be used when fitting the PDFs to give more importance to the observations where that abundance is highest. This data frame is only necessary if the users use a personal calibration dataset instead of the gbif4crest dataset. It will enable the use of the climate space weighting (Section 2.1.1). Its structure is straightforward, with the first two columns containing longitudes and latitudes and the subsequent columns the climate variables to reconstruct. The spatial resolution and the ordering of the climate variables should be identical to the distributions table (Table 4).

The selectedTaxa data frame 265
The last data frame that may be used to inform the reconstruction is a data frame of ones and zeros called selectedTaxa.
This data frame has as many rows and columns as there are taxa and climate variables, respectively. Each entry should be either 1 or 0 and is used to indicate if the taxon should be used to reconstruct the climate variable (value = 1) or not (value = 0). If it is not provided, a default selectedTaxa data frame with all entries set to 1 is added to the crestObj.
Users can modify this information at any point using the includeTaxa() and excludeTaxa() built-in functions. The crest.get_modern_data() function also modifies this data frame by setting the value to -1 when the PSE classification failed for a taxon or when the amount of data in the study area is insufficient to fit a reliable PDF (see the parameters description in section 4.2).

Package dependencies
The crestr package is built in R (R Core Team, 2020) using the devtools package (Wickham et al., 2020). crestr 4 Step-by-step user guide for crestr

Example application: Pollen-based Mean Annual Temperature reconstructions from marine core MD96-2048
To illustrate the different ways of using crestr and its graphical diagnostic tools, I use pollen data that were recently analysed with the original CREST software to reconstruct mean annual temperature (MAT) from marine core MD96-2048. The core

Formatting the calibration data in a crestObj
As the gbif4crest dataset was used to fit the PDFs, using the function called crest.get_modern_data() is required to following points describe the different parameters of the crest.get_modern_data() and how they relate to the data extraction and modelling.
-The parameter taxaType is used to choose the type of proxy used and takes a value between 1 to 6 for plants (i.e. pollen and plant macroremains), beetles, chironomids, foraminifers, diatoms and rodents, respectively.
-The name(s) or code(s) of the climate variables to study should be provided here (see Tables 1 and 2 or use the function accClimateVariables() for a list of accepted values). Here, 'bio1' means mean annual temperature.

310
-Geographical parameters can be provided to create a regional subset of the gbif4crest dataset that is adapted to the studied data. These can include minimum and maximum longitude and latitude (xmn, xmx, ymn, ymx), continent, ocean or country names (see accCountryNames() and accBasinNames() for a list of accepted names), and also some ecological classifiers, such as realms, biomes or ecoregions (see accRealmNames() for a list of accepted names).
Only the occurrences that respect all the specified constraints will be returned.

315
-To estimate reliable PDFs, it is recommended to use at least 20 distinct occurrences for each species, but different values can be specified with the minGridCell parameter.
-Optional information about the site, such as a name and coordinates, can be provided and, where possible, this information will be represented on the different graphical diagnostic tools created by crestr (e.g. the location of the record can be added to the maps if the coordinates are provided). If the calibration data needed for this study were available as a distributions and climate_space data frames (Section 3.3), similar results would be obtained using the following command:

Estimating the climate responses (the PDFs)
The probabilistic proxy-climate relationships, i.e. the PDFs, are estimated from the presence data using the crest.calibrate() function. As described in section 2.1.1, the different parameters should be carefully considered to produce reliable PDFs. These include specifying:

365
the shape of the species PDFs, which should be either 'normal' or 'lognormal', depending on the variable to reconstruct.
the width of the climate bins ( the number of intervals required to divide the studied climate range and fit the PDFs npoints. This will ultimately 370 define the climate resolution of the reconstructions. crestr runs faster with a lower resolution but this can lead to a problem of aliasing of the reconstructions.
set geoWeighting to TRUE if the species PDFs of the different composing species should be weighted according to the square-root of the extent of their modern distribution.

Assessing the coherency of the climate space and climate responses
In every study involving estimating relationships between biological entities and environmental parameters, the first step is 385 always to ensure that the defined calibration dataset is as coherent as possible. This includes ensuring that 1) all the important taxa are present, and their distribution is not truncated, 2) the climate values to reconstruct are likely to be in the study area (the reconstructions are bounded by the lowest and highest values observed in the modern climate space) and 3) there is no large sampling or representativity bias (e.g. along country borders due to different sampling efforts). Ideally, the climate sampling should be as homogeneous as possible, even if the extreme climate values will always be under-represented compared to 390 the median ones. Equally importantly, the occurrence data across that space should be similarly distributed to ensure a proper sampling of the climate space. However, deviations from a theoretical one-to-one (or at least proportional) equivalence between climate and occurrence data abundance are not necessarily a bad characteristic. In our case study, the spatial variability represent true patterns in regional species diversity with the presence of several biodiversity hotspots across eastern and southern Africa (Myers et al., 2000). All these elements should be checked and accounted for while designing the final calibration dataset.

395
The 'climateSpace' graphical diagnostic tool was thus designed for a rapid assessment of all these characteristics (Fig.   6). This diagnostic figure is also very important to identify potential local or global correlations between different climate variables and assess the risks of confounding variables (Juggins (2013), Chevalier et al. (2020b)). Any change to any of the parameters related to the definition of the climate space (study area, climate variables to reconstruct) will require to re-run crest.get_modern_data() or crest.set_modern_data() with updated parameters and/or data. With the study area and climate space defined, the next step is to search for taxa that show specific relationships with climate. While all species eventually respond to all climate variables, within a given region they can be more sensitive to one over another. The low taxonomic resolution of some fossil proxies, such as pollen data, can also mask strong species-climate 410 relationships. Looking at each individual climate response(s) and assessing their significance within the boundary conditions of the study is thus critical. The 'taxaCharacteristics' diagnostic plot was designed for this task (Fig. 7). One summary plot can be generated for each taxon, where the geographical distributions and climate responses can be assessed and inter-compared.
As illustrated in Fig. 7, Ericaceae is preferentially observed in the colder environments of the study area, its higher percentages occur during glacial periods, and a coherent response of all its composing species can be observed despite a high 415 diversity (141 species with at least 20 unique occurrences across the study area). All these elements indicate that Ericaceae can be considered as an indicator of colder environmental conditions in eastern and southern Africa. Sensitivities to other variables can also be expected but are not considered in this study. Making such sensitivity inferences can be key to define a list of taxa that can be used to reconstruct temperature (Chevalier and Chase (2015), Chevalier et al. (2021)), but also to support qualitative interpretations of palaeoecological datasets (Chevalier et al. (in press), Quick et al. (2021)). Another useful diagnostic plot is the 'violinPDFs' (Fig. 8) that represents the PDFs of a selection of taxa on the same scale, which helps comparing the different responses. All the violins have the same area (all the probabilities sum to 1) and the taxa are ranked by increasing value of their temperature optima (i.e. the temperature corresponding to the peak of the PDF).

430
However, due to the possible multimodality of the PDFs and differences in tolerance ranges, this ranking does not mean that a taxon on the left always represents colder conditions than a taxon on the right. This is illustrated in many ways on Fig. 8 with, for example, Cassia-type that is estimated to experience warmer conditions than Coffea-type~61% of the time (based on 100,000 random draws from their respective PDFs) despite having a 'colder' climate optimum, or Diospyros that can tolerate much warmer conditions than most taxa with warmer optima. This type of representation can be particularly helpful to have 435 objective interpretations of ecological changes from pollen diagrams (Chevalier et al., in press, Quick et al. (2021)).
20 https://doi.org/10.5194/cp-2021-153 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. Figure 6. 'climateSpace' graphical diagnostic tool to evaluate the calibration dataset. The map on the right represents the density of unique species occurrences in each grid cell, here highlighting a certain bias towards South Africa. The lower abundance of plant data available from Angola and the Democratic Republic of Congo (data not shown) is the reason why these two countries were excluded from the study area. The map on the left represents the studied climate variable (MAT) across the study area. The double histogram in the middle represents the distribution of MAT (in grey, top), while the black histogram (bottom) represents how this distribution is sampled by the calibration data. Differences between the two histograms can be used to identify biases in the calibration dataset. Here, the small shift of the black histograms towards colder values (towards the left) is another way of seeing that more data are available from South Africa than other countries and reflects the regional patterns of biodiversity. If more variables had been selected in this study, additional rows would be added to the figure with a similar climate map and histograms and a scatterplot of the climate variables to highlight potential local or regional modern correlations.

Reconstructing climate
Along with the df data frame provided to the crest.get_modern_data() or crest.set_modern_data() functions, a set of reconstruction parameters have to be chosen to combine PDFs and estimate climate parameters. The selectedTaxa data frame that is stored in the crestObj (see Fig. 5) defines the taxa that will be used to reconstruct each climate variable 450 (by default, all taxa are included if sufficient data are available to fit a PDF). This selection can be modified by using the includeTaxa() and excludeTaxa() functions. Here for instance, both Aizoaceae and Chenopodiaceae/Amaranthaceae were excluded because they are not primarily sensitive to temperature in southern Africa. which the taxa will be always considered absent can also be provided to reduce the noise of the fossil dataset (e.g. pollen percentages lower than 1 or 2% are commonly excluded from reconstructions, Chevalier et al. (2020b)). If presenceThreshold is set to zero as is the case here, all the strictly positive pollen percentages are considered as true presences and used accordingly to reconstruct MAT. To weight the taxa as described in Eq. 7, four options are available in crestr: -The data can be converted to presence/absence with all the values above and below presenceThreshold being changed to ones and zeros, respectively. This option is recommended for data such as macrofossils for which relative 465 abundances cannot be reliably estimated.
-The data can be converted to percentages to weight the taxa according to their relative abundance. This option is recommended for data where reliable and direct proportions can be estimated.
-The data can be normalised following the method proposed by Chevalier et al. (2014) and described here by Eq. 8. This option is recommended for pollen data for instance.

470
-The data can be directly weighted by the values provided in df, which implies that users can define their own specific weighting strategy. Characterising the most important factors underlying a reconstruction is often complex. In this section, I present three different 480 graphical diagnostic tools that provide different perspectives on the reconstructed data and help opening the statistical black box in a visual way. First, the full probabilistic breadth of the reconstructions can be represented by using the standard R plot() function, which has been adapted to plot data stored in a crestObj (Fig. 9). While the probabilistic representation of the data should be preferred because it represents all the information available, a simpler version of this plot with the climate Figure 9. MAT probabilistic reconstructions for marine core MD96-2048. The yellow-green-blue colour gradient represents the uncertainties associated with each sample. This reconstruction is identical to the reconstruction presented in Chevalier et al. (2021) and the reader is referred to this publication for an in-depth validation and discussion of these results.
optima and more common uncertainty ranges represented as colour bands can be obtained by specifying simplify=TRUE 485 (not shown). The plot_combinedPDFs() function can then be used to identify the taxa that are driving the reconstruction by focusing on the sample level (Fig. 10). In this plot, the PDFs of all the taxa present in the sample and selected to reconstruct the climate Figure 10. 'combinedPDFs graphical diagnostic tool that shows the combination of the PDFs of all the taxa recorded in a sample (in colour).
The thickness of the lines is proportional to the weight of the taxa in the sample (absolute value indicated next to each taxon name). The black curve represents the posterior MAT reconstruction, from which a 'best' climate estimate can be estimated from the maximum of the curve and uncertainties derived by calculating the area under the curve.
variable are represented along with the reconstruction. This type of plot can help identify if there is a particular PDF that is at odds with the general assemblage, which can be indicative of a confounding factor. It is also useful to visualise the full spread 500 of the uncertainties and, by extension, highlight that reconstructions can be multimodal. The presence of multimodality can be the underlying cause of apparent noise in the reconstructions because minor changes in the taxa composition or percentages can force the system to oscillate between two maxima and thus 'appear' noisy, even if the background rate of change is minor.
This can also be seen from the reconstruction plot if the full uncertainties are represented (simplify=FALSE). The LOO analysis is a powerful tool to understand the taxa that primarily drive the reconstruction and the LOO results can be represented as a common stratigraphic diagram, where each row represents the effect of removing a taxon from the reconstruction (Fig. 11). For example, the expected net cold effect of including Ericaceae in the MAT reconstructions from MD96-2048 pollen record is immediately visible on Fig. 11. Similarly, the effect of removing either Ericaceae or Caryophyllaceae undiff.
from the sample dated at 64 kyr (represented on Fig. 10  This function takes the same parameters described for the 'step-by-step' functions with the same default values and may be preferred, for instance, when reconstructing several records at the same time.

Exporting the reconstructions
All the data stored in the crestObj can be easily exported from the R environment as spreadsheets and RData files using either export_pdfs() to save the climate responses of the studied taxa, or export() to save the reconstructions and many associated data in a publishable format. The latter also saves the crestObj as an RData file for easy reloading and/or sharing of the data.

Citing building elements
Finally, all the reconstructions derived from crestr are built on numerous, independent research efforts, including data compilations, modelling projects, statistical developments and/or software engineering. To support the long-term development 565 and/or support of all these building elements, it is crucial to always acknowledge all of them, even if their processing is mostly invisible to the users. As such, the list of references that must be cited for each use of crestr is automatically included in the summary tab of the spreadsheet generated by the export() function. In addition, the citation information can also be directly obtained from R using cite_crest(rcnstrctn). The function looks at the type of data that were used (e.g. the subset of the GBIF data were used, which climate variables) and returns a corresponding list of references. For example, a simple way of 570 crediting all the contributors for the MAT reconstruction from marine core MD96-2048 presented here could be the following: "To create this MAT reconstruction from the pollen record from marine core MD96-2048 (Dupont et al. (2019)), we employed the CREST method (Chevalier et al. (2014), Chevalier (2019)). The PDFs were estimated by combining the MAT field of Fick and Hijmans (2017) and the plant occurrence data of GBIF (GBIF (2020m), GBIF (2020k), GBIF (2020g)). The numerical analyses were realised with the crestr R package v1.0.0 (this reference)."

Perspectives and conclusion
CREST is a statistical method designed 1) to model probabilistic relationships between climate proxies and climate from modern observations and 2) to use these relationships to reconstruct past climate quantitatively. Initially, its use was enabled by a python point-and-click graphical user interface to favour accessibility to a diverse spectrum of palaeoclimatologists, palaeoceanographers and palaeoecologists, including those with limited coding experience (Chevalier et al., 2014). However, 580 such accessibility limited flexibility and how the method could be most effectively employed. Therefore, this article introduces an open-source R package to use CREST in a programmatic environment. The benefits of this transition are significant and include scriptability (i.e. the possibility of analysing many records automatically and sequentially), reproducibility (i.e. the capacity to reproduce an analysis) and better inter-operability (i.e. R packages are compatible with all computer systems). However, maintaining the highest level of accessibility remained at the core of the development process and is illustrated in The package is aimed at all researchers that are interested in using CREST to analyse palaeoecological datasets. In addition, its broad applicability will allow taking advantage of the recent growth of curated, open-access fossil datasets that created unprecedented opportunities to reconstruct climate from a wide range of proxies, particularly in regions where they are urgently 590 needed. This package will also contribute to the current transition from single-site to multi-site studies that is necessary to better understand past climate dynamics. However, it is essential to remember that running such techniques on several datasets should be done carefully, as many factors can impact the reconstruction process. Calibration datasets and reconstructions should always be assessed, possibly against independent evidence when they are available, even if there is, unfortunately, no single way to validate a reconstruction (but see Chevalier et al. (2020b) for some discussions on the generic principles).

595
Over time, the crestr package will be enriched with new functionalities to facilitate reconstruction validation. The online documentation will also be updated with diverse examples and tutorials based on real applications and assessments (https://mchevalier2.github.io/crestr/). Finally, bug reports, feedback, and suggestions for newer functionalities and graphical diagnostic tools are encouraged and can be transmitted to the author directly or through GitHub's bug report portal (https://github.com/mchevalier2/crestr/issues).