Inferring Past Environments from Biological Data: Progress – Problems - Potentialities

1. Inferring Past Environments from Biological Data: Progress – Problems - Potentialities John Birks University of Bergen, University College London, and University of Oxford INQUA July 2011

2. INTRODUCTION Early attempts at quantitative environmental reconstructions used presence of one or more ‘indicator species’ (e.g. Andersson, Samuelsson, Iversen, Grichuk, Coope). Major development in Quaternary science occurred in 1971 with publication of the classic paper by Imbrie & Kipp. Paper laid the foundation of calibration functions (transfer functions) as a tool for the quantitative reconstruction of past environments using the whole fossil assemblage, not just a few indicator species. Paradigm shift.

3. General Theory of Quantitative Reconstruction Presented by Imbrie & Kipp (1971) Y - biological responses (‘proxy data’) X - set of environmental variables that are assumed to be causally related to Y(e.g.SST) B - set of other environmental variables that together with X completely determine Y (e.g. trace nutrients) If Y is totally explicable as responses to variables represented by X and B, we have a deterministic model (no allowance for random factors, historical influences, etc.) Y = XB Imbrie & Kipp (1971)

4. Y = XB If B = 0 or is constant, we can model Y in terms of X and Re, a set of ecological response functions Y = X (Re) In palaeoecology we need to know Re. We cannot derive Re deductively from ecological studies and thus cannot build an explanatory model from our currently poor ecological knowledge. Instead we have to use direct empirical models based on observed patterns of Y in modern surface-samples in relation to X, to derive U, our empirical calibration function. Y = XU Imbrie & Kipp (1971)

5. Basic Approach to Quantitative Environmental Reconstruction – Calibration-in-Space 1, , m taxa Y f t samples Fossil data (e.g. pollen) ‘Proxy data’ Environmental variable (e.g. temperature) 1 variable X f Unknown To be estimated or reconstructed t samples To solve for Xf, need modern data about species and climate from n samples

6. 1, , m taxa Ym n samples Modern biology (e.g. pollen) Modern environment (e.g. temperature) 1 variable Xm n samples Model Ym in relation to Xm to derive modern calibration function Ûm Apply Ûm to Yf to estimate past environment Xf - Potential problems – discuss later

7. Calibration-in-space Xm Ym transfer function Ûm Yf Xf Based on an unpublished diagram by Steve Juggins

8. Basic Biological Assumptions Marine planktonic foraminifera - Imbrie & Kipp 1971 Foraminifera are a function of sea-surface temperature  Foraminifera can be used to reconstruct past sea-surface temperature Pollen is a function of regional vegetation Regional vegetation is a function of climate  Pollen is an indirect function of climate and can be used to reconstruct past regional climate Chironomids (aquatic non-biting midges) are a function of lake-water temperature Lake-water temperature is a function of climate  Chironomids are an indirect function of climate and can be used to reconstruct past climate Freshwater diatoms are a function of lake-water chemistry  Diatoms can be used to reconstruct past lake-water chemistry

9. Biological Proxy-Data Properties • Contain many taxa (200-300) • Contain many zero values (absences) • Commonly expressed as proportions or percentages - "closed" compositional data • Quantitative data are highly variable, invariably show a skewed distribution. Few common taxa, many rare taxa • Multicollinearity between variables • Can show spatial autocorrelation e.g. forams, dinocysts, pollen • Taxa generally have non-linear relationship with their environment, and the relationship is often a unimodal function of the environmental variables

10. Species Response Models A straight line displays the linear relation between the abundance value (y) of a species and an environmental variable (x). Modelled by linear regression. LINEAR A unimodal relation between the abundance value (y) of a species and an environmental variable (x). (u=optimum or mode; t=tolerance; c=maximum). Modelled by Gaussian logit regression (GLR) UNIMODAL

11. Environmental Data Properties • Generally few variables, often show a skewed distribution • Strong multicollinearity (e.g. July mean temperature, growing season duration, annual mean temperature) • Often difficult to obtain (few modern climate stations, corrections for altitude of sampling sites, etc.) • Strong spatial autocorrelation (tendency of values at sites close to each other to resemble one another more than randomly selected sites). Values at one site can be partially predicted from its values at neighbouring sites. Problem of nearly all data in real world. Recognised by Francis Galton in 1889. First methods to eliminate spurious correlation due to spatial position developed by ‘Student’ in 1914.

