Ecology

Dealing with Pseudo-Replication in Linear Mixed Models

Pseudo-replication is associated with the incorrect identification of the number of samples or replicates in a study. This can occur at the planning or execution stage of a study, or during the statistical analyses. Almost always this incorrect identification will make the statistical analyses invalid. To illustrate this, we will consider an experiment that compares two treatments applied to eight plots, each containing four plants (see figure below). In this experiment one treatment (_e.g._, control) is applied to four plots at random, and the other treatment (_e.g._, fertilizer) to the remaining four plots.  What is the replication of a given treatment in this experiment? The replication per treatment is four, but in total 16 plants are treated. The experimental units (EU) are the plots, while the individual plants, or measurement units (MU), are the pseudo-replicates. Even though we might have a total of 16 observations per treatment, these have a structure, and the treatment was randomized to the plots not to the plants. In statistical terms if pseudo-replication is not properly accounted for the degrees of freedom (_df_) are incorrect, and generally inflated. This inflation of the _df_ leads to p-values that are lower than what they should be, and therefore, with pseudo-replication we are more likely to reject our null hypothesis. For example, in the experiment above, a t-test comparing the two treatment means should have a total of 6 _df_, but under pseudo-replication we will incorrectly consider that we have 30 _df_. In this case our critical t-value at a 5% significance level will change from 2.447 to 2.042, making our p-value and confidence interval incorrect. You can imagine, that pseudo-replication has important implications for statistical inference and scientific research, however, it is common to find this mistake in research papers. In this topic, there are two important elements to consider: 1. _How to identify pseudo-replication in a study,_ and 2. _How to correctly deal with this pseudo-replication within a statistical analysis._ We will discuss both aspects with emphasis on the use of linear models (LM) and linear mixed models (LMM). ## Identifying pseudo-replication Pseudo-replication is not always easy to recognize, particularly for experiments that contain multiple layers or strata. However, the most important aspect is to clearly distinguish between experimental units (EU) and the measurement units (MU). Treatment levels should always be applied randomly to different EU, as in the above example where a treatment level (_e.g._, fertilizer) was applied to a plot, and thus all four plants within the plot have the same treatment level. The difficulty comes when we take measurements (_e.g._, plant height) from each of the plants, and therefore the EU differs from the MU. You can identify similar cases, for example where we obtain several measurements from a single EU, or even several measurements from a given MU. For example, in our experiment we could have counted insects in a sample of two leaves per plant. Here, these leaves will also be MU, but in this case, they correspond to a subsampling within each plant. Careful identification of the design structure of an experiment, or equivalently of the arrangement (or layers) of the observations will allow us to identify the existence of pseudo-replication, and also it will help us to correctly specify this multi-stratum structure into a statistical model. ## Dealing with pseudo-replication Linear mixed models (LMM) are particularly useful in dealing with the multi-stratum structure of an experiment, as they will recognize these layers and identify at what stratum there is pseudo-replication or true replication. Once the design and treatment structure are recognized the correct analysis follows easily. Further details on these structures and multi-stratum Analysis of Variance can be found in Welham _et al._ (2014). For example, we will consider a dataset from a _Eucalyptus_ hybrid study. The original study evaluated clones from several intra- and inter-specific crosses between different _Eucalyptus_ species and their hybrids (Alves _et al._ 2018). Here, we will focus only on the hybrids between two species (_Eucalyptus urophylla_ and _E. resinifera_) corresponding to a total of 17 crosses or families. The study had a randomized complete block design (RCBD) with eight blocks, where propagules (_i.e._, trees) were arranged by family in plots of six plants. Hence, there is a total of 136 plots (8 blocks x 17 families) and 816 plants (136 plots each with 6 plants) in the dataset. For this experiment, we can identify the structure (or stratum) as having replicates (or blocks) that are partitioned into several plots, and finally, we have individual trees within these plots. Hence, similar to the previous example, the plots are our EU and the trees are the MU. Under this structure the treatments (_i.e._, families) were assigned at the plot level and not at the plant level. The data for our analyses can be found in the dataset called $\\texttt{EUCAL}$ and the first few lines of this data are:  | id | rep | plot | tree | family | male | female | height | dbh | volume | | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | | 1 | 1 | 181 | 1 | 18 | 14 | 3 | NA | NA | NA | | 2 | 1 | 181 | 2 | 18 | 14 | 3 | 25.4 | 16.39 | 0.2680 | | 3 | 1 | 181 | 3 | 18 | 14 | 3 | 26.5 | 17.86 | 0.3318 | | 4 | 1 | 181 | 4 | 18 | 14 | 3 | 20.6 | 17.57 | 0.2498 | | 5 | 1 | 181 | 5 | 18 | 14 | 3 | 27.8 | 18.46 | 0.3721 | | 6 | 1 | 181 | 6 | 18 | 14 | 3 | 20.5 | 11.62 | 0.1087 | | 7 | 1 | 281 | 1 | 28 | 11 | 4 | NA | NA | NA | | 8 | 1 | 281 | 2 | 28 | 11 | 4 | 26.1 | 18.18 | 0.3386 | We observe that the first six observations correspond to plot 181 from replicate 1 However, this plot has only five trees alive at the time of measurement (data for tree 1 is missing, denoted with $\\texttt{NA}$). Hence, as with many experiments, we have a few missing observations. Now, we can proceed to write the correct model that considers the full structure of this experiment: $$ y\_{ijk} = \\mu + r\_i + f\_j + p\_{ij} + \\epsilon\_{ijk} $$ where $y{ijk}$ is the phenotypic observation of the kth plant from the jth family on the ith replicate (or block); $\\mu$ is the overall mean;  $ri$ is the fixed effect of the ith replicate; $f\_j$ is the fixed effect of the jth family; $p\_{ij}$ is random effect of the ijth plot, with $p\_{ij} \\sim N(0, \\sigma^2\_p)$; and $\\epsilon\_{ijk}$ is the random residual with $\\epsilon\_{ijk} \\sim N(0, \\sigma^2\_e)$. In the above model we have defined the replicate and family terms as fixed effects. But depending on the objective of the analyses these could be considered random effects. The plot effect, $p\_{ij}$, was assumed to be random, because it is used here to describe the design structure of the study. The above model is easily fitted in much of the statistical software that will fit LMM. In this example, we will use ASReml-R as implemented in R. ## Fitting a LM to plant-level data To demonstrate the effects of pseudo-replication in this experiment, we will be evaluating a few different models. First, we fit an incorrect model: one that ignores the plot effects, and therefore has pseudo-replication. The code for this model in ASReml-R is: ```plaintext modINC&lt;-asreml(fixed=height~rep+family, data=EUCAL) ``` Recall that both $\\texttt{rep}$ and $\\texttt{family}$ are factors considered fixed effects. Our dataset called $\\texttt{EUCAL}$ has a row for each tree, and therefore a total of 816 records. Here, we are interested in the ANOVA table and the predicted means for the 17 families. These can be obtained using the commands: ```plaintext wald.asreml(modINC,denDF='numeric')$Wald predsINC&lt;-predict(modINC,classify='family')$pvals ``` resulting in the ANOVA table: |   | Df | denDF | F.inc | Pr | | --- | --- | --- | --- | --- | | (Intercept) | 1 | 686 | 20280.000 | 0.0000000 | | rep | 7 | 686 | 1.062 | 0.3864759 | | family | 16 | 686 | 2.237 | 0.0036815 | The predicted means for the first six families are: | family | predicted.value | std.error | status | | --- | --- | --- | --- | | 18 | 24.06383 | 0.6799762 | Estimable | | 28 | 25.72138 | 0.7499339 | Estimable | | 102 | 23.48852 | 0.7052685 | Estimable | | 104 | 22.91479 | 0.7311813 | Estimable | | 158 | 26.25878 | 0.6799763 | Estimable | | 172 | 24.60833 | 0.6582289 | Estimable | In our analysis, we obtained an estimate of the residual variance, $\\hat{\\sigma}^2\_e$ of 20.80 (not shown), and the ANOVA table has 686 denominator _df_ for each fixed term. To better understand this ANOVA table we can write it as: | Source | df |   | | --- | --- | --- | | Replicate | r – 1 | \= 7 | | Family | f – 1 | \= 16 | | Error | r x f x m – r – f + 1 | \= 792 | | Total | r x f x m -1 | \= 815 | where r = 8 is the number of replicates, f = 17 is the number of families, and m = 6 is the number of plants per plot. The significance of replicate and family are assessed using the F-values of $F\_{7,792}$ and $F\_{16,792}$, respectively. Therefore, we have 792 _df_ for the denominator for these tests for any confidence interval we might calculate and for any contrast or multiple comparison we might desire to construct from this analysis. The reason why our ANOVA table from ASReml-R has only 686 _df_ is because of the missing values in the dataset. For the response height, there are 106 missing values: 792 – 106 = 686. The ANOVA denominator _df_, 686, is inconsistent with the total number of plots in the experiment: f x r = 136, a much smaller number. ## Fitting a LM to plot-level data Probably one of the easiest ways to deal with multi-stratum data that avoids complex models, but that still fits a correct model, is to collapse the data by obtaining means, and then use these means as our response variable instead of the individual observations. The dataset called $\\texttt{EG\_PLOT}$, has a total of 136 observations (_i.e._, plots), and the first few lines of this data are: | id | plot | rep | family | male | female | height | dbh | volume | | --- | --- | --- | --- | --- | --- | --- | --- | --- | | 1 | 181 | 1 | 18 | 14 | 3 | 24.1600 | 16.3800 | 0.2661 | | 2 | 281 | 1 | 28 | 11 | 4 | 24.1600 | 15.1160 | 0.2388 | | 3 | 1021 | 1 | 102 | 13 | 9 | 22.8167 | 13.7450 | 0.1912 | | 4 | 1041 | 1 | 104 | 14 | 9 | 21.7667 | 13.2533 | 0.1785 | | 5 | 1581 | 1 | 158 | 14 | 14 | 24.7500 | 15.2267 | 0.2374 | | 6 | 1721 | 1 | 172 | 13 | 15 | 24.7667 | 16.0083 | 0.2535 | | 7 | 1931 | 1 | 193 | 14 | 16 | 21.7750 | 13.2775 | 0.1900 | | 8 | 1991 | 1 | 199 | 11 | 16 | 25.3000 | 15.9567 | 0.2631 | The first observation for height corresponds to plot 181 from replicate 1, with a value of 24.16 m. This is the arithmetic average of five trees (recall that there is one dead plant for this plot). Now, we are ready to fit the correct model to the plot means using the following code: ```plaintext modPLOT&lt;-asreml(fixed=height~rep+family, data=EG_PLOT) ``` The above model does not contain any random effects, as plot is now absorbed by the residual term. That is, the residual term represents the error associated with the plot mean for height. Before, we look into the output, a generic ANOVA table for this model is: | Source | df |   | | --- | --- | --- | | Replicate | r – 1 | \= 7 | | Family | f – 1 | \= 16 | | Error | (r – 1) x (f – 1) | \= 112 | | Total | r x f – 1 | \= 135 | Here, we have that the significance of the replicate and family factors are assessed with F-values of $F\_{7,112}$ and $F\_{16,112}$, respectively. Hence, the denominator _df_ here is 112 instead of 792 (or 686 due to missing values) from the previous model. As you can expect, this reduction in _df_ will have important consequences. We can see this from requesting from ASReml-R the ANOVA table and the predicted means: |   | Df | denDF | F.inc | Pr | | --- | --- | --- | --- | --- | | (Intercept) | 1 | 686 | 20280.000 | 0.0000000 | | rep | 7 | 686 | 1.062 | 0.3864759 | | family | 16 | 686 | 2.237 | 0.0036815 |   The predicted means for the first six families are: | family | predicted.value | std.error | status | | --- | --- | --- | --- | | 18 | 24.05501 | 0.8324633 | Estimable | | 28 | 25.77875 | 0.8324633 | Estimable | | 102 | 23.08334 | 0.8324633 | Estimable | | 104 | 22.99771 | 0.8324633 | Estimable | | 158 | 26.27583 | 0.8324633 | Estimable | | 172 | 24.60834 | 0.8324633 | Estimable | The ANOVA table from ASReml-R provides us with the correct _df_, and we can see that the significance of the family factor has a p-value of 0.0433 for this mean model analysis, in contrast with the p-value of 0.0037 for the incorrect model. If we are using a significance level of 1%, then our conclusions will be quite different! This is a clear consequence of pseudo-replication, where we have arrived at a different conclusion due to our incorrect consideration of the true-replication of an experiment. The predicted values for the means in both analyses are very similar, with differences only due to how the missing data were treated (more on this below). But, the standard errors of the mean (SEM) from this analysis are 0.83, a value larger than the ones obtained earlier that ranged from 0.66 to 0.81; again, this is due to the fact that we have a different number of replications considered between the analyses. ## Fitting a LMM to plant-level data Having identified the differences between an incorrect model that uses the plant-level data against a correct model that uses the plot-level data, now we are in a position to fit the correct model to the plant-level data based on our generic model defined before: $$ y\_{ijk} = \\mu + r\_i + f\_j + p\_{ij} + \\epsilon\_{ijk} $$ and this can be fitted with the code: ```plaintext modCC&lt;-asreml(fixed=height~rep+family, random=~rep:plot, data=EUCAL) ``` Here, we are specifying the plot term to be random, but this is written as $\\texttt{rep:plot} to ensure that each plot is uniquely identified (for example, we want to differentiate between plot 1 – replicate 1 from plot 1 – replicate 2). As before, we will have a look at the generic ANOVA table for this model: | Source | df |   | | --- | --- | --- | | Replicate | r – 1 | \= 7 | | Family | f – 1 | \= 16 | | Plot | (r – 1) x (f – 1) | \= 112 | | Error | r x f x (m – 1) | \= 680 | | Total | r x f x m – 1 | \= 815 | The above ANOVA table looks very similar to the one from our incorrect analyses, but it incorporates an additional row for the plot factor. In this case, the significance of the replicate and family factors will be assessed using the F-values of $F\_{7,112}$ and $F\_{16,112}$. These are the same tests used for the plot-mean analysis, and they are using the correct _df_: 112 instead of 792. The results from fitting our LMM to this dataset using our full model produces the following results: |   | Df | denDF | F.inc | Pr | | --- | --- | --- | --- | --- | | (Intercept) | 1 | 107.5 | 1.537e+04 | 0.0000000 | | rep | 7 | 107.8 | 8.001e-01 | 0.5889034 | | family | 16 | 109.2 | 1.714e+00 | 0.0541757 | \<br/> | family | predicted.value | std.error | status | | --- | --- | --- | --- | | 18 | 24.06071 | 0.7860554 | Estimable | | 28 | 25.73520 | 0.8441194 | Estimable | | 102 | 23.43381 | 0.8127377 | Estimable | | 104 | 22.93619 | 0.8313801 | Estimable | | 158 | 26.26424 | 0.7860553 | Estimable | | 172 | 24.60833 | 0.7678512 | Estimable | The first thing to note is that the denominator _df_ for replicate and family is ~108 instead of 112. This value is approximated using the Kenward-Roger correction for _df_, which aims to help with testing the significance of a model factor. This correction considers that under LMM often the null distribution of the test is unknown, and it accounts (penalizes) for the additional uncertainty introduced by the fact that we are using estimates of the variance component and not its true values. For this reason, the denominator _df_ approximated under Kenward-Roger correction tends to be lower than the ones expected from our ANOVA table. For this analysis, based on a significance level of 1%, we have a p-value for the family term of 0.054, a value that will lead us to similar conclusions as the one from the plot-level data (p-value of 0.043). However, the p-value here is slightly larger because the Kenward-Roger correction has been used. The above output resulted in predicted values for each family that are similar to the other analyses; however, the SEMs are much closer to the ones from the plot-level analysis. Again, differences are mostly due to some dead trees (_i.e._, missing data). One additional piece of information that we can obtain from our correct LMM are the estimates of two variance components. These are obtained with the following command: ```plaintext summary(modCC)$varcomp ``` resulting in the output: |   | component | std.error | z.ratio | bound | %ch | | --- | --- | --- | --- | --- | --- | | rep:plot | 1.45570 | 0.731542 | 1.989907 | P | 0 | | units!R | 19.56639 | 1.152809 | 16.972794 | P | 0 |   The plot variance estimate $\\hat{\\sigma}^2\_p$, which corresponds to the between-plot variance, has an estimated value of 1.456. For this dataset, the estimated between-plot variance is much smaller than the residual or within-plot variance estimate, $\\hat{\\sigma}^2\_e$, with a value of 19.57. These estimates tell us that this experiment has considerably more variability between trees in a given plot (approximately 93%), than between plots. This information can be useful to improve the design of future experiments. For example, given that in this case most of the variability is at the tree-within-plot level, it might be recommended to have larger (_i.e._, more than six trees) but fewer plots than used at the present. But this might require a reduction in the number of EU if we want to keep the total number of trees constant. It is also interesting to note that the plot-level residual variance corresponded to 5.54, a smaller value than the one obtained from our LMM, with a value of 19.57. This difference occurs because the residual variances represent different units: the former is from plot-means, and the latter is from individual plants. It is possible to connect them using the expression: $$ \\sigma^2\_{e.plot} = \\sigma^2\_p + \\sigma^2\_{e.plant}/m $$ where $\\sigma^2\_{e.plot}$ and $\\sigma^2\_{e.plant}$ are the residual variances for the plot-level and plant-level analyses, respectively; and $\\sigma^2\_p$ is the plot variance from the plant-level analysis. Finally, mm is the number of MU per EU. In our experiment this $\\texttt{m}$ will correspond to 6 trees/plot, but given that we have some missing values we must estimate $\\texttt{m}$ using the harmonic mean on the number of non-missing records per plot. This corresponds to $\\hat{m}$ = 4.89. Therefore, we obtain: $$ \\hat{\\sigma}^2\_p + \\hat{\\sigma}^2\_{e.plant}/m = 1.46 + 19.57/4.89 = 5.46 $$ a value quite close to the original estimate of $\\hat{\\sigma}^2\_{e.plot}$ = 5.54. It is possible to assess different experimental designs with varying numbers of plants per plot mm, with the use of the above equation, followed by power calculations of the statistical models. This is an additional advantage of LMM: they not only provide the correct analysis for any multi-stratum structure, but they also provide estimates of the variance components for each stratum, allowing us to understand the partition of variability of the response variable of interest. ## Closing Remarks We have explored the consequences of pseudo-replication in some statistical models, in that they can lead to incorrect conclusions. In addition, we have fitted correct models using plot-level or plant-level data, which both led us to the expected ANOVA tables and _df_. But, which of these two is a better approach? Clearly the answer is the LMM as it uses all data at its appropriate level and also it provides us with an estimation of the contribution of the variability for each stratum. One additional advantage, which was not discussed above, is that the LMM allows us to treat missing data at random easily. In our example, there were a different number of live plants per plot. This was not an issue for our LMM methodology as it used all raw information with its appropriate weight to obtain unbiased family mean estimates. Finally, note that our approach to deal with pseudo-replication, and indirectly with multi-stratum structure, using LMM can be easily extended to other more complex cases. For example, those where we have sub-sampling within measurement units. In all cases, using random effects that identify each stratum will take care of the complex nature of our experiment providing us with the correct degrees of freedom and therefore, our main goal of obtaining the correct statistical inference is achieved. ### **Author** Salvador A. Gezan ### **References** Alves, RS; de Carvalho Rocha, JRdAS; Teodoro, PE; de Resende, MDV; Henriques, EP; Silva, LA. 2018. Multiple-trait BLUP: a suitable strategy for genetic selection of _Eucalyptus_. Tree Genetics & Genomes 14, 77 (2018). Welham, SJ; Gezan, SA; Clark, SJ; Mead, A. 2014. Statistical Methods in Biology: Design and Analysis of Experiments and Regression. Chapman and Hall/CTC. **Notes**: SAG April-2020.

