Variance inflation factors and ordination model selection

Variance inflation factors (VIF) give a measure of the extent of multicollinearity in the predictors of a regression. If the VIF of a predictor is high, it indicates that that predictor is highly correlated with other predictors, it contains little or no unique information, and there is redundancy in the set of predictors.

We don’t really want to have redundant predictors in a constrained ordination: they make the analysis more difficult to interpret, and the more predictors we have, the less constrained the ordination is, the more it resembles an unconstrained ordination.

Model selection in ordinations is tricky. Techniques such as forward selection can guide the choice of predictors to include, but also bias the results (see Juggins 2013).

Recently I have seen a paper and reviewed a manuscript using VIF to determine which predictor to drop from a constrained ordination in a backwards selection.

The procedure they were using was

  1. Generate a constrained ordination with all available predictors.
  2. Calculate the VIF of each variable.
  3. If any variable has a VIF over a threshold (typically 20), drop the variable with the highest VIF
  4. Repeat until all remaining variables have a VIF below the threshold.

I want to show that this is a procedure that can go badly wrong using an example based on the SWAP diatom-pH calibration set.

First I load the data from the rioja package and make two fake pH variables that are highly correlated with the observed pH by adding Gaussian noise to the observed pH.


env<-data.frame(pH=SWAP$pH, fakepH1=rnorm(length(SWAP$pH), SWAP$pH, .2), fakepH2=rnorm(length(SWAP$pH), SWAP$pH, .2))
pairs(env)#highly correlated.

Now I fit a constrained correspondence analysis using cca, and calculate the VIF.

mod<-cca(sqrt(SWAP$spec)~., env)
#      pH  fakepH1  fakepH2
#26.82960 12.63575 14.49851
CCA of the SWAP calibration set with the predictor pH between two fake predictors correlated with pH

CCA of the SWAP calibration set with the predictor pH between two fake predictors correlated with pH

The VIF for pH is greater than the threshold 20, and so under the above procedure is dropped. If we fit a new model without pH we find the VIF of the remaining predictors has fallen and all are now below the threshold. An anova shows that both predictors are statistically significant

mod2<-cca(sqrt(SWAP$spec)~.-pH, env)
# fakepH1  fakepH2
#6.822985 6.822985

anova(mod2, by="terms")
#Permutation test for cca under reduced model
#Terms added sequentially (first to last)
#Model: cca(formula = sqrt(SWAP$spec) ~ (pH + fakepH1 + fakepH2) - pH, data = env)
#          Df   Chisq       F N.Perm Pr(>F)
#fakepH1    1  0.3406 11.0050     99   0.01 **
#fakepH2    1  0.0530  1.7122     99   0.01 **
#Residual 164  5.0754
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

What’s happening is that pH is highly correlated with both fake predictors, whereas the fake predictors are less strongly correlated with each other so have a lower VIF.

Any procedure with a high risk of throwing out the true predictor and instead accepting two fake predictors is not a useful method.

Is my toy example realistic? I would argue it is. One of the studies I saw used the temperature of all the months of a year as predictors for a calibration set. Imagine that June is the true predictor, but June temperature will be highly correlated with May and July temperatures. The procedure risks returning three or four months spaced evenly around the calendar rather than just June.

VIF is a very useful indicator that there is multicollinearity in a data set. It is a poor indicator of which predictor should be dropped from a model.

Posted in EDA, R | Tagged | 3 Comments

On the niche of Thalassiosira faurii, perils in palaeoecology

Palaeoenvironmental conditions can be reconstructed from microfossils preserved in sediments using the relationship between species and the environment. Usually the species’ relationships with the environment — their niches — are insufficiently constrained by experimental data, so we are forced to rely upon observations of species abundances and environmental conditions in a modern calibration set and assume that we can make useful inferences of the niche from these data.

Sometimes this assumption fails. This post tells of one such case.