12. PROGRESS Since 1971, calibration functions widely used in palaeoceanography, terrestrial palaeoecology, and palaeolimnology Used with wide range of biological proxies • foraminifera, radiolaria, marine diatoms, coccolithophores • pollen, testate amoebae, mollusca, bryophytes, plant macrofossils • diatoms, chrysophytes, chironomids, ostracods, cladocerans Now many different numerical methods – at least 26 methods published

13. Reconstruction methods can be divided into three main types(Birks et al. 2010) • Indicator-species approach – one or many taxa considered as presence/absence • Similarity-based assemblage methods involving a quantitative comparison between past assemblages Yf and modern assemblages Ym (e.g. MAT, response surfaces) • Multivariate calibration methods involving a quantitative calibration functionÛm estimated from Xm and Ym, modern calibration or training data-set (e.g. weighted averaging regression and calibration) Concentrate on calibration-function approach

14. Approaches to Estimating Calibration Functions 1. Basic Numerical Models • Classical Approach (1) Y = f(X) + error Biology Environment (2) Estimate f by some mathematical procedure and 'invert' estimated (f) to find unknown past environment Xffrom fossil data Yf Xf f-1(Yf)

15. Inverse Approach In practice, for various mathematical reasons, do an inverse regression or calibration (3) X = g(Y) + error (4) Xf = g(Yf) Obtain 'plug-in' estimate of past environment Xf from fossil data Yf f or g are calibration functions Easier to compute and nearly always performs as well as classical approach

16. 2. Assumed Species Response Model • Linear or unimodal • No response model assumed (linear or non-linear) • 3. Dimensionality of Model • Full (all species considered) • Reduced (selected components of species used) • 4. Estimation Procedure for Model • Global (estimate parametric functions, extrapolation possible) • Local (estimate non-parametric functions, extrapolation not possible)

17. Commonly Used Methods I = inverse; C = classical L = linear; U = unimodal; NA = not assumed; R = reduced dimensionality; F = full dimensionality; G = global parametric estimation; Ln = local non-parametric estimation CF = calibration-function based; S = similarity-based

18. Good reasons for preferring methods with assumed biological response model, full dimensionality, and global parametric estimation(ter Braak (1995), ter Braak et al. (1993), etc.) • Can test statistically if taxon A has a statistically significant relation to particular environmental variables • Can develop ‘artificial’ simulated data with realistic assumptions for numerical ‘experiments’ • Such methods have clear and testable assumptions – less of a ‘black box’ than e.g. artificial neural networks • Can develop model evaluation or diagnostic procedures analogous to regression diagnostics in statistical modelling • Having a statistical basis, can adopt well-established principles of statistical model selection and testing. Minimises ‘ad hoc’ aspects of MAT “To make sense of an observation, everyone needs a model … whether he or she knows it or not” Marc Kéry (2010)

19. Basic Requirements in Quantitative Palaeoenvironmental Reconstructions • Need biological system with abundant fossils that is responsive and sensitive to environmental variables of interest.  • Need a large, high-quality training set of modern samples. Should be representative of the likely range of variables, be of consistent taxonomy and nomenclature, be of highest possible taxonomic detail, be of comparable quality (methodology, count size, etc.), and be from the same sedimentary environment.  • Need fossil set of comparable taxonomy, nomenclature, quality, and sedimentary environment. 

20. 4. Need robust statistical methods for regression and calibration that can adequately model taxa and their environment with the lowest possible error of prediction and the lowest bias possible and sound methods for model selection.  5. Need means of establishing if reconstruction is statistically significant. 6. Need statistical estimation of standard errors of prediction for each reconstructed value.  7. Need statistical and ecological evaluation and validation of the reconstruction and of each reconstructed value.

21. PC1 PC2 Ym Xm PC3 Major Methods Used Principal components regression (PCR) = Imbrie & Kipp (1971) approach Multiple linear regression or quadratic regression of X on PC1, PC2, PC3, etc, to derive Ûm. Express Yf as principal components and apply Ûm to estimate Xf Principal components maximise variance withinYonly Selection of PCA components done visually until recently. Now cross-validation is used to select model with fewest components, lowest root mean square error of prediction (RMSEP), & lowest maximum bias. ‘Minimal adequate model’ in statistical modelling Inverse, linear, reduced dimensionality, global estimation. Linear response model is assumed, although non-linear responses are possible.