Food Standards Agency: safer food for the nation

Everything we do reflects our vision of ‘Safer food for the nation’. We aim to ensure that food produced or sold in the UK is safe to eat, consumers have the information they need to make informed choices about where and what they eat and that regulation and enforcement is risk-based and focused on improving public health.  Genstat has been deployed in many important areas of FSA activity, including the following:- 1\. Exposure Assessment There are many types of chemical whose presence in food can pose some level of risk for consumers. These include environmental contaminants (e.g. lead or cadmium), as well as certain pesticides and food additives. The FSA needs to estimate the likely level of human exposure to such chemicals, either from individual foods, or from the diet as a whole. Our exposure assessment models combine occurrence data (on the level of a given contaminant in foods) with extensive human intake data (from individual food diaries). Genstat provides excellent tools for the necessary data import, indexing, manipulation and analysis and produces fast, reliable results. It has the flexibility to perform more complex modelling when required: e.g.  probabilistic assessments of risk that allow for variability and uncertainty in some of the key input variables. 2\. Risk Assessment The FSA uses Hierarchical Generalized Linear Models to model the probability of elevated levels of naturally occurring biotoxins in shellfish. These estimates are used to prioritise biotoxin sampling activity around the shores of the UK, indicating where, and at what times of year, the level of sampling should be raised, and also where and when it can be reduced. Genstat has excellent procedures for performing models of this kind and for extracting model coefficients in a form that can be used to make a range of policy-relevant inferences.

