The Humpty Dumpty theory of palaeoecology

“When I use a proxy,” Humpty Dumpty said, in rather a scornful tone, “it means just what I choose it to mean – neither more nor less.”

“The question is,” said Alice, “whether you can make a proxy mean so many different things.”

With apologies to Lewis Carroll

Just as I was beginning to run out of things to write about in the chironomid reconstruction from Lake Żabińskie (I still want to write about figure 2 from Larocque-Tobler et al and ask why it is missing 18 lakes), along comes another chironomid reconstruction from Lake Żabińskie, this one spanning the last 1000 years (Hernández-Almeida et al 2016).

You might expect that a 1000-year August air temperature reconstruction would be published in its own paper, as it is sure to be included in compilations of palaeo-climate reconstructions, but no, it is tagged onto the end of a multiproxy study where most of the time the chironomid stratigraphy is used as an anoxia indicator.

Chironomids are of course sensitive to anoxia (most oxygen breathers are), but it is not clear why the authors believe that the first principal component of the chironomid stratigraphy is an anoxia indicator. They write

In Lake Żabińskie, between AD 1896 and 2011, the chironomids PC1 had a relationship with changes in anoxic taxa

but don’t inform the reader that the relationship is very weak (Pearson’s correlation coefficient = 0.23) and non-significant. Not much of a relationship. The second PC axis, in so much as you would trust any ordination axis of these data, has a somewhat higher (0.37) and significant correlation.

Since the chironomid assemblage is a mixture of littoral chironomids (73%) that live in the oxygenated water above the thermocline and profundal chironomids that live in deeper, perhaps oxygen-depleted water, the relationship between anoxia and chironomid assemblages over time is likely to be complex and noisy.

1000yr reconstruction

Chrysophyte-inferred winter temperatures and chironomid-inferred summer temperatures

Despite the importance of the chironomid record, the stratigraphy is not shown and no reconstruction diagnostics are given. There is no way for the reader to evaluate whether the reconstruction is any good (it is possible that some of this essential information is included in the supplementary online material which does not seem to be available at the moment).

I do hope that the authors have used the correct version of the chironomid stratigraphy,  and that the data will be archived promptly. I wonder what the count sums were?


1000yr PC1.png

First principal component of chironomids, diatoms and chrysophytes. Note the poor correspondence with the reconstructions.

Hernández-Almeida et al also include a chrysophyte-based reconstruction of winter temperature (first presented in Hernández-Almeida et al 2015a), but spend most of the paper using the first principal component of the chrysophyte data as a nutrient indicator. In contrast, Hernández-Almeida et al 2015b, interpret the same chrysophyte stratigraphy as an indicator of calcium concentrations and relate this to May-October zonal wind speed. These two previous papers by Hernández-Almeida et al make no attempt to acknowledge, still less address, these radically different interpretations of the same data. The winter temperature reconstruction follows previous literature and seems to have been planned in advance, whereas the zonal wind reconstruction reads like a fishing expedition.

Of course I understand the motivation for producing more that one reconstruction from a proxy record that took months of microscope time to generate. And of course I fully accept that biotic proxies are influenced by multiple environmental variables and that at different sites different variables might be possible to reconstruct. The problem is in knowing which. This is why well designed palaeoecological experiments are so powerful.

Reconstructing multiple variables from the same data inevitably means that the assumption of transfer functions that

Environmental variables other than the one of interest have negligible influence, or their joint distribution with the environmental variable does not change with time.

must be violated. Juggins (2013) showed that ecologically important secondary environmental variables can severely bias reconstructions. The chironomid anoxia response will bias the temperature reconstruction and vice versa; the chysophyte winter temperature response will bias the May-October zonal wind speed record and vice versa.

Considerable work is needed to demonstrate that a proxy can be used to make meaningful reconstructions of a single variable. Far more is need to make a convincing case for multiple reconstructions: Hernández-Almeida et al don’t even try.

“The question is,” said Humpty Dumpty, “is it publishable – that’s all.”

Posted in Peer reviewed literature, transfer function | Tagged , , , | 1 Comment

Nei, de er ikke geiter

P1050069I knew they were not ibex or ibis, so I assumed they were goats, high (at least it felt high) in the Tatra Mountains. Later I was asked if I had seen any chamois. No. But I I had seen… Wait, what does a chamois look like?

The Tatra chamois is critically endangered, but the population has recovered substantially over the last couple of decades.

These are what they probably eat.


Posted in Uncategorized | Tagged | Leave a comment

Mountain detour on the way to Białowieża

I am on my way to Białowieża Forest in eastern Poland to meet with collaborators on the KlimaVeg project and start analysing some of the project data. I have taken a detour to the Tatra Mountains, a compact calcareous range on Poland’s southern border.