22. Gaussian logit regression (GLR) and maximum likelihood (ML) calibration ter Braak & van Dam (1989) b0, b1, b2 ML calibration Ym + Xm Yf Xf b0, b1, b2 environmental reconstruction modern data fossil data b0, b1, b2 taxon GLR regression coefficients for all taxa Ûm Classical, unimodal, full dimensionality, global estimation. Robust to spatial autocorrelation. Can be computationally difficult. ML finds the most likely value of Xf that maximises the likelihood function given Yf and Ûm

23. Two-way weighted averaging regression and calibration (WA) ter Braak & van Dam (1989); Birks et al. (1990) U1 WA regression WA calibration Ym + Xm Yf Xf U2 modern data environmental reconstruction fossil data Ut taxa WA optima ‘calibration function’ Ûm Inverse, unimodal, full dimensionality, global parametric estimation. Robust to spatial autocorrelation. First used in Quaternary science by Lynts and Judd (1971) Science 171: 1143-1144

24. Ecologically plausible – based on unimodal species response model. • Mathematically simple but has a rigorous mathematical theory. Properties fairly well known now. • Empirically powerful: • does not assume linear responses • not hindered by too many taxa, in fact helped by many taxa! Full dimensionality • relatively insensitive to outliers • Tests with simulated and real data – at its best with noisy, taxon-rich compositional percentage data with many zero values over long environmental gradients. • Because of its computational simplicity, can derive error estimates for predicted inferred values by bootstrapping. • Does well in ‘non-analogue’ situations as it is not based on the assemblage as a whole but on INDIVIDUAL taxa optima and/or tolerances. Robust to spatial autocorrelation. Globalparametric estimation. • Ignores absences of taxa.

25. WA WA GLR GLR pH • Weaknesses • Sensitive to distribution of environmental variable in training set, leading to ‘edge effects’ where responses are truncated. J. Oksanen (2002) 2. Disregards residual correlations in biological data. Can extend WA to WA-partial least squares to include residual correlations in biological data in an attempt to improve estimates of taxon optima

26. Weighted averaging partial least squares regression and calibration (WA-PLS)ter Braak & Juggins (1993) and ter Braak et al. (1993) WA-PLS regression PLS1 WA-PLS calibration Ym Xm βm Yf Xf PLS2 coefficients (Ûm) PLS3 Components selected to maximise covariance between taxon weighted averages and environmental variable X Selection of number of PLS components to include based on cross-validation. Model selected should have fewest components possible and low RMSEP and maximum bias – minimal adequate model. Inverse, unimodal, reduced dimensionality, global parametric estimation. Can be sensitive to spatial autocorrelation.

27. Comparison of different methods Imbrie & Kipp (1971) data Model performance statistic is root mean squared error of prediction (RMSEP) based on leave-one-out cross-validation Linear Unimodal Shows importance of using a unimodal-based method(ter Braak et al. (1993))

28. Other Areas of Progress Besides the development of new methods for deriving calibration functions and of modern calibration data-sets, there have been major developments in model evaluation and selection and in reconstruction assessment, namely statistics of calibration functions and in understanding the strengths and weaknesses of different methods and in their underlying theory

29. 1. Model evaluation and selection Tendency to use several different methods and to select so-called ‘best’ method. Resulted in a shift from an obsession with the model with lowest RMSEP or, even worse, the highest r2. More concern with model performance statistics including estimates of bias and number of components fitted (e.g. in WA-PLS). Model performance usually based on some form of internal cross-validation (leave-one-out, n-fold cross-validation, or bootstrapping) or external cross-validation with independent test-set.

30. van der Voet (1994) randomisation test of models helps find ‘minimal adequate model’ (MAM). Model with good performance statistics and fewest number of fitted parameters. May be more than one MAM. More work needed on model selection using criteria like Akaike Information Criterion (AIC) where unnecessary parameters are penalised. Active research area in ecology and evolutionary biology today. Of course, performance of modern model is being assessed with other modern data, not with fossil data! Major problem.