MIKE: tracking roaming behaviour in elephants

It’s not news that elephants are endangered; for literally thousands of years ivory has driven hunters to kill these magnificent beasts. Combine this deliberate destruction with increased land use change and habitat transformation by private enterprise, international aid and government expansion across Africa and these magnificent beasts are still under threat even within the national parks and other protected areas set aside for them.  However, before we get too carried away with horrific idea that elephants may well become extinct; there is another side to this. Two separate taxa, or races, of elephants exist in Africa today: the savannah elephant, which is well known and has been studied and valued by the international tourist trade for decades, and the diminutive forest elephant, which until recently has remained in obscurity beneath the trees. In the Congo Basin, an area noted for being one of the great wilderness areas with the second largest rainforest on earth – it is home to an outstanding range of biological diversity; forest elephants have until recent decades thrived. These types of wildernesses have always been recognised for their high conservation value, being home to a wide variety of flora and fauna, and as such provide invaluable resources of timber, minerals and other natural resources. The Congo Basin is a typical case where the resources available are in high demand in rich, industrial nations and are also critical to the economies of poor nations. And it is at this point we in the Western developed world need to be very careful about conservation and conservation techniques; the international conservation efforts are aimed at managing and promoting sustainable natural resource exploitation and functional ecosystems. It is incredibly important that we recognise the needs of people in the developing world, and try to ensure they are able to exploit their natural resources in ways that do not fundamentally alter the ecology of the environment, upon which millions of rural people still depend. The high value of resources in these wilderness areas means that it is only fair that Africans have an opportunity to exploit and benefit from them and are encouraged to do so legally, ethically and with an eye on their own local ecology. Complications occur as industrial expansion pushes into the wildernesses to obtain these natural resources, an expansion that relies on a network of roads and rail infrastructure. It is this increase in roads that appears to be causing some of the biggest threats to the forest elephants. Illegal elephant killings tend to occur close to roads, since these areas are easily accessible, so the density of forest elephants in the protected areas is determined, not by the size of the area but by how many roads criss-cross through the wilderness. Elephants are by their nature, shy and avoid human contact as much as possible, so how are they reacting to the increasing number of roads across their home ranges? Will they adopt a “siege” mentality and reduce their home range to avoid roads, or a “skirmish” strategy, where they will continue to range widely, irrespective of the number of roads? I spoke with Stephen Blake of the Wildlife Conservation Society and Max Planck Institute of Ornithology who has recently headed up a study into this very issue. His team fitted a large sample of forest elephants in the Congo Basin with GPS telemetry collars, so that they could track their roaming behaviour. They collected data from 28 elephants located in conservation areas in the Congo, Central African Republic and Gabon. Figures 1 & 2 below show data collected from 28 elephants located in the Congo, Central African Republic and Gabon using Genstat. Figure 1: Relationship between area of MCP home range and the area of road-less polygon; y axis – Home range size (kmsq), x axis roadless polygon size (kmsq).  Figure 2: Maximum linear displacement versus area of road-less polygon; y axis maximum linear distance moved (km), x axis roadless polygon size (kmsq).  “The data we collected covered the sex of the elephant, the site (the protected area in which the elephants were collared), the number of roads in the area, the number of days the elephants were collared,” explained Stephen. “I initially analysed the result collected in Genstat, using linear regression for data exploration followed by fitting standard curves to the data (Figures 1 and 2). Genstat was vital in the research programme as I can trust the statistical procedures. I was able to talk directly with developers when I came up against any queries. I have found the Genstat support team to be extremely helpful, I usually get a response to a query within 24 hours and often much sooner. We could have gone much further with the analysis, but in order to understand the influence of roads on forest elephant ranging, the regression options in the Genstat menu showed us clear statistical patterns.” Results indicated that the forest elephants routinely crossed protected roads, however only one elephant crossed an unprotected road, and of these one did so, as far away from the nearest village as she could, and her speed increased as she crossed the road. “Our study showed that elephants seem to adopt a siege mentality with the encroaching roads”, said Stephen Blake, “and the issue with this is that over time the elephants will have less and less access to food, in quantity and variety, leading to food competition and aggression; both of which are likely to lead to a reduction in population.” “The main issue seems to be that elephants are risk averse, hence why they avoid roads where they can. So here we have a problem; central African states need the income earned from the natural resources found in their countries, but not only do we not want to see the elephants themselves disappear, a change like that to any ecosystem is likely to have more wide ranging impacts.” “There are options for policies to overcome the issues, before it is too late,” says Stephen. “Based on our research, we recommend that no new roads are built on existing forest elephant conservation areas, but to access the natural resources near these sites temporary roads should be built, so that once the resources are exhausted the area can easily return to natural wilderness. The other option is to try to reduce the elephants’ fear of roads, by preventing settlements on the roads, enforcing anti-poaching. These are of course expensive options, but without addressing this we face losing these beautiful creatures.”