Climate envelopes for Madagascan butterflies

Cite as: Lees, D.C. 2002. Modelling climate envelopes for rare Madagascan butterflies. Final report with atlas. Unpublished report to Conservation International, Washington.Biogeography and Conservation Laboratory,Dept. of Entomology,Cromwell Rd.,S. Kensington,SW7 5BD

Modelling climate envelopes for rare Madagascan butterflies

Prepared by: Dr D.C. Lees


Despite mounting knowledge of the importance of Madagascar’s highly endemic biota, there are as yet no published examples of bioclimatic modelling for suites of species on the island, and so no assessments have been possible of the potential impact of climate change. Now that high resolution datasets are emerging, whilst habitat loss continues apace in an uncertain political climate (e.g. Bradt, 2002), time is ripe for such studies. Evaluating how organisms may respond to the environments and landscapes of the near future is vital to inform conservation and restoration planning, and ultimately to avert an extinction catastrophe on the island. Bioclimatic and habitat models can be applied to future climate and forest fragmentation scenarios, as soon as they become available.

Butterflies are a particularly appropriate group of organisms to examine, because they are highly sensitive to climate change (work done in temperate regions: Parmesan, 1996; Hill et al., 1999). Moreover, it is likely that butterflies are somewhat more representative of overall biodiversity (as most species are insects) than, let us say, birds or mammals and are likely to show more fine-grained responses to the environment (Kremen, 1992). As a group that shows a wide range of dispersal abilities, but with the majority of species forest-restricted, Madagascan butterflies constitute an excellent group to examine the potential to track future climate space in the context of habitat fragmentation models. On Madagascar, butterflies exhibit a 74% level of species endemism when both described and undescribed species are considered (Lees et al. in press), increasing to 100% in some large evolutionary radiations on the island (Raharitsimba, 1997). Also, their alpha-taxonomy has been subject to recent (Lees, 1997) and ongoing (Lees et al. in press) revision. Not least, a high quality point dataset for Madagascar butterflies has recently been compiled that builds on previous distributional models for all ca. 300 described species on the island (Lees et al., 1999).


  • To validate available climatic surface models for the Madagascar biodiversity hotspot
  • To evaluate techniques for producing species bioclimatic models
  • To model present day climate envelopes using the selected technique for a set of Madagascar butterflies of potential conservation concern
  • To explore how such models might be applied to future climate change in the context of anthropogenic habitat fragmentation

These four objectives are refinements of the work that was proposed in Phase I of the proposal of Lovett and McClean (2001). Although originally proposed to model a few species, it has proved possible here to model 48 species (1/6 of the described butterfly fauna, including many of the rarest species). This report is self-contained, embedding a representative sample of the climate surfaces and a complete atlas of all the species model surfaces produced for this pilot study.

Existing bioclimatic modelling tools

There are an increasing range of available methodologies and programs for bioclimatic modelling which can be used to produce bioclimatic envelopes. The current list is based only on those come across during this work, and is far from exhaustive. These include:

(1) Bioclim (Busby, 1991; Nix, 1986l)
(2) Domain (Carpenter et al. 1993);
(3) Distmodel (Skov, 2001; Skov and Borchsenius, 1997)
(4) GARP and its successor WhyWhere (see also Stockwell and Noble, 1992)
(5) Species Range Analyst (dedicated software for Madagascar; not yet public domain);
(6) Floramap (Jones and Gladkov, 1999)

Each tool includes some advantages and disadvantages, and some were very briefly (and thus probably unfairly) compared during this work. Species Range Analyst (hitherto known as Platforme d’Analyse) was extensively used in this work for visualising species ranges based on points, and this toolset can make predictions based on both grid surfaces such as Digital Elevation Maps and categorical shapefiles such as the Cornet (1974) dataset.

However, because this extension to Arcview presently forces a convex hull around the points, it is not yet suitable for climate envelope modelling. Bioclim was beyond our budget for this work. Skov’s program was freely available and also provides an interface with Arcview, but has not yet been properly evaluated for Madagascar work; like Domain, it includes a similarity-based approach (see also botanik.uni-bonn), but also allows for standard climate envelope mapping and includes some nice distance-weighting features.

I examined the program Domain32.exe, which I found too restrictive by not allowing evaluation of best fit surfaces and similarity matrices (a black box problem), although it will accept points. The GARP approach has recently been used in a major paper predicting climate for a large number of Mexican organisms (Peterson et al., 2002), but when I evaluated WhyWhere on the web on a point file, the program took some time to return an apparently uninterpretable image that was evidently a composite of numerous surfaces with no indication of the importance of each. As GARP is an extremely complex algorithm, and may also suffer from the black box problem, it was decided to use a simpler approach. Initially, I tested a custom-built genetic algorithm.

Genetic algorithm

Customised for use of this project was a genetic algorithm programmed in C++ by Colin McClean, which I attempted to use in conjunction with the Madagascar climate grids. In its first incarnation (GAFIN4.EXE) the program was a little tedious to use, as it required keyboard input of parameters (a temporary workaround was found by writing DOS batch files to call the program, which worked). In its second incarnation (LARGERGA.EXE) the program was modified to call its own parameter file. This was a great time-saver, and moreover the program was expanded to accept up to ten, rather than four, grids at a time (however, this still meant that a subjective preselection of ten grids out of twenty was necessary). Whilst GARP (mentioned above) provides a considerably more sophisticated genetic algorithm incorporating rule set selection, GAFIN was far more useful to us as it is far more focused in its application. Because it is useful to understand for its potential importance for climate modelling, the genetic algorithm approach is reproduced below from Lovett and McClean (2001, unpublished proposal).