I did my PhD, and some subsequent work, trying to infer climate from diatom assemblages in sediment cores from closed-basin lakes in Ethiopia, and later Tanzania. Closed-basin lakes are lakes that don’t have a river flowing out of them (or groundwater leakage), and should be sensitive to climate change:  deep and relatively fresh when the climate is humid; and shallow and saline during arid periods. Diatoms, algae with beautiful siliceous cell walls or frustules that preserve well in most lake sediments, are sensitive to salinity and water depth and so should be possible to use to constrain past changes in hydrology. Françoise Gasse and coworkers (Gasse et al 1995) developed a transfer function to reconstruct conductivity (a variable related to salinity) from diatom assemblages. I used this transfer function in my thesis and subsequent work.

The first paper I wrote was on a 6500-year long diatom stratigraphy from Lake Awassa, a large but fairly shallow caldera lake in the Ethiopian rift valley. Even though the lake has no out-flowing rivers, it is fairly fresh, suggesting that there is groundwater leakage, probably to the north. I was expecting to find evidence that Lake Awassa was deeper and fresher during the early Holocene, which was more humid in most of North Africa, with lakes scattered across the Sahara. I was hoping to be able to make a continuous climate reconstruction across the end of the termination of the early Holocene humid period – was it an abrupt or gradual change? Instead, the diatom-conductivity model suggested that the lake was more saline during two phases of the mid-Holocene, before the expected termination of the humid period. Either the timing and nature of the termination inferred from several sites across North Africa were incorrect, or there was something funny going on in Lake Awassa. I decided that the latter was more plausible.

Both of the apparently saline phases in Lake Awassa occurred immediately after very fine-grained tephra layers, so I investigate the potential for a volcanic origin for the saline phases. I found sufficient evidence from other sites to make what I thought was a plausible case for pulses of saline groundwater somehow related to the volcanism making the lake more saline.

One of the main salinity indicators in my Lake Awassa core was Thalassiosira faurii. Most Thalassiosira are marine species, but some live in, often brackish, lakes.

Thalassiosira faurii (Roubeix et al 2014)

Thalassiosira faurii (Roubeix et al 2014)

Gasse et al (1995) estimated the conductivity optimum of T. faurii to be ~9000 μScm-1, calculated by weighted-averaging. With the European Diatom Database (EDDI) we can look at the data behind this estimate, and plot them (you might need to enable java for unsigned applets – sorry I wrote/borrowed the java code for these plots long before security holes in java were known).

T. faurii relative abundance against salinity in the African diatom calibration set

T. faurii relative abundance against salinity in the African diatom calibration set

The marked sample and the other two samples with a high relative abundance of T. faurii are all from a salt swamp in Niger. Without these samples, the optimum would be somewhat lower. This dependence on one site, where T. faurii does not reach the maximum relative abundance (60%) found in Late Awassa, is not ideal.

Roubeix et al (2014) investigated the salinity niche of T. faurii collected from Lake Langano, another rift-valley lake 60 km to the north of Awassa, by measuring its growth rate when grown at different salinities under controlled conditions in a laboratory.

Growth rates of Thalassiosira faurii (black diamonds) and Anomoeoneis sphaerophora (open diamonds) versus conductivity of the culture medium. The dashed line shows the variations of pH between cultures. Vertical bars represent the standard deviations of the triplicates. A star means that the cells were still alive after 15 days, whereas a cross means that they were all dead. The two triangles indicate the conductivity and pH of lake water from which the diatoms were isolated.

Growth rates of Thalassiosira faurii (black diamonds) and Anomoeoneis sphaerophora (open diamonds) versus conductivity of the culture medium. The dashed line shows the variations of pH between cultures. Vertical bars represent the standard deviations of the replicates. A star means that the cells were still alive after 15 days, whereas a cross means that they were all dead. The two triangles indicate the conductivity and pH of lake water from which the diatoms were isolated.

At least the clone grown by Roubeix et al (2014) did not have an conductivity optimum of ~9000 μScm-1; it died at conductivities this high. The optimal conductivity was estimated at just  400 μScm-1. This is the optimum in the absence of any competition, the realised optimum in lakes many be somewhat higher, but is unlikely to be much above ~1000 μScm-1, as by ~2000 μScm-1 there is no growth. Roubeix et al use this new information to reinterpret the conditions in Lake Abiyata, a saline lake adjacent to Langano, during the Younger Dryas. Previous work by Chalié and Gasse (2002) had inferred saline conditions during this interval.