I am used to the English Lake District, but the Tatras are much busier. Zakopane (which has night trains from the far corners of Poland), makes Bowness look like a ghost town.

Posted in Uncategorized | Tagged , | Leave a comment

Replicating the Lake Żabińskie reconstruction

Having shown that the archived chironomid data from Lake Żabińskie are strange in several ways and have bad reconstruction diagnostics, I want to see how well I can replicate the August air-temperature reconstuction.

First, I need to build the transfer function. I’m using WAPLS on square-root transformed species data. All taxa with more than one occurrence are included, and the lakes with declared low count sums are omitted.

keep <- colSums(spp > 0) > 1
mod1 <- crossval(WAPLS(sqrt(spp[, keep]), env), cv.method = "bootstrap", nboot = 5000, verbose = FALSE)
knitr::kable(performance(mod1)$crossval[1:3, 1:4], digits = 2)
RMSE R2 Avg.Bias Max.Bias
Comp01 2.68 0.65 -0.06 9.97
Comp02 2.33 0.76 -0.04 6.75
Comp03 2.38 0.76 -0.08 6.27

The bootstrap cross-validation performance is very similar to that reported by the corrigendum (WAPLS component 2: RMSEP = 2.3°C, r^2^ = 0.76).

The reconstruction is similar to that archived but not identical.