The genetic algorithm starts with a “population” of possible climate variable ranges, known as “chromosomes”. To speed the algorithm, the starting values can either be the possible extremes for the total area under consideration or the possible extremes for the total extent of occurrence of the species being modeled. For example, one chromosome from the population might have a minimum value for mean maximum temperature of 25 oC, a maximum value for mean maximum temperature of 32 oC, a minimum value for mean annual precipitation of 600 mm, a maximum value for mean annual rainfall of 850 mm and so on, until ranges are defined for each of the climate variables included in the analysis. In this case, each minimum or maximum value represents one “gene” on a chromosome. The areas accounted for by the set of variable ranges contained in the genes can be mapped easily and compared to the species range identified. This gives a score, or fitness value, of how well the variable ranges given by that chromosome of the population account for the spatial distribution of the species range. The starting population could be made up of any number of individual chromosomes, but a prototype algorithm developed on comparable data sets performed well with 100 individuals. The maximum and minimum values for the variable ranges, the genes, are allocated randomly to the initial population.

On the basis of their fitness values (see below), members of the population are allowed to reproduce. That is, if a chromosome has a high fitness score it has a higher probability of reproducing with another chromosome from the population. Reproduction occurs by taking a number of genes from one chromosome and a number of genes from the other chromosome. This is known as a cross-over operation, where two offspring are produced. The “parent” chromosomes are split at random along their lengths. The portion of genes before the split on the first chromosome provides the first genes for the first offspring, the portion after the split on the second chromosome provides the rest of the first offspring’s genes. The second off-spring is created from the first portion of genes from the second parent chromosome and the second portion of genes on the first parent chromosome. In this way “good genes” are taken from successful individuals into the next generation in a search for the optimal solution.

However, the genes of the initial population may not be the genes that could provide the best possible fit between the climate envelope and the actual distribution points. In an attempt to avoid local maxima resulting from the limited nature of the genes of the initial population an operation known as “mutation” is added to the algorithm. A mutation operation simply allows the value of a randomly selected gene on a randomly selected chromosome to be changed by a random number. Mutations can survive into future generations if the fitness scores of the mutated chromosomes are high enough to give them a good probability of reproducing.

The algorithm advances through generations of solutions until no noticeable improvements are made to the fit of the climate envelope to the species distribution points. This heuristic optimisation approach has certain advantages over more traditional statistical approaches to creating a predictive model of species presence/absence. Notably, possible regression techniques such as logistic regression may be affected by over-dispersion caused by model mis-specification. This can result from the spatial autocorrelation found in the climatic independent variables and the dependent variable itself. It also makes intuitive ecological sense to consider ranges of climate variable values that may be suitable for the occurrence of species, rather than finding single optimal variable values associated with species occurrence.

Fitness was measured at each step using the Dice-Sorensen index (see Magurran, 1988; basically the same as Jaccard’s index of similarity, but with congruence, the number of gridcells in common or shared, emphasised by multiplying by two). The basic Jaccard formula is:

Shared/(Shared + Omissions + Commissions)