Roubeix et al also question my reconstruction from Lake Awassa. If T. faurii has a much lower conductivity optimum, then the conductivity spikes I reconstruct in the mid-Holocene may be an artefact. While I recognise that my proposed mechanism of pulses of saline hydrothermal water is speculative, I’m not entirely convinced that I was wrong as the salinity spike is not based on the presence of T. faurii alone. Instead, I found T. faurii in association with the saline indicators Navicula elkab (optimum ~14000 μScm-1) and T. rudolfii (optimum ~11400 μScm-1). Both these conductivity optima estimates seem to be well supported by data from a variety of sites.

The discrepancy between the old observational and new experimental optimum of T. faurii may reflect taxonomic problems. Perhaps T. faurii is not a single species but a suite of morphologically similar cryptic species with different salinity niches, with the diatoms from the salt swamp in Niger representing one end of the spectrum and the clone from Langano representing the other. Unless morphological characteristics that distinguish salt-tolerant from non-tolerant variants of T. faurii, it will not be possible to guarantee which variant is present in the palaeoecological records and any reconstruction will be rather uncertain.

Cryptic species are known in several of the taxonomic groups used for reconstructing climate, so there is a risk of cryptic species with different ecological preferences may be a widespread problem in palaeoecology.

Part of the problem relates to the lack of external data. If I examine a pollen calibration set and find the optimum temperature for lime (Tilia) to be 5°C, there is enough literature on the ecology of Tilia to recognise that there is a problem. With diatoms, dinoflagellate cysts and many other groups, there is little or no literature on the ecology of the different species, so we cannot easily recognise when the optima we calculate are spurious. Yet another reason why care is required when using transfer functions to reconstruct palaeoenvironmental conditions.

Gasse, F., Juggins, S. & Ben Khelifa, K. 1995. Diatom-based transfer functions for inferring past hydrochemical characteristics of African lakes. Palaeogeogr. Palaeoclim. Palaeoecol. 117: 31–54.

Roubeix, V., Chalié, F. & Gasse, F. (2014) The diatom Thalassiosira faurii (Gasse) Hasle in the Ziway–Shala lakes (Ethiopia) and implications for paleoclimatic reconstructions: Case study of the Glacial–Holocene transition in East Africa. Palaeogeogr. Palaeoclim. Palaeoecol. 402: 104–112.

Posted in climate, Palaeohydrology, Peer reviewed literature, transfer function | Tagged , | Leave a comment

Smilodon and the climate skeptic

I had been meaning to write about La Brea Tar Pits, having spent a morning with the mammoths last month, paying homage to the bones of the Pleistocene megafauna, the relics of extinction.


Mammoth columbi  (my photo)

Researchers from the museum have published two papers on the evolution of two of the predators found at La Brea, Smilodon fatalis the sabre-toothed-cat, and Canis dirus the dire wolf, in response to climate change in the late Pleistocene. The papers are based on the cranial morphology of the two species, which, especially the dire wolf, are present in large numbers in the tar pits. For Smilodon, they find a large morphotype during warmer periods of the late Pleistocene and a small morphotype in the colder periods, which they suggest may relate to the choice of prey species.


Smilodon fatalis (my photo)

Anthony Watts has written an incisive and intelligent commentary about the press release announcing these two papers in a post “The La Brea Tars Pits gets themselves in a sticky wicket over climate change and adaptation“.

One of the most shrill arguments from alarmists is the idea that climate change will wipe out species because they can’t adapt. The claims run from polar bears to tortoises, to plants and coral. Yes, if we listen to these arguments, Nature so poorly equipped it’s creatures that they can’t adapt to a slightly warmer future.

Except when the last ice age ended, and it got warmer, and the saber-toothed cats and wolves got bigger because the prey got bigger…instead of disappearing due to “climate change”.