31. 2. Effects of spatial autocorrelation Estimating model performance in terms of RMSEP, r2, maximum bias, etc, assumes that the test-set is statistically independent of the training-set. Cross-validation in presence of spatial autocorrelation violates this assumption as test samples are not spatially and statistically independent. Spatial autocorrelation property of almost all environmental data and much ecological and biological data. Telford & Birks 2005 Quat. Sci. Rev. 24: 2173-2179 Telford 2006 Quat. Sci. Rev. 25: 1375-1382 Telford & Birks 2009 Quat. Sci. Rev. 28: 1309-1316

32. Results show the apparent performance of some models is enhanced as a result of spatial autocorrelation in oceans and on land Problems in finding spatially independent test-sets. Telford & Birks (2009) have developed methods for cross-validating a calibration function in presence of spatial autocorrelation, h-block cross-validation

33. 3. Partitioning Root Mean Square Error of Prediction Model uncertainty commonly expressed as RMSEP

34. 4. Testing the statistical significance of a quantitative palaeoenvironmental reconstruction All calibration-function programs will produce output or ‘reconstruction’ Does the resulting reconstruction explain more of the variance in the fossil data than most (say 95%) reconstructions derived from calibration functions trained on random environmental data? If it does, then it is statistically significant. Telford & Birks 2011 Quat. Sci. Rev. 30: 1272-1278

35. 5. Evaluation of individual reconstructed estimates Assuming overall reconstruction is statistically significant, some individual estimates may be less reliable than others (poor preservation, unusual composition or peak, etc). Need to evaluate individual reconstructed values. • Goodness-of-fit measures for each individual fossil sample, as in regression modelling (Birks et al. 1990) • Analogue statistics (Birks et al. 1990; Simpson 2007) • Proportions of taxa in fossil assemblage absent or rare in modern training data with no or poorly estimated taxon parameters (Birks 1998) • Sample-specific errors for reconstructed values estimated by bootstrapping or Monte Carlo simulation (Birks et al. 1990)

36. What to do with sample-specific errors? Has a statistically significant (p=0.009) reconstruction but there is also a continuous overlap in RMSEP. Problems of temporal autocorrelation in assessing RMSEP for samples. Birks & Peglar (unpub.)

37. 6. Highlighting ‘signal’ from ‘noise’ in reconstructions Use of LOESS smoother a great help Sample-specific errors or LOESS smoother Seppä & Birks (2002) Brooks & Birks (2001)

38. 7. Ecological validation Compare reconstructed values with historical data. Rarely possible as few historical data exist. Renberg & Hultberg (1992) But when done, sometimes the model that gives the closest correspondence is not the model with lowest RMSEP or maximum bias! Conflict between model performance and selection based on cross-validation of modern data and validation results using independent historical test-sets

39. 8. Palaeoecological validation by multi-proxy data Birks & Ammann (2000) Similar trends, different absolute values. Not surprising, given different biology of different groups of organisms

40. PROBLEMS • Violation of assumptions • Multiple –variable reconstructions • Spatial autocorrelation and non-independent test-set • Confounding effects of correlated environmental variables • Assumption of uniformitarianism • Different results from different proxies

41. The biggest set of problems is that the calibration-function approach, like any other quantitative procedure, makes assumptions, as originally stated by Imbrie & Kipp (1971), Imbrie & Webb (1981), and Birks et al. (1990). • These assumptions are being increasingly violated, especially in the last 5-10 years. • What are these assumptions?

42. Assumptions in quantitative palaeoenvironmental reconstructions • Taxa in training set (Ym) are systematically related to the physical environment (Xm) in which they live • Environmental variable (Xf , e.g. summer temperature) to be reconstructed is, or is linearily related to, an ecologically important variable in the system • Taxa in the training set (Ym) are the same as in the fossil data (Yf) and their ecological responses (Ûm) have not changed significantly over the timespan represented by the fossil assemblage • Mathematical methods used in regression and calibration adequately model the biological responses (Um) to the environmental variable (Xm) • Otherenvironmental variables than, say, summer temperature have negligible influence, or their joint distribution with summer temperature in the fossil set is the same as in the training set • In model evaluation by cross-validation, the test-data are independent of the training data • Imbrie & Kipp (1971), Imbrie & Webb (1981), Birks et al. (1990), Telford & Birks (2005)