where shared = # of species gridcells for validated points
and with omissions and commissions defined as in Eeley et al. (1999): commissions are cells predicted but not observed to be present, whereas omissions are observed to be present but not predicted. Note, however, that that for point or museum-type data, the terms “false positives” and “false negatives” respectively, are only relevant in terms of observations in the database, rather than in potential reality. In the case of the Dice-Sorensen index, the number shared is multiplied by 2, emphasising the importance of congruence. The genetic algorithm was designed for cases where entire areas (such as biomes) could be modeled. and in these cases pixel values can be found for any cell within the region using remote sensing. Prior to this work, the algorithm was not specifically tested on scattered point data (such as museum data) characteristic of data-poor areas such as Madagascar. In such a case, Omissions = 0 (if all points are taken into consideration in generating the envelope) and potential commissions are minimised to find fitnesses as close to 1 as possible.
After considerable experimentation with the required nature of the ArcInfo ascii format gridfiles in Arcview (no-data values and datavalues should be in the format -9999 and 1 respectively for this program to work), I found that LARGERGA runs on points as well as range fill polygons saved as ascii integer grids (temperature surfaces were multiplied by 1000 to retain floating points and avoid small value errors that occur otherwise (see also below for methodology), but:

(1) final values for all grid surfaces input are returned; none are rejected
(2) there is no way to evaluate one grid’s fit against another from the model output; the program just runs towards a final fitness value
(3) for points predicted areas may in some cases be vanishingly small (i.e., zero area)
(4) a purely cosmetic limitation – the output produced presently requires recoding in order to be usable in the Map Calculator of ArcView Spatial Analyst.

At the time of writing, none of these potential bugs had been fixed. It was therefore decided to abandon for present purposes the genetic algorithm approach, in favour of a somewhat simpler method further customised to requirements. Extra resources will be needed to support algorithmic advances including combined approaches (e.g. in Phase II of the proposed work). For this pilot work then, the following algorithm was used to model the rare butterfly species (for descriptive convenience this algorithm is termed a “point-driven tournament”), for which the same indices as were used by GAFIN/LARGERGA can be also be used to optimise fit (in this case, the simplest one was used, Jaccard). In practice this quite simply meant dividing the number of cells in the bioclimate envelope by the total cells in Madagascar and searching for the climate grid that provides the tightest solution.

Point-driven tournament

» Identify and remove any climatically extreme points across the sum of point values for all [20] surfaces examined whose location is uncertain under pre-defined criteria

» Find surface out of 20 for observed remaining point ranges that minimises commissions to greatest degree (# or % gridcells in prediction).

» Pair winner of (2) with each of 19 non-duplicate point climate ranges and identify new minimum.

» Pair winner of (3) with each of 18 non-duplicates and identify new minimum.

» Pair winner of (4) with each of 17 non-duplicates and identify new minimum…..etc., until commissions level off, or the requisite number of surfaces are reached

In practice, this approach has a number of clear advantages over many other approaches:

1. Outliers are identified and removed according to required criteria (high spatial error as coded in the database, suspect collector, old collection, expert decision, etc.) with respect to the surfaces examined (in this case, entirely climatic).

2. Any number of grids can be used (here, 20).

3. There is no prior subjectivity in the choice of surfaces. This is a great advantage in areas like Madagascar when nothing or next to nothing is known about the physiology of the species (and in the case of butterflies, there is the physiology of different stages to consider, apart from any requirements that the physiology of specific foodplants might impose on their own bioclimatic range).

4. Minimising potential prediction errors produces more conservative solutions that need less “explaining away”.

5. Each grid is evaluated on the basis of the total number of gridcells in the prediction.

6. The actual number of sampled cells for any taxon on average is roughly proportional to the surface area of the prediction. Initial tests show that area is a reasonable surrogate for actually sampled points – but the area, as opposed to the cells sampled approach, provides sharper discrimination amongst competing climatic explanations in cases of few points.

7. The grid resolution used being fine (20343 gridcells across Madagascar’s surface), there were few ties to deal with, and thus few bifurcations requiring extra processing time.

8. The algorithm tends to choose complementarity classes of data (e.g. moisture-related paired with temperature related-sets) at consecutive rounds. Although by no means a rule of thumb, closely spatially autocorrelated sets such as the DEM-modelled variables were not chosen consecutively, despite being available to the algorithm.

9. Not all surfaces fed to the algorithm have to be retained (for some real range of values or as a single limit) within the output from the algorithm. Species can be compared for their fit to n grids. In this case, it was decided that a maximum of 4 surfaces should be used, because the tightness of the fit tends to level off after this stage, and the more surfaces are used, the higher the chance that surfaces will be included that really have no biological meaning for the species. Note particularly that whilst point-driven “carpet bombing” using multiple climate surfaces still produces predictions outside known geographic range, explanations involving fewer surfaces are more parsimonious.

10. Observed reliable point minima and maxima are forceably retained (these are more likely to have some bearing to potential range-boundary constraining climatic contours, than a single statistical value)

The disadvantages of this algorithm include:

1. The final solution for n grids may be suboptimal, because the winner of one round may not necessarily be the optimal one when measured in combination with the consecutive round (non-duplicate candidates). It was found, encouragingly, that there is only a strong likelihood that the winner of one round is also optimal in combination with the winner of the next. Meanwhile, neither does use of a genetic algorithm on its own guarantee optimality. Potential solution: the algorithm can be adapted to search all possible combinations of grids, but total processing requirements rather quickly becomes unmanageable for non-heuristic approaches (less of a problem with efficient programming and where few grids are required in the solution).

2. It may be that only the minimum or maximum of one surface is important in determining (parts of) the species range boundary in reality (e.g. minimum temperature, rather than maximum). Potential solution: combining the above type of algorithm with a genetic one could provide a powerful solution which would further minimise potential prediction errors.

3. The algorithm needs to be scripted, most conveniently within an ArcView environment (for this analysis, it was implemented within Excel).

Strabena sufferti

The basic rationale underlying the point-driven approach for minimising commission errors is illustrated by Fig. 1. The distribution of points, shown in black for the rare satyrine butterfly Suffert’s Hook-tip, is here modelled by a climate envelope using a point-driven tournament for 4 grids. All points recorded for any of the Satyrinae (just an example of an appropriate higher taxon on which collectors might have concentrated) are shown (those outside the envelope in blue). A tight predictive envelope, as shown here, is much more conservative than a loose one in that in only needs to explain away, in this case, 7 cells (in red) where Satyrinae have been recorded and Strabena sufferti not. Because on average, the number of such cells will be proportional to area, we can use area of the bioclimatic envelope as a surrogate for measuring the fit of each surface to the points. This prediction does not, of course, take account of sampling intensity in those cells, but it could do. Note also that the prediction of absence by the climate envelope (all blue cells) is much more reliable than the prediction of presence, which in this case is less than 50% confirmed by the data.

Climatic datasets

In this study I examined three different datasets:
(1) the bioclimatic model created by Schatz and Weiss from the Cornet 1974
(2) two custom built models created specially for this project using thin plate splines and coastline clipping, from the Oldeman (1990) dataset for rainfall and temperature stations, and
(3) the climate models created by Hutchinson, known hereafter as the “CRES” dataset, a set of 20 nested climate surfaces at 0.05 degree resolution. For this evaluation, all datasets were converted to Arcview grids at the same resolution (0.05 degree, or about 5 km sided squares).

I evaluated the “Cornet” GIS model visually. Because pycnophylactic interpolation (e.g. Tobler, 1979) would be needed to make sense of this categorical dataset for continuous point values, and because the dataset would not provide enough independent climate surfaces (there is basically the dry-day and categorical maps) without much work, I rejected this dataset for use in climate envelope modelling.

Next, using a pre-compiled dataset from Oldeman (160 stations dotted across Madagascar for mean annual rainfall, and 115 stations for temperature), climate data range 1931-1970, I created two spline surfaces. These were created basically for cross-validation purposes against the CRES dataset, rather than in combination. Unfortunately, the precise climate station data on which the latter dataset is based is not precisely known to me, but may be assumed for now to be ca. 1920-1970, as for other Africa region data.

Butterfly dataset

One sixth of the butterfly fauna of Madagascar (50 species here using updated taxonomy: Lees et al. in press), including many of the species of most conservation interest and focusing on the substantial endemic rainforest radiations, were originally selected for initial climate modelling. Distribution point and polygon files and corresponding ascii grid surfaces with the same spatial resolution and extent as the climate data were created for all these species. For very few species originally selected, there remained insufficient distributional gridcells to allow distributional modelling. (e.g. Heteropsis wardi), or in one case (Charaxes phraortes), points fell just out to sea with respect to the climatic grid surfaces, and on an island outside the range of the climate data. In future then, the models should be extended to cover islands and coastal towns.

The area given for each species is that for convex hull of the points, and can be used to estimate rarity and conservation status using IUCN criteria of occupancy where other data is lacking (see Lees et al., in press). The elevational range limits are the observed ones for the data, rather than from a digital elevation map intersection. Note how significantly, even at the fine resolution of the analysis (5 km.), the number of actual points in the butterfly databases cut down.

Spatial and temporal coverage

To put the coverage of the butterfly dataset into spatial and temporal perspective, the sampling level for the full butterfly database is shown in Figure 2. Fig. 2A shows all records since 1949, when high quality data particularly from French collectors started to become available, and is compared with that since 1988 in Fig. 2B, when modern sampling (including extensive use of GPS) started.

Fig 2: Spatial and temporal coverage

This comparison makes the clear point that museum data cannot be ignored when performing climate modelling. for most groups of organisms in Madagascar. This means that to achieve best possible results, the butterfly data has to be compared with roughly equivalent climate models covering a similar and significant period. The CRES data is here considered to be roughly 1920-1971, but this should be checked. Considering that a substantial amount of butterfly data precedes 1949 (indeed good data emerges from some 1880’s collectors such as the Perrot Frères), 1950 would seem to be a reasonable average pivot date for comparison. It must be regarded as a substantial caveat, therefore, that stricter temporal partitioning has not been possible for this pilot study.


Climate surface cross-validation

The results of the climate surface cross-validation are shown in Fig. 3. Two CRES surfaces, mean annual rainfall and mean annual temperature, were validated with surfaces created during this work from the Oldeman (1990) dataset. An identical classification was used, with integer grids, at 0.05 degree resolution. Figs 3B and 3D were interpolated from the Oldeman data (using a tension spline, minimum 12 points, >98% correspondence in values with the original data points; temperature based on 115 climate stations; rainfall based on 160 climate stations; draped over the digital elevation map). Note that both CRES surfaces were originally modeled by CRES in ANUCLIM using the half-minute digital elevation map, but the Oldeman data has NOT been remodelled for the elevation map. White areas shown are spline over or undershoots, due to local absences of climate stations. Note that bands are broadly similar in extent for independent datasets, but the lack of high montane climate stations accounts for some of the anomalies particularly in the north of the island. There are also differences in the contours that cannot be addressed, without direct comparison of the baseline climate datasets used in this case, and remodelling of the Oldeman data using identical algorithms such as ANUCLIM.

Fig. 3. Cross-validation of climate datasets.
Here, two CRES climate surfaces are cross-validated visually using the Oldeman climate surfaces. MAR= mean annual rainfall; MAT= mean annual temperature.

It is reasonable to conclude from this fairly crude comparison that the CRES surfaces seem to be good working models, with the proviso that those variables do not come with metadata specifying the exact climate station data on which they were based. Significant differences in the maps in Fig. 3 suggest merely that somewhat different raw climate data were indeed used, but there is probably substantial overlap. Another important caveat to consider when using the CRES surfaces for butterfly modelling, is that they will show strong correlation with the elevation surface.

Therefore, during initial testing it was also important to examine the inter-relationships of the variables. The results of a correlation matrix analysis of the surfaces are summarised in Fig. 4.

Fig. 4. Inter-relationship of CRES climate variables and DEM.

The 20 CRES variables used for the modelling are compared in Fig. 4. The digital elevation map (DEM) is shown at the left, and mean annual temperature surface as an example to the right. Note the very close correlation coefficients for the temperature-related variables including potential evapotranspiration. These surfaces are strongly spatially auto-correlated with the DEM (and therefore also among themselves). This is not surprising since the DEM was clearly used in their modelling. Note also how the moisture related variables, to the left of the graph, are little correlated with the DEM, and little correlated among themselves. This implies that it would not be a good idea to preselect a number of temperature-related variables as grids for input to a genetic algorithm (for example), when a parsimonious explanation is required for the distribution of points. Fortunately, the method used here tackles this problem.

Exploration of approach to, and tests for, climate models

What measure of distribution to use?

During testing of the genetic algorithm, three type of inputs were investigated to represent the butterfly distribution, all exported as grids at the resolution and spatial extent of all climate surfaces for each butterfly surface
(1) filled convex hull of the point distribution
(2) filled minimal enclosing polygon of the point distribution and (3) the points themselves.

Algorithms for the first two solutions were obtained from the ArcScripts website on For the third output, the points were represented by the enclosing 0.05 degree grid-cell.

These results are not given in detail here, as in the end, due to shortcomings already mentioned for the butterfly modelling problem, the genetic algorithm was not used to produce the butterfly models presented here. Suffice it to say here that the polygon approach does not fairly represent the data, as is obvious from an examination of Fig. 5. Even the minimally enclosing polygon has edges that overlap certain climatic contours, outside any represented by the points themselves.

Fig. 5. Which representation of points to use?
1. Convex hull (black line)
2. Minimally enclosing polygon (red line)
3. Closest approximation to original points (white squares). Distribution is that of “Heteropsis” (Admiratio) paradoxa, illustrated, against spring mean temperature

Therefore, the filled polygon approach will invariably generate an overprediction (in many cases a very substantial one) of the bioclimatic envelope. The only possible approach is to use the points themselves within the algorithm, either with spatial error radii where they are encoded in the database using the Create Buffers option in ArcView Spatial Analysis (a variant also tried that also may cause similar problems where a more complex error shape is in fact appropriate), or as here at the minimal resolution of the grids. However, due to undetermined factor(s) it was empirically observed that the use of such points in the genetic algorithm led to under-prediction, one reason it was abandoned for the purpose of this study. For the case of points, it is necessary to choose either the minimum and maximum for each surface (raising the problem of how to deal with outliers), or some measure of optimality plus or minus some statistical deviation (raising problems as to which if any statistical deviation is appropriate, and meaning that species with very few points cannot be used. It was decided to use the minimum/maximum, together with a method to deal with outliers.

Meaning of the term “climate envelope”

The above consideration of what measure of distribution to use for modelling climate envelopes, begs the question – what do we really mean by a climate envelope? This is important to define from the outset. Applying the present methods within a GIS system, for each best fit climate surface to the points, we are able to plot climate surfaces between the minimum and maximum contours bounded by those points.

The same goes for the next best surface, and in the GIS we then can look for the overlap between two surfaces using a logical “and” statement to find all cells for which that statement is true (and so on for further surfaces, minimising the total number of cells in the prediction at each step). That is how the phrase “climate envelope” should be interpreted for this report. However, traditionally, one may think of a climate envelope as an envelope around the observed points plotted around the “cloud” of observed points for 2 or even 3 climate variables. Two logical questions then arise:

1 how to delimit the points in climate space;
2 how to represent more than three variables in climate space and appropriately delimit them.

Confining the problem to thinking of a simple 2-dimensional graph, say rainfall against temperature, because we are using the minimum and maximum for each variable, in effect we are plotting a bounding box (rectangle) in climate space, whereas the points might themselves plot as (say) an elliptic shape, or even near a straight line in case of high correlation. Should we not be using a polygon though – e.g. a convex hull, rather than a bounding box, for the graphical climate description? Ideally yes, and this leads to the suspicion that the GIS solution may not be implementing a sufficiently tight-fitting climate envelope for the points (it is certainly not using a graphical method).

A full answer to this question goes somewhat outside the scope of this analysis for now, including the complication of how technically to map the more complex envelope back into geographic space, and how to describe the best envelope for the points statistically as well as mathematically. However, the illumination can be offered for now that if we search for random points within the bounding box of climate space, many if not most of those outside falling outside (let us say) a convex hull in the bounding box are not geographically feasible (i.e. not autocorrelated at all) in any case.

Therefore, there is some reassurance that in fact the GIS system is reproducing a climate model not so distant from a polygon fitted to the cloud of graphical points (at least in the 2 dimensional case) than may appear at first sight. A corollary is that (considering especially the multidimensionality problem) the graphical representation as a bounding box is misleading in that it does not display the geographic constraint envelope, when graphically feasible.

Testing the Models

Tests by removing outliers

What is the effect of outliers on the bioclimatic model predictions? Should they be kept in or thrown out? Their coordinates might anyway be correct, or at least correct to a certain spatial error.

First we must again define what we mean. By outliers in this case, we mean extreme points with respect to their distribution across all climate variables examined. Because we do not know in advance which climate surfaces to use, we can use the climate surface intersection values for each point summed across all surfaces. Those such points that were also unreliable beyond a certain spatial error, or collected by a collector who rarely labelled his specimens accurately, were rejected, for one test, and kept in, for another, in a sensitivity analysis that asked whether removing them made a significant difference (junk in may mean junk out).

Clearly, there is a potential role for expert decisions here, either in a priori validation or a posteriori manual rejection. An example of this type of test is given for the Malagasy Red Admiral which may be threatened with extinction (as has happened for its genus Antanartia in the Mascarenes). This case is illustrated in Fig. 6. Here, the same number [arbitrarily 9] surfaces were used for each case, but the point for Maromandia that lies close to sea level was taken out (in fact, there are other butterflies for which the exact same block label by Decary is suspect). Notice how much tighter the climate envelope becomes. The point-driven tournament was run independently in both cases (with and without Maromandia record), using the same surfaces.

Fig. 6. Test by removing outliers:
Malagasy Red Admiral
Antanartia hippomene madagassorum

There was also a complete change of ranking of choice among the nine surfaces, for example wb (water balance) was the first choice in fig 6A, whilst sprmat (spring mean temperature) was the first choice in fig. 6B. This result shows that the retention of outliers can have a massive effect on the results. Note also though, that there is some suitable climate space in Fig. 6B within about 30 km. (i.e. walking distance) of Maromandia (where shown in yellow), so the record may not be totally without value. One can imagine a more complex algorithm that assumes each point is an outlier and takes each out in turn to examine each point’s contribution to any overprediction.

Tests by random subpartitioning

Next, we examine another approach to building confidence in the bioclimate models. This is a subsampling approach which may be called random subpartitioning. An example is given for a rare (in the field) but fairly widespread species, Heteropsis vola. First, half of the points are randomly drawn (ideally one would first cut the points in the database down to the resolution of the grid, but for the sake of illustration, this was not necessary).

Heteropsis vola

This random set of half the points are plotted, as well as the residue. The climate envelopes picked objectively by the point-driven tournament for each point set, together with their observed climatic ranges and the cumulative geographic area of the envelope (% of Madagascar) for the 1-4 surfaces in the prediction, are shown below. Note the similarity in the bioclimatic model in each case. This approach lends a certain amount of confidence to the prediction, and the approach might be added to a future algorithm.

Fig. 7. Tests by random subpartitioning
7A. Random draw – half 60 points.
7B. Residue from random draw.

Tests by taxonomy

We might expect on average closely related species to be more similar in their response to climate than distantly related taxa. For example, the genus Strabena is known to show high levels of elevational endemism and is more diverse in montane areas (Kremen, 1992), and so we might expect that temperature might limit its distribution more than moisture.

Here for example, two species in this genus are compared that show little overlap in actual points for the climate envelopes selected by a point-driven tournament for three grids. Notice that the same first best fit variable is picked for each species, and that the second and third variables in combination are of a similar nature (rainfall, aridity) although their order is transposed. Does then a pattern emerge then across all taxa examined, or is the climate envelope picked by the algorithm unrelated to taxonomy?

Fig. 8. Example of test by taxonomy

Tests by taxonomy and habitat using single best climate model

The above question was examined both for taxonomy and habitat, for the single best fitting climate surface, by tabulating the number of “hits” for each of the 20 surfaces by genus. The total number of hits for each surface across all 48 species was ranked according to moisture-related variables and temperature related variables . Several points emerged from this analysis:

1. Moisture-related variables predominate slightly as best fits across the 48 species (29:19).
2. Some taxa (e.g. Perrotia) show largely moisture-related best fits, whilst others, such as Strabena and Heteropsis, may include groups adapted to either.
3. A pattern according to habitat preference is not obvious.
4. Particular variables score highly (e.g. water balance, spring mean temperature), especially for shorter periods of months. Surprising, estimated evapotranspiration (etp) only scores once.
5. Five variables score as best hits for no taxa, and these include mean annual rainfall and mean annual temperature. This is reassuring as well as expected because at certain times of year extremes in climate should show more discrimination as potential boundaries for point distributions.

Ideally, we would use a phylogeny (e.g. Torres et al. 2001) to test in a more robust fashion for climate response effects in taxa using the method of independent contrasts (CAIC program). Nevertheless, the analysis suggests new hypotheses (e.g.differential effects on latitudinal or elevational distribution, of dependence on moisture or temperature-related variables) that might explain the observation of Lees (2000) that lowland narrow-ranging endemics for some animal groups are concentrated nearer the equator than are upland narrow-ranging endemics (perhaps more limited by moisture).

Having explored some methods of testing the methodology for producing bioclimate envelopes (which may well benefit from further fine tuning in future), and found some potential pitfalls as well as built up some confidence in the modelling approach, the models for each species are generated. These models for the 48 species, for a 4 grid analysis in each case are shown in the Atlas. The technique is fully comparable and objective across species. In order to indicate successive layers of confidence in the model, the 1; 1 and 2; 1, 2 and 3; and 1-4 grid models are overlaid in these maps using successive shades of colour (pink, blue, light green, and dark green, with primary forest habitat added as another layer). The probability of occurrence in relation to each combination of climate grids is not calculated here (it would depend on levels of sampling intensity in sampled gridcells), but it is easy to see that for the pink area, the conditional probability of occurrence is lowest).

Atlas of predictive bioclimate models for rare Madagascan butterflies.

Application: climate change

What about application of the models? An obvious application is in predicting future climate. There was not the scope in this work to explore in any detail future climate scenarios for Madagascar, not least because none were available at the time of the work at a suitable resolution. The beauty of this system, though, is once the models have been built for each species, and the future climate scenarios are available, the best fit equations can simply be dropped into the future models under the assumption of no adaptation.

Because we lacked suitable resolution future climate surfaces, the example below merely provides an illustration of the predictive potential. Simplistic assumptions underlie this example: that for example a 1.4 deg. C global increase predicted by the Hadley CM3 model for 2050 might be sustained in the Madagascar region and evenly across the island (possibly in fact, the island’s position in the Indian ocean might buffer the effects of climate change in some way), and that it would also apply to a quarterly period. By way of an illustration, we examine what would happen to the potential climate envelope of the upland endemic Strabena mandraka based on the observed best fit envelope for all its known points. We also have to ignore for now the important question of how feasible (or necessary) it would be for the metapopulation to track its climate envelope. We cannot predict extinction without considering these questions.

Fig. 9 shows the effect of a 1.4 deg.C increase in the climate surface for spring mean temperature for 2050, assuming that the model in Fig. 9A is for 1950, and allowing for a 0.5 deg. C increase since then typical of surface air temperatures in the southern hemisphere (Jones & Briffa, 1992; Jones et al. 2001), plus a conservative prediction of + 0.9oC for 2050 (approximating to the Hadley CM3 models with the inclusion of damping sulphate aerosols, as well as CO2).

Fig. 9. One climate scenario for 2050.
This climate scenario shows a +1.4deg. C increase in surface air temperature (spring mean temperature) for 2050. Note the complete loss of some climate space. Temperature scale multiplied by 1000 [device to retain floating point in integer grids]

Strabena mandraka

Fig. 10 shows the effect of this temperature increase, uniformly applied across Madagascar, on the single surface model for Strabena mandraka. For this species, the mid-twentieth century climate envelope has an occupancy of 1493 gridcells, whilst those gridcells actually forested in 1999 numbered 347 gridcells, both out of 20343 for all Madagascar.

Fig. 10. Single surface climate envelope scenario for 1950 (left) and 2050 (right) for the butterfly Strabena mandraka
Climate envelope in green, 1999 primary forest in dark green (climatically suitable) or grey (climatically unsuitable), and all point observations in black.

The 2050 scenario gives a climate envelope occupancy just larger than the predicted current forest occupancy of 378 gridcells, but the forest occupancy (assuming no further habitat removal) declines to just 29 0.05-degree gridcells. This one test suggests that for such a montane forest species, the historical removal of high plateau forest on Madagascar in the past 2-3000 years removed a massive potential climate buffer for temperature-constrained organisms that existed during previous warming periods, whereas the principal last climate and forest refuge with sufficient elevational flexibility for this species under such a future climate scenario would be in the northern mountain forests around Tsaratanana (where however, it has never been recorded), with few climatically suitable fragments in central parts and at the western edge of the eastern rainforest chain.

This apparently grim picture needs of course to be rigorously tested on many other organisms and on more realistic climate scenarios. It is likely that for moisture-constrained organisms, the predicted response will be very different, a very interesting area for investigation.

Discussion, prospects and challenges

There are already many methods and programs for predictive climate modelling of species distributions, but most incorporate the simple idea of mapping climatic envelopes in geographic space. The idea that climate limits constrain or even drive species boundaries is almost becoming a paradigm rather than an explicitly tested hypothesis, and there is obviously much more work to be done. What is the role of geology, soil type, coastal boundaries, foodplants, hostplant specificity and tracking, metapopulation dynamics, and historical events, in species distribution, for example, and what are the relative importance of such factors to climate?

This work has developed a simple approach to climate modelling which is readily applicable (across a spectrum of species range sizes) to point data; increasingly databases are of this form. However, it should not be considered in exclusion to other approaches. Statistical methods such as autologistic regression (Augustin et al., 1996) and generalised additive models (Hastie and Tibshirani, 1990) have obvious power, whilst genetic algorithms have been shown to have great heuristic processing value in this field (Peterson et al. 2002). As shown, sensitivity analyses and randomisation procedures also have great value in evaluating model robustness. Rather, there will be great advantages in comparing or combining methods. For example, whilst observed extremes may delimit the climatic contours that are consistent with the distribution of other points, could information from all the observed points be used more powerfully? Also, whilst one extreme of a variable may be constraining within geographically feasible values, the other may not; here the use of a genetic algorithm should be helpful. With a multiplicity of possible explanations, methods should address which are the best or most parsimonious in terms of the data; maximum likelihood and bayesian methods could also very effectively be applied in this field. Some of these methods are potentially useful in generating probability surfaces. That is an area where more development of is needed for Madagascar models, considering the low number of sites recorded for absence (methods include generating response surfaces based on presence and absence: Hill et al. 1999; and programs such as CLIMEX: Sutherst and Maywald, 1999).

Whilst more methodological development is needed, much of this will be done and already is being done elsewhere. There is an urgent need now to examine climatic potential ranges of species in terms of conservation planning, and Madagascar is an ideal test case, for the many groups where good point data now exists and which may be responsive to climate change. There is also a need to use the best possible datasets for past, present and future climate scenarios and forest change analyses, on which models can be applied, and to test and compare these datasets. For example, the Hadley models are presently being remodelled at a fine resolution by the University of Cape Town (G. Midgely, pers. comm.), and such models may soon be available to which observed species climate envelopes can be fitted. Conservation International has just prepared a forest change map for 1990-2000 (L. Hannah, pers. comm.) . The advantages of predicting the future using the past and present are obvious, and if potential geographic range shifts even in a half century for some species are likely to be massive assuming species could track their envelopes, there is a need to reexamine historical data (e.g. using new ground-truthing) to see if shifts predicted for documented climate shifts have already started to happen .

How should conservationists respond to such predictions? It may be that for most organisms there is inbuilt resilience and that the lag times to local population extirpation are great. Nevertheless, the need for ecosystem restoration in a large and tremendously biodiverse island that has considerable remaining latitudinal and ecosystem connectivity has never been greater, and yet only in the last year or two are organisations seriously considering funding that kind of activity on the island (e.g. WWF’s Forest Landscape Restoration program). The old patchwork reserve model for conservation, whilst arguably very successful in arresting deforestation in many localities, will become increasingly less dependable as climate modelling leads to insights as to where most turnover in species, or most local extirpations and colonisations, are likely to occur. Only ecosystem restoration, on a large scale and in an appropriate social context, can avert a major or incipient extinction event on the island in the coming century. As the present political crisis comes to resolution, there are likely to be major opportunities to build climate modelling into forest conservation and restoration planning, and resources have to be used as efficiently as possible, taking account of natural processes and models for them. There is a chance now to plan for at very least the average-case scenario.


Guy Midgely kindly supplied the CRES climate data; Dinah Millar (also Univ. Cape Town) provided helpful advice; Si Tokumine kindly helped with the initial ArcView setup; Colin McClean generously provided and modified a genetic algorithm. Jon Lovett facilitated this work and international collaboration. Funding for this preliminary study on butterflies was provided by Conservation International through Lee Hannah. Space and resources were provided at the Natural History Museum (Biogeography and Conservation Lab.). Some of the butterfly data was compiled or supplied by Claire Kremen (Univ. Princeton) and Tiana Raharitsimba.


Augustin, N.H., Mugglestone M.A. and Buckland, S.T. 1996. An autologistic model for the spatial distribution of wildlife. Journal of Applied Ecology 33: 339-347

Bradt, H. 2002. Loggers exploit disputed vote. BBC Wildlife 20: 25.

Busby, J.R. 1991: BIOCLIM – A Bioclimatic Analysis and Prediction System. In: Margules, C.R. & Austin M.P. (eds.)

Nature Conservation: Cost Effective Biological Surveys and Data Analysis. Canberra: CSIRO: 64-68.

Carpenter, G., Gillison, A.N. and Winter, J. 1993. DOMAIN: a flexible modelling procedure for mapping potential distributions of plants and animals. Biodiversity and Conservation 2:667-680.

Coope, G.R. 1995. Insect faunas in ice age environments : why so little extinction? Pp. 55-74 In Lawton, J.H. and May, R.M. Extinction Rates. Oxford : Oxford University Press.

Cornet, A., 1974. Essai de cartographie bioclimatique à Madagascar. Notice explicative Nº 55, ORSTOM, Paris.

Eeley, H.A.C., Lawes, M. J. and Piper, S.E. 1999. The influence of climate change on the distribution of indigenous forest in KwaZulu-Natal, South Africa. Journal of Biogeography 26: 595-617.

Hastie, T. and Tibshirani, R. 1990. Generalized Additive Models. Chapman and Hall, London.

Hill, J.K., Thomas, C.D., and Huntley, B. 1999. Climate and habitat availability determine 20th century changes in a butterfly’s range margin. Proceedings of the Royal Society of London B 266: 1197-1206

Hughes, L. 2000. Biological consequences of climate change: is the signal already apparent? TREE 15: 56-61.

Jones, P.D., and K.R. Briffa. 1992. Global surface air temperature variations during the twentieth century: Part 1, spatial, temporal and seasonal details. The Holocene 2:165-179.

Jones, P.D., D.E. Parker, T.J. Osborn, and K.R. Briffa. 2001. Global and hemispheric temperature anomalies–land and marine instrumental records. In Trends: A Compendium of Data on Global Change. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tenn., U.S.A.

Jones, P.G., and Gladkov, A. 1999. FloraMap: a computer tool for the distribution of plants and other organisms in the Wild. CIAT, Cali, Colombia.

Kremen C. 1992. Assessing the indicator properties of species assemblages for natural areas monitoring. Ecological Applications 2: 203-217

Lawton, J.H. and May, R.M. Extinction Rates. Oxford : Oxford University Press.

Lees, D.C. 2000. Latitudinal distribution of elevationally endemic fauna in Madagascar and conservation planning. Pp. 273-284 In Lourenco, W.(ed.). Diversité et Endémisme à Madagascar. Mémoires de la Société de Biogéographie, Paris.

Lees, D.C., Kremen, C. & Andriampianina, L. 1999. A null model for species richness gradients: bounded range overlap of butterflies and other rainforest animals in Madagascar. Biological Journal of the Linnean Society 67: 529-584.

Lees, D.C., Kremen, C. and Raharitsimba, H. Classification, diversity, and endemism of the butterflies (Papilionoidea and Hesperioidea) of Madagascar: a revised species checklist. In: The Natural History of Madagascar (S.M. Goodman and J. Benstead, eds.). University of Chicago Press (in press).
Magurran, A.E. 1988. Ecological diversity and its measurements. Croom Helm, London.

Mayaux P., Valéry G. and Etienne, B. 2000. Mapping the Forest-Cover of Madagascar with SPOT 4-VEGETATION data. Int. J. Remote Sensing 21 (16): 3139-3144;
see also

Lovett, J. and McClean, C. 2001. Predicting the effect of climate change on Madagascan biota. Phase I (2001): Preliminary present and future modelling of selected butterfly ranges using validated climate surfaces. Phase II (2002-): Climate modelling for larger datasets. 15 pp. Unpublished proposal to CI.

Myers N., Mittermeier, R.A., Mittermeier, C.G., da Fonseca, G.A.B., and Kent, J. 2000. Biodiversity hotspots for conservation priorities. Nature 403: 853-858

Nix, H. 1986. A biogeographic analysis of Australian elapid snakes. Atlas of elapid snakes of Australia. Australian Flora and Fauna Series, Number 7. pp. 4-15. Bureau of Flora and Fauna. Canberra: Australian Government Publishing Service.

Oldeman, L.R. 1990. An agroclimatic characterisation of Madagascar. ISRIC Technical Paper 21. Wageningen: International Soil Reference and Information Centre
Parmesan, C. 1996. Climate and species’ range. Nature 382, 765-766.

Parmesan, C. Ryrholm, N. Stefanescu, C. Hill, J.K., Thomas, C.D., Descimon, H., Huntley, B., Kaila, L., Kullberg, J., Tammaru, T., Tennent, W.J., Thomas, J. and Warren, M. 1999. Poleward shifts in geographical ranges of butterfly species associated with regional warming. Nature.399: 579-583.

Peterson, A.T., Ortega-Huerta, M.A., Bartley, J., Sánchez-Cordero, V, Soberón, J., Buddemeier, R.H. and Stockwell, D.R.B. 2002. Future projections for Mexican faunas under global climate change scenarios. Nature 416: 626-629.

Raharitsimba, N.H. 1997. Alpha taxonomy of Madagascar endemic butterfly [genus] Strabena (Nymphalidae, Satyrinae) and its implication in conservation. MSc thesis, University of Kent, Canterbury.

Skov, F. 2001. Potential plant distribution mapping based on climatic similarity. Taxon 49: 503-515.

Skov, F. & Borchsenius, F. 1997. Predicting plant species distribution patterns using simple climatic parameters: a case study of Ecuadorian palms. Ecography 20: 347-355.

Stockwell, D.R.B. and Noble, I.R. 1992. Induction of sets of rules from animal distribution data: a robust and informative method of data analysis. Mathematics and Computers in Simulation 33: 385-390.

Sutherst, R.W. and Maywald, G.F. 1999. CLIMEX: predicting the effects of climate on plants and animals. Program and manual. CSIRO, Melbourne.

Tobler, W. R. 1979. Smooth pycnophylactic interpolation of geographic regions. Journal of the American Statistical Association. 74: 519-530.

Torres, E., Lees, D.C., Vane-Wright, R.I., Kremen, C., Leonard, J.A. and Wayne, R.K. 2001. Examining monophyly of a large radiation of Madagascan butterflies (Lepidoptera: Satyrinae: Mycalesina) based on mitochondrial DNA sequence data. Molecular Phylogenetics and Evolution 20 (3): 460-473.

Last Updated: 18 December 2018