OK, so that wasn’t particularly incisive or intelligent. Watts  is comparing the evolutionary response of large, free-roaming populations of carnivores over millennia, to what might be expected for the remnant populations of wildlife, constrained by agriculture and urbanisation to small and fragmented areas, over a century.

Posted in Fake climate sceptics, Peer reviewed literature, Silliness, WUWT | Tagged , , | 3 Comments

INTIMATE training school 2014

Don’t try searching for “INTIMATE dating” on Google; you won’t find what you are looking for.

This is what you want – a week long training school in palaeoecological and palaeoclimate techniques, organised by the INTIMATE  network (INTegrating Ice core Marine and TErrestrial records) in Romania. The training school will be held by the volcanic crater Lake St. Anne between May 30 and June 5, 2014. I’ll be giving a session on developing transfer functions and interpreting palaeoecological data using R.

The course is open to end-MSc students, PhD fellows, and early-stage postdocs. The application deadline is 27th April. 

Here is a video of last year’s INTIMATE training school.

Posted in climate, transfer function | Tagged , | Leave a comment

Clumped Isotopes and the Pacific Warm Pool at the Last Glacial Maximum

How warm were the tropical oceans at the Last Glacial Maximum (LGM), 21 thousand years ago? A good estimate would be valuable for several reasons. First, it would help us to understand LGM climates and the ecological response to this climate. Second, it would provide a constraint by which to judge climate models when they are used to simulate LGM climates – if the climate model output resembles the past conditions, it increases our confidence that they give good projections of future climate with increased greenhouse gases. Third, the change in temperature between the LGM and the pre-industrial can be used to estimate climate sensitivity as the forcings (greenhouse gases, albedo and orbital)  are known.

The first global estimate of SST changes since the LGM, which came from the CLIMAP project in the 1970s and 1980s, mainly based on planktonic foraminifera, showed limited cooling in the tropics and warming in the subtropical Pacific.

CLIMAP estimates of the SST anomaly between the LGM and the pre-industrial

The MARGO project, a more recent compilation of global SST proxies for the LGM, found more extensive tropical cooling, typically of 1-2°C, with an east-west gradient in anomalies that climate models do not replicate.

MARGO estimates of LGM tropical SST

MARGO estimates of the tropical SST anomaly between the LGM and the pre-industrial

Some other proxy data suggest that the tropics were much cooler than reconstructed by MARGO or CLIMAP. Strontium/calcium ratios of corals suggest SSTs were up to 6°C cooler, which might imply a much higher climate sensitivity than the 2.4°C (1.4-3.5°C) estimated using the MARGO data.  However there are few coral data, perhaps partly because the oceans were too cold for widespread coral growth, and partly because most LGM corals are difficult to sample beneath 120 m of water).

This discrepancy has interested me and is part of the motivation for the reanalysis of the MARGO foraminiferal assemblage data I have started.

A paper in a recent issue of  Nature Geosciences by Tripati et al adds some more data that suggests that CLIMAP and MARGO underestimated tropical cooling. Tripati et al use clumped isotope palaeothermometry to reconstruct sea surface conditions off Papua New Guinea since the LGM. Clumped isotope palaeothermometry is a relatively new stable isotope method. Traditional methods have used the ratio of 18O to 16O in calcite, for example from foraminifera tests, but this ratio is sensitive to both temperature and the isotopic composition of the source water, which complicates interpretations. Clumped isotope palaeothermometry uses the temperature sensitive-tendency of the heavier carbon and oxygen isotopes to clump together in the same calcite molecule. Carbonate molecules with a heavy carbon atom and an heavy oxygen atom,  13C18O16O, are more common than expected by chance, but become rarer as the temperature at which the calcite formed increases. About half a percent of oxygen atoms are 18O, making the ratio of 18O to 16O much easier to measure than the excess of 13C18O16O, which is present in 44 ppm of carbonate molecules. Importantly, clumped isotope palaeothermometry, abbreviated as Δ47 (the sum of the isotopic weights), is independent of the source water composition and seems to show few species-specific vital effects.