43. Multiple-variable reconstructions • Increasing tendency to reconstruct 2 or 3, even 7-8, environmental variables that on the basis of current ecological knowledge of, e.g., vegetation, chironomids, or diatoms, cannot all be ‘ecologically important’ (assumption 2) • e.g. mean January, mean July, mean annual temperature, growing degree days above 0C and above 5C, annual precipitation, and evaporation : potential evaporation. • Ecological data are not usually influenced by 8 independent ‘ecologically important’ variables. Usually only 1-3 significant ordination axes. • All variables may be statistically significant in a RDA or CCA when considered individually(‘marginal’ effects) but almost certainly not significant when considered together (‘conditional’ effects, high multicollinearity, variance inflation factors). Many reconstructions of, for example, ‘distance to littoral vegetation’ suspect.

44. Problems of spatial autocorrelation and lack of independence in cross-validation test-data (assumption 6) • Spatial autocorrelation results in highly optimistic model performance statistics when data are spatially autocorrelated, especially when MAT (or ANN) are used (local non-parametric estimation procedures) • h-block cross-validation – allows for spatial autocorrelation, unlike leave-one-out c-v • Increase of RMSEP in h-block c-v for all data except diatoms & pH. Greatest increase in RMSEP always in MAT and in training-sets with the highest spatial autocorrelation in the X variable (forams & salinity; forams & SST; pollen & July sunshine)

45. 4. Confounding effects of correlated environmental variables Present in all studies, starting with Imbrie & Kipp (1971) with reconstructions of summer and winter sea-surface temperature and salinity. Covarying environmental variables e.g. temperature and lake trophic status (e.g. total N or P) or temperature and lake depth and chironomids. Is the fossil chironomid signal temperature or trophic status? Broderson & Anderson (2002)

46. In almost all ecological systems, assemblages are a complex function of multiple climatic, edaphic, land-use, biotic, and historical factors. First part of assumption 5 (environmental variables other than the variable being reconstructed have negligible influence) is therefore almost never met. Need very careful design of modern training-set and rigorous statistical analysis to establish what can reliably and significantly be reconstructed. Second part of assumption 5 (the joint distribution of additional variables with the variable of interest does not change with time) is also violated in many cases.

47. Climate model and glaciological results suggest that the joint distribution between summer temperature and winter accumulation has not been the same in the past 11,000 years. Good evidence to suggest that lake-water pH has decreased naturally (soil deterioration) whilst summer temperature rose and then fell in the last 11,000 years. pH-climate relationship changed with time. In Norway today, lake-water pH is negatively correlated with summer temperature because lakes of pH 6-7.5 are on basic rock and this happens in Norway to occur mainly at high altitudes and hence at low temperatures. In the past after deglaciation, almost all lakes had a higher pH than today, so the pH-temperature relationship in the past was different than today.

48. 5. Assumption 3 “Taxa in the training-set are the same as in the fossil data and their ecological responses have not changed significantly over the timespan represented by the fossil assemblage” Assumption not unique to calibration functions. Basic assumption of all Quaternary palaeoecology, namely uniformitarianism. Considerable interest in niche-conservatism amongst biogeographers and conservation and evolutionary biologists. Increasing evidence for conservatism of ecological niche characteristics in the timespan of last 20,000 years. Problems of ‘cryptic’ species and of taxa like Saxifraga oppositifolia-type in environmental reconstructions currently unresolved.

49. 6. Use of different proxies can give different reconstructions Mean July temp, Bjørnfjell p = 0.001 p = 0.183 ns Validate using another proxy – e.g. macrofossils of tree birch Validate using second proxy – e.g. chironomids Importance of independent validation and establishing what is statistically significant

50. POTENTIALITIES Quantitative palaeoenvironmental reconstructions in the context of Quaternary palaeoecology are not really an end in themselves (in contrast to Quaternary palaeoclimatology) but they are a meansto an end. Use the reconstructions based on one proxy (e.g. chironomids) to provide an environmental history against which observed biological changes in another, independent proxy (e.g. pollen) can be viewed and interpreted as biological responses to environmental change.