wapls.sqrt <- predict(mod1, sqrt(fos))$fit[, "Comp02"]
reconstruction <- data.frame(
chron = chron,
Instrumental = instrumental$temperature,
Archived = recon$temperature,
Replication = wapls.sqrt
reconstruction2 <- gather(reconstruction, key = "Series", value = "Temperature", -chron)
ggplot(reconstruction2, aes(chron, Temperature, colour = Series)) +
geom_line(alpha = 0.5) +
scale_colour_manual(limits = c("Instrumental", "Archived", "Replication"), values = c("black", "red", "blue")) +
labs(x = "Year CE", y = "Temperature °C", colour = "Series")


The mean of the archived and replication reconstructions are very similar (16.92 vs. 16.88), but the variance of my replication is about 20% higher (2.19 vs. 2.67).
One of the most noticable differences is that the archived reconstruction has a value for 1925, whereas the chironomids have data for 1927. This 1925/1927 switch also occurred during the evolution of the fossil data.

The other differences between the archived reconstruction and my replication might be because of different sites included in the calibration set (LT15 omit nine lakes on the basis of a PCA, but its is not clear which lakes these are and whether they are also omitted for the corrigendum) or different species inclusion rules. The bootstrap that LT15 use will also cause variability.

Ideally anybody who has the raw data should be able to replicate the results of any paper exactly. Given the vague description of the methods in LT15 and the corrigendum, the replication is as good as can be expected.

This is probably my last post detailing oddities in the data archived by Larocque-Tobler et al (2015). I think I have done enough to demonstrate that the data have unexpected properties that need explaining. In my next couple of posts in this series, I’ll describe my quest to get that explanation, or how I mistook a yo-yo for a die.

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

Low calibration set counts in LT15

Larocque-Tobler et al (2015) promised that, for their chironomid calibration set

  At least 50 head capsules (Heiri et al., 2001 and Larocque et al., 2001) were mounted.

A cursory examination of the data shows that this is unlikely to be true: many of the calibration set samples have a minimum relative abundance above the 2% expected for a singleton in a count of 50 heads.

The corrigendum admitted as much

Lakes with low numbers of head capsules were included in the transfer function used to create the published temperature reconstruction… The following lakes have now been deleted for the corrected temperature reconstruction: GOR, KOS, LEK, SAL, SZE, SZOS, TRZ, WAS and ZAB.

If a few samples were a few (i.e. can count on fingers of one hand) chironomids short of the promised fifty, I would not object. In the Polish calibration set, nine of the forty-eight lakes (two of the fifty lakes sampled have no chironomid data), are now reported to have fewer chironomids than promised, but how many fewer?

Since most samples should have at least one singleton, we can estimate the probable count sum from the reciprocal of the relative abundance of the least abundant taxa.

This plot1 shows the minimum percentage and the probable count sum for the Polish calibration set, highlighting the lakes that the corrigendum declared to have low counts.


Two aspects of this figure are striking.

First that some of the probable counts are very low. Lake GOR appears to have a count of just five heads, based on a minimum relative abundance of 20% (the missing 20% might be an unidentified chironomid).

Cricotopus Endochironomus albipennis Microtendipes pedellus
GOR 40 20 20

Second, in addition to the reported nine Polish lakes with counts below fifty head capsules, a further ten have a minimum relative abundance greater than 2%, indicative of a count lower than fifty heads.

What about the Canadian portion of the calibration set? There are 73 Canadian lakes (not 72 as reported in LT15 – the data source, Larocque (2008), also reports 73 lakes). Four of these lakes appear to have a minimum count below fifty.


This possible undercount in a few samples in Larocque (2008) is interesting but will not materially affect the results of that paper.

There are some additional curious aspects to the sums of the relative abundances for each lake in the calibration set. This sum can be below 100% if rare or unidentified taxa have been omitted. Values slightly above 100% are possible with rounding errors. But I struggle to understand how a sum of 105.1% is possible for Polish lake GIL. The sum of 175% for Canadian Lac H is beyond comprehension.

  1. All code used in this post is available on github 
Posted in Peer reviewed literature, Uncategorized | Tagged , | 1 Comment

The missing rare taxa at Żabińskie

“rarity is the attribute of a vast number of species of all classes, in all countries.” Charles Darwin

In any census of any species-rich community, the rarest taxa are likely to be represented by a single individual. This can be visualised with rank abundance curves, shown here for Barro Colorado Island tree counts (jittered for clarity).


In the BCI data, every sample contains many taxa represented by a single individual. Across all samples, 38% of taxon occurrences are of a single individual.

This is not the case for the Lake Żabińskie chironomid counts where 34 of 89 samples lack any taxa represented by a single (or a half) head.


Not only do many samples lack singletons, one lacks any taxa occurring with fewer than five head capsules. There is a curious trend in the proportion of samples without singletons, with a higher prevelance of such samples in the second half of the record.


These are the counts for the two most extreme cases.

pos(fos_counts[minc == 5, ])
##      Chironomini Microtendipes pedellus Cladotanytarsus mancus1
## 1967           7                      6                       6
##      Tanytarsus sp Tanytarsus lactesens Tanytarsus mendax
## 1967            11                    5                 6
pos(fos_counts[minc == 4.5, ])
##      Polypedilum nubeculosum Cladotanytarsus mancus1 Paratanytarsus
## 1970                     4.5                       8             12
##      Corynoneura Cricotopus Nanocladius branchio Orthocladius
## 1970           5          5                  4.5          4.5


How unlikely is it to have so many samples without singletons?

One way to test this is to fit a rank abundance model to the count data and then simulate assemblages from the model. There are a variety of rank abundance models that can be fitted with the vegan package. Here they are fitted to the first BCI sample.

rf <- radfit(BCI[1, ])


Few singletons are expected with the pre-emption model, so I am going to apply this model to the chironomid counts and then simulate assemblages and test if they have singletons.

fos_counts1 <- fos_counts
# promote half heads to full
fos_counts1[fos_counts1 %% 1 == 0.5] <- fos_counts1[fos_counts1 %% 1 == 0.5] + 0.5
#fit model and extract coefficients
alpha1 <- apply(fos_counts1, 1, function(r)  coef(rad.preempt(r)))

simulateMinAbun <- function(alpha, maxRank = 50, n = 10000, J = 30){
#maxRank is number of taxa considered, J is number of individuals, n is number of trials
rank <- 1:maxRank
abun <- J * alpha * (1 - alpha)^(rank - 1)
sims <- rmultinom(n = n, size = J, prob = abun)
table(apply(sims, 2, min0))/n
high20 <- simulateMinAbun(max(alpha1), J = 20)
high30 <- simulateMinAbun(max(alpha1), J = 30)
high70 <- simulateMinAbun(max(alpha1), J = 70)

med30 <- simulateMinAbun(median(alpha1), J = 30)
phigh <- pbinom(q = sum(minc > 1), size = nrow(fos_counts), 1 - high30[1], lower.tail = FALSE)
pmed <- pbinom(q = sum(minc >  1), size = nrow(fos_counts), 1 - med30[1], lower.tail = FALSE)

For the sample with the largest coefficient for the pre-emption model (i.e. the steepest slope of the rank-abundance relationship), the probability of the count having at least one singleton is 0.91. This result is not sensitive to the number of heads counted (at least over the range reported for the chironomid data).

With the median pre-emption model coefficient, the probability of the count having at least one singleton is 0.9984.

Even in the most generous case, the probability of having the observed number of samples without singletons is 5 x 10-14. With the median case, the probability is 9 x 10-74.

This calculation does not account for the taxon inclusion rules (which were not adhered to) which would have removed taxa occurring in fewer than three samples. This would adjust the probability by several orders of magnitude. However, the calculation does not consider that in some samples neither singletons nor doubletons occur, which would adjust the probability by orders of magnitude in the other direction.

Whatever adjustments you want to make, the lack of rare taxa in the chironomid data is remarkable.

In the first version of the dataset, samples apparently missing singletons (identified by the minimum percent being above 2%) tended to have all their values as integer multiples of the minimum percent. This strongly suggested (as is now admitted) that the counts were not fifty as claimed, but much lower for some samples. Since then, the data have “evolved” and this integer multiple only holds for a few samples.

It is not obvious why there are so many samples without singletons. Perhaps the taxon inclusion rule (“3 percent in at least 3 samples”) was misapplied, smiting singletons in some levels only until the observed pattern emerged. At the calculated probabilities, almost anything is more likely than the data being correct.

It is trivial to infer that the currently archived data are definitely not the original data. The original data (all 76 taxa) needs to be archived.

Posted in Peer reviewed literature, Uncategorized | Tagged , | 1 Comment

ggplot2 maps with inset

Lots of students want maps showing their field sites for their thesis. Several have come to me with code they found on the internet but couldn’t get to work. The problem is that the plotting package ggplot2 has evolved since that code was written and some things no longer work (for example opts instead of theme), ggplot2 now supports map projections and Norway is more complicated than Canada in that it has remote territories.

This post updates the code from code to make a map with an inset showing the location against a regional map.

I’m going to plot the seedClim locations.

Load the required packages and some data.


dat siteID           latitude longitude Temperature Precipitation
Alrust           60.8203    8.70466           2             1
Arhelleren       60.6652    6.33738           3             3
Fauske           61.0355    9.07876           3             1
Gudmedalen       60.8328    7.17561           1             3
Hogsete          60.876     7.17666           2             2
Lavisdalen       60.8231    7.27596           1             2
Ovstedal         60.6901    5.96487           3             4
Rambera          61.0866    6.63028           2             3
Skjellingahaugen 60.9335    6.41504           1             4
Ulvhaugen        61.0243    8.12343           1             1
Veskre           60.5445    6.51468           2             4
Vikesland        60.8803    7.16982           3             2
dat$Temperature dat$Precipitation precipLab tempLab ```

Find the latitude and longitude ranges required for the high resolution map. Here I've added a buffer round the locations.

xlim ylim ```

### Get the map data

# Low resolution map of Norway for the inset map

# High resolution map of Norway for the main
norwaymapHires ```
Unfortunately, even the high resolution map is not high-resolution enough - lots of important islands off the Hordaland coast are missing, for example Askøy near Bergen. I need to find a higher resolution map.

### Make the maps

The inset map showing all of Norway

a geom_map(data = norwaymap, aes(x = long, y = lat, map_id = region),
map = norwaymap, colour = NA, fill = "grey60") +
geom_rect(data = data.frame(),
aes(xmin = xlim[1], xmax = xlim[2], ymin = ylim[1], ymax = ylim[2]),
colour = "red", fill = NA) +
coord_map(xlim = c(3, 33), ylim = c(57, 72)) +
labs(x = NULL, y = NULL)

#print(a)# print this map until you are happy with it.

geom_map draws the map.

geom_rect shows the limits of the high-resolution map.

coord_map projects the map. Any of the projections in the mapproj package can be used – there are probably better ones than the default Mercator projection. I had to set the x and y limits to exclude Svalbard, Jan Mayen and Bouvetøya (a subantarctic island).

Zoomed in map, with the site locations

b shape = Temperature, fill = Precipitation)) +
geom_map(aes(x = long, y = lat, map_id = region), data = norwaymapHires, map = norwaymapHires, colour = NA, fill = "grey60", inherit.aes = FALSE) +
geom_point(size = 3) +
coord_map(xlim = xlim, ylim = ylim) +
labs(x = NULL, y = NULL) +
scale_shape_manual(breaks = 1:3, labels = tempLab, values = c(24, 22, 25)) +
guides(shape = guide_legend(override.aes = list(fill = "black")))+
scale_colour_brewer(breaks = 1:4, labels = precipLab, palette = "RdBu") +
scale_fill_brewer(breaks = 1:4, labels = precipLab, palette = "RdBu")

I’m using various scale_?_??? commands to set the shape and colour of the points. I need to use the guides function as otherwise the temperature legend has empty rather than filled shapes.

A theme

maptheme axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
panel.border = element_rect(fill = NA, colour = "black"),
panel.background = element_blank()

This theme can be added to the plots. It sets a dark rectangle round the map and removes the grey background and axis annotation. An alternative would be ggthemes::theme_map() but it moves the legend.

Printing the maps

Print the two maps together, one an inset of the other using the grid package for flexible alignment of ggplot graphs.

vp_b vp_a print(b + maptheme, vp = vp_b)
print(a + maptheme, vp = vp_a)


Posted in R | Tagged | 7 Comments