Tripati et al find that the SST of Papua New Guinea was 4-5°C cooler than modern, and that the depression of LGM snowlines in the mountains of Papua New Guinea by almost a kilometre is not consistent with the smaller SST cooling reconstructed by MARGO and by many of the climate models.

Many sites and samples were used to develop statistically robust regional estimates. Values shown are calculated to be consistent with MARGO (ref.2), which reported Δ SSTs= mean annual SSTclimatological −LGM proxy SSTLGM (Δ47 difference is 4.1 ± 0.4 °C (2 s.e.m.), SSTcoretop −SSTLGM is4.7 ± 0.8 °C (2 s.e.m.)). 1 s.e.m. Δ47, Mg/Ca, alkenone, transfer function error bars shown (Δ47 black line = 2 s.e.m.); Sr/Ca: range. Model ΔSSTs (refs 9, 10) show LGM − pre-industrial annual mean (15° N–15° S, 140–165° E); vertical bars indicate spatial variability at core locations. b, LGM temperatures are 25.3 ± 0.4°C (2 s.e.m.). Temporal variation is similar to Mg/Ca records4. c, Ice volume-corrected δ18Oseawater.

Figure 2 from Tripati et al a) Proxy and model estimates of the anomaly between the LGM and modern SST. Note the sign of the anomaly is different from that given by MARGO and CLIMAP.
b, LGM temperatures are 25.3 ± 0.4°C (2 s.e.m.). Temporal variation is similar to Mg/Ca records.
c, Ice volume-corrected δ18O seawater.

If the work of Tripati et al is supported by new reconstructions, and perhaps by a reanalysis of the MARGO foraminifera assemblage data, then questions will need to be asked about the climate model performance. Any discrepancy could be because the models have underestimated climate sensitivity, alternatively the increases in dust and aerosols at the LGM, which have not been incorporated into the climate models, may have helped cool the climate. If the dust can explain the discrepancy, then perhaps current estimates of climate sensitivity do not need to be greatly increased.

Tripati et al 2014. Modern and glacial tropical snowlines controlled by sea surface temperature and atmospheric mixing. Nature Geoscience, 7205–209 doi:10.1038/ngeo2082

Posted in climate, Novel proxies, Peer reviewed literature | Tagged , , , | 2 Comments

Keenan’s accusations of research misconduct

NIH, the National Institutes of Health, defines research misconduct as

fabrication, falsification, or plagiarism in proposing, performing, or reviewing research, or in reporting research results

with the proviso

Research misconduct does NOT include honest error or differences of opinion

The emphasis is in the original. Other funding agencies have similar definitions.

Some climate change sceptics would appear to believe that a different definition should apply to climate scientists: being a climate scientist.

Allegations of research misconduct have been made against several climate scientists. Michael Mann and Phil Jones have, together, endured at least five inquiries into academic misconduct and been exonerated by each. I’m not aware of any climate scientist who has had an inquiry into academic misconduct make material findings against them.

Perhaps it is the hope that they will eventually that drives Doug Keenan, who has just submitted his fourth (I believe) formal allegation of research misconduct (in addition to several other claims of misconduct that he has not pursued formally – including against the Met Office). His previous allegations against Professor Wang (University at Albany) for his work on urban heat islands and Professors Manning (University of Reading) and Kuniholm (Cornell University) for their work on dendrochronology, were dismissed. This time he is making allegations of misconduct against Christopher Bronk Ramsey on the Bishop-Hill blog and reporting that he has submitted the allegations to Oxford University. Keenan’s allegations of misconduct stem from Bronk Ramsey’s radiocarbon work which Keenan claims is in error.

Keenan takes issue with two aspects of Bronk Ramsey’s work. The first is in how radiocarbon dates are calibrated. Keenan has argued in a paper published in Nonlinear Processes in Geophysics that the standard method of calibrating radiocarbon dates is wrong. Bronk Ramey continues to use the standard methods, as implemented in his OxCal software. Keenan believes this is grounds for a misconduct complaint – I fail to see how it can be categorised, at this stage, as more than “differences of opinion” even if Keenan was correct. And he certainly is not.

I wrote about the problems with his paper last year. The basic problem with Keenan’s proposed method is that it assumes that all radiocarbon ages are equally likely – they are not. Radiocarbon dates that correspond with plateaux in the radiocarbon calibration curve are more likely than dates that correspond to steep parts of the curve. Today, others have explained the same in the comments  at Bishop-Hill (see comments by Radford Neal and Nullius in Verba) . Keenan has yet to either accept that he is wrong or to try to demonstrate that he is correct.

The second allegation concerns the correct method for combining radiocarbon dates. Keenan believes that Bronk Ramsey used the wrong method in some papers and asks if corrigenda will be published. Despite Keenan’s combative email, Bronk Ramsey’s response is constructive and explains that the choice of model is “ultimately is a matter of opinion”. Neither in the email nor the paper does Keenan demonstrate that the choice of model makes a material impact on the results.

Keenan’s allegations of misconduct against Bronk Ramsey are hopeless – both complaints can be ascribed to “differences of opinion” and thus be dismissed. Keenan must know this - some of his previous complaints were dismissed on the same grounds. Any honest errors are Keenan’s.  These allegations, like the allegations against Mann and Jones before, will consume the time of the researchers and the committee assembled to investigate the complaint. When his complaint is inevitably dismissed, Keenan will no doubt complain about biases in the system, and other climate sceptics will nod vacuously.

Given its likelihood of success (almost exactly zero), it is difficult to see Keenan’s current complaint as anything other than vexatious.  Hopefully, he will realise there are errors in his calibration procedure, withdraw the complaint, issue a corrigendum to his paper, and apologise to Christopher Bronk Ramsey. Well one can hope.

Posted in Fake climate sceptics | Tagged , , | 6 Comments

The depressingly bright green lawns of Tucson

As a boy, I kept cacti. A hook-spined Mammillaria, an Opuntia with tufts of small orange spines which became near-invisible once embedded in my delicate fingers, and a few others shared the greenhouse with tomato and cucumber plants and my father’s collection of Streptocarpus. The cacti were sparingly watered in the vain hope that they would flower. I remember my surprise when some cacti seed I had sown sprouted with spineless cotyledons; I think the mould killed them before the wonderful yet monstrous morphologies of the adult plants were revealed. Outside the greenhouse, a lawn grew fast and green, watered by the copious Cumbrian rains.

You cannot do this with a lawn

You cannot do this with a lawn

The rains are conspicuously less copious in Tucson, Arizona. Most of the gardens reflect the desert environment with planting of drought-tolerant plants, magnificent tall saguaro, opuntia, yucca and aloe, set in gravel. But some of the gardens deny their location. Bright green lawns, fertilised and watered, stretched down to the pavement.

Never has a lawn looked more ridiculous

Desert denial in Tucson

Much of the western US is in the grip of a drought, yet in Tucson some people are pouring precious water on grass. These green patches have no obvious utility; they need frequent mowing, weeding and watering. If these people of Tucson are not prepared to adopt more sustainable gardens, which would save them money and inflict no loss of utility, how much harder is it going to be to persuade them to make other changes towards a sustainable lifestyle that will require even a small loss of utility?

US drought severity March 2014

Compared with tackling greenhouse gas emissions, tackling wastage of water should be easy. The time scales are short, severe drought may be affecting the south western US within a few months, rather than the multidecadal timescales for the worst of the effects of global warming to become apparent. The problem is local: what the Chinese do in tackling their water shortage problems are irrelevant to the problems in the US and so cannot be used as an excuse to do nothing. There are no feedbacks or complex physics for drought deniers to obfuscate. And some of the solutions are easy: desert gardens or astroturf; or at least not sprinklers on the lawn at noon; covering swimming pools, low water use toilets.

California is also in the grip of drought. In LA sprinklers keep the precious green from turning gold and men hose down the pavement. No hosepipe bans here, no infringement on the personal liberty to fritter away resources, but yet people are happy to, by order of the city council, abandon the sacred right to reverse park.

No profit in reverse parking, so nobody lobbies to protect our rights

No profit in reverse parking, so nobody lobbies to protect this right

Posted in climate | Tagged , , | Leave a comment