Conservation by Consensus: Reducing Uncertainty from Methodological Choices in Conservation-Based Models
by
Mark S. Poos
A thesis submitted in conformity with the requirements for the degree of Doctorate of Philosophy Department of Ecology and Evolutionary Biology University of Toronto
© Copyright by Mark S. Poos 2010
Conservation by Consensus: Reducing Uncertainty from Methodological Choices in Conservation-Based Models Mark S. Poos Doctorate of Philosophy Department of Ecology and Evolutionary Biology University of Toronto 2010
Abstract Modeling species of conservation concern, such as those that are rare, declining, or have a conservation designation (e.g. endangered or threatened), remains an activity filled with uncertainty. Species that are of conservation concern often are found infrequently, in small sample sizes and spatially fragmented distributions, thereby making accurate enumeration difficult and traditional statistical approaches often invalid. For example, there are numerous debates in the ecological literature regarding methodological choices in conservation-based models, such as how to measure functional traits to account for ecosystem function, the impact of including rare species in biological assessments and whether species-specific dispersal can be measured using distance based functions. This thesis attempts to address issues in methodological choices in conservation-based models in two ways. In the first section of the thesis, the impacts of methodological choices on conservation-based models are examined across a broad selection of available approaches, from: measuring functional diversity; to conducting bio-assessments in community ecology; to assessing dispersal in metapopulation analyses. It is the goal of this section to establish the potential for methodological choices to impact conservation-based models, regardless of the scale, study-system or species involved. In the second section of this thesis, the use of consensus methods is developed as a potential tool for reducing uncertainty with methodological choices in conservation-based models. Two separate applications of consensus methods are highlighted, including how consensus methods can reduce uncertainty from choosing a modeling type or to identify when methodological choices may be a problem. ii
Acknowledgments No accomplishment is ever singular, and my doctoral work is no different. For me, I have had the great fortune of a loving family, a group of wonderful friends, supportive colleagues, a thorough and considerate academic committee, and a brilliant advisor. I think Marston Bates was right when she said “Research is the process of going up alleys to see if they are blind.” In my personal journey through this doctoral work, I have gone through many blind alleys. Without the help and support of all the people in my life, this thesis would not have been possible. First, I am grateful for the opportunity to have conducted this research under the supervision of Dr. Don Jackson. I have learned a great deal about a great many things from Don, including: multivariate statistics, sampling aquatic systems and natural ecology. Don was always willing to share his thoughts and ideas, leave his door open door for questions, and let me steal a cup of coffee; you couldn`t ask for a better combination. I am also very thankful for the mentorship of Dr. Harold Harvey. Harold is my academic grandfather, a mentor and I hope to say friend. I am thankful for Harold’s advice, his stories, and providing laughter and support. I will greatly miss our daily conversations (and the recession cookies). I hope one day I am half as wise as Harold is. I am deeply indebted to Dr. Nicholas Mandrak for all his help and advice. Not only was Nick my M.Sc. supervisor, he was a member of my Ph.D. academic committee, and he was integral part from the start of my thesis. Nick’s insights into biology of fishes are second to none, and I owe a great deal of gratitude for all his mentorship throughout the years. Nick was also instrumental in helping to obtain funds (IRF #1410) to keep aspects of this project going; which saved the project. I am also greatly indebted to my academic committee members: Dr. MarieJosee Fortin, Dr. Brian Shuter and Dr. Keith Somers, whose expertise and knowledge greatly improved this thesis and my own research. All of my academic committee members were wonderful in guiding me through their areas of expertise, from spatial ecology (Thanks MarieJosee), to fisheries techniques (Thanks Brian), to multivariate statistics and bioassessments (Thanks Keith). Finally, I am thankful to Dr. Bill Matthews for acting as my external examiner. Bill has always been one of my academic idols, and it was an honor and pleasure to have him as an external. I will never forget it. From my advisor, to my academic committee, to those involved on the appraisal/defense; this thesis was clearly built on the shoulders of giants. They say it takes a village to raise a child, for this thesis, it took a small army. None of the field work would have been possible without the dedication and support of many people. In iii
particular I wish to thank Andrew Drake (aka D.A.R. Rubington III) and Cavan Harpur. So frequently I’d found myself in dire straits, without field crew or transportation, a rag tag assemblage of gear, missing GPS coordinates, lack of fish (or fill in any number of other circumstances); yet Andrew and Cavan were always there to help out and make it work. Their perseverance and friendship turned what should have been an unquestionable failure into what I hope is a great success. I can’t thank them enough. An enormous thank you goes to the Toronto Region Conservation Authority, and in particular the Aquatics group: Christine Tu, Dave Lawrie, Trevor Parker, and Tim Rance. Through a simple handshake, we became partners, and their help and support ensured that much of this research would get done. When no-one was there to help, TRCA was always ready to come to the rescue. Regardless of the obstacles, their support never faltered and I am so happy that we worked together. Dave Lawrie deserves special recognition for being the absolute best colleague you could ask for. Dave was always willing to give you the shirt off his back, I know of few people with his dedication and passion for preserving aquatic species. Too all the volunteers who helped sample and collect data; thank you. Of course, sampling endangered species comes with its own set of surprises. The fire ants were my favorite. My guess is that Dave was right when he observed a negative correlation between returning volunteers and fire ants; yet so many of you returned for multiple feedings. So thank you all the dedicated volunteers from- the University of Toronto - Andrew Drake, Cavan Harpur, Brie Edwards, Meg St. John, Maggie Neff, Paul Venturelli, Alex Manning, Monica Granados, Jonathan Ruppert, Sapna Sharma, Nicole Puckett, Cristal Hart, Chris Howard, Moe Luksenberg, John Brett, Don Jackson, Steve Walker – and the Toronto Region Conservation Authority - Dave Lawrie, Christine Tu, Trevor Parker, Tim Rance, Brennan Paul, Elyssa Elton, Brad Stephens, Cristal Hart, Laura DelGiudice, Brian Moyle, and Maria Parish – And: Peter Ng, Michell Wong Ken, Connie Zehr, Bev Edwards, Yuko Nozoe, Derek Trim, Kenny Lee, Kari Jean (ABCA), Davin Heinbuck (ABCA), Doug Forder (Ontario Streams), Daniel Morodvanschi, and Kenny Lee. I also wish to thank the Jackson lab past and present (Sapna Sharma, Steve Walker, Maggie Neff, Meg St.John, Brie Edwards, Monica Granados, Riku Pavola, Lifei Wang, Jean Bernard Caron, Angela Strecker, Karen Wilson, and Theo Willis), and Harvey labs (Andrew Drake, Cavan Harpur) for their enormous support of my project. Most (if not all) members of these labs participated in some form or another in my project, and if not, their influence from discussions can be found throughout the following pages. In particular thanks to Steve Walker for his help in iv
developing tests for measuring sensitivity in functional diversity. I never thought it would take as long or work out nearly as well as it did; my sincere thanks to Steve. I thank the following agencies and individuals for provided data for portions of this research: Royal Ontario Museum (ROM), Fisheries and Oceans (DFO), Ontario Ministry of Natural Resources (OMNR), Toronto Region Conservation Authority (TRCA), Conservation Halton (CH), Lower Lake Simcoe Conservation Authority (LLSCA), Credit Valley Conservation, Dr. Nicholas E. Mandrak (DFO), Doug Forder (Ontario Streams), Jeff J. Anderson (LLSCA), Les Stanfield (OMNR), Erling Holm (ROM), David Lawrie (TRCA), Scott Jarvie (TRCA), Sherwin Watson-Leung (CH), Dr. Scott Reid (OMNR), John Pisapio (OMNR), and Mark Heaton (OMNR). To my extended academic family, the department of Ecology and Evolutionary Biology, thank you. Dr. Lisa Manne was very kind in providing desk space at UTSc in the later stages of my thesis for writing. I am ever so grateful as it saved me many hours of commuting and kept me close to my family (who also thank you). I also really enjoyed coming in and talking ecology with her and her lab (Caroline Tucker, Christopher Grouios, and all the undergraduate students). Dr. Helen Rodd has also been a great sounding board for ideas and has been wonderfully supportive. Dr. Locke Rowe and Dr. Spencer Barrett provided laboratory space for imaging and genetics work, which didn’t make its way into the thesis, but will hopefully be out in publication shortly and I thank them nonetheless. Bill Cole and Jen Perry were also wonderful in helping me navigate through their worlds of genetics and imaging. To my friends and family ... thanks. Mom, dad, Nicole, Dave, Mau, Evan, Hayley, Shannon, my extended family, and friends; thanks for keeping me sane and always bringing a smile to my face. Last but not least, none of this would be possible without the constant love and encouragement from my wife Jessica and son Jacob. My absolute favourite part of everyday is coming home to see you. Thank you just isn’t enough. Thanks for accepting the fish smell in the car (you just can’t get that stink off), for missing your own events, or getting dropped off early/late to accommodate mine, for the missed evenings and weekends, the added stress, the crappy pay, the long hours, and everything else. To Jacob, I hope one day you can take some inspiration from these pages. Not for its content (I would never do that to you), but from reading in-between these lines and noting a worthy lesson: that even when something is difficult, perseverance and hard work can help you find your way through any blind alley.
v
Table of Contents Acknowledgments.......................................................................................................................... iii Table of Contents ........................................................................................................................... vi List of Tables .................................................................................................................................. x List of Figures ............................................................................................................................... xii List of Appendices ....................................................................................................................... xvi Chapter 1 General Introduction ...................................................................................................... 1 I) Impacts of Methodological Choices in Conservation-Based Models .................................... 3 Functional Diversity............................................................................................................ 4 Bioassessments ................................................................................................................... 5 Metapopulation Viability Models ....................................................................................... 6 II) Using Consensus Methods to Reduce Uncertainty from Methodological Choices .............. 7 Scope .......................................................................................................................................... 8 Statement of Contribution .......................................................................................................... 9 Publication of Thesis Material ................................................................................................... 9 References ................................................................................................................................ 10
Section I: ...................................................................................................................................... 18 The importance of methodological choices in ecological models ................................................ 18 Chapter 2: The importance of methodological choices in influencing the measure of functional diversity across ecological communities ..................................................................... 19 Abstract .................................................................................................................................... 19 Introduction .............................................................................................................................. 20 Methods .................................................................................................................................... 22 Results ...................................................................................................................................... 26 Discussion ................................................................................................................................ 29 vi
Acknowledgements .................................................................................................................. 32 References ................................................................................................................................ 33 Appendices ............................................................................................................................... 36 Chapter 3: Addressing the removal of rare species in bioassessments with other choices in multivariate analyses ..................................................................................................................... 43 Abstract .................................................................................................................................... 43 Introduction .............................................................................................................................. 44 Methods .................................................................................................................................... 47 Evaluating Decisions in Multivariate Bioassessments ..................................................... 48 Assessing the Statistical Argument ................................................................................... 51 Assessing the Biological Argument .................................................................................. 52 Results ...................................................................................................................................... 54 Discussion ................................................................................................................................ 58 Acknowledgements .................................................................................................................. 61 References ................................................................................................................................ 61 Appendices ............................................................................................................................... 71 Chapter 4: Contrasting direct versus indirect dispersal in metapopulation viability analyses ..... 73 Abstract .................................................................................................................................... 73 Introduction .............................................................................................................................. 74 Methods .................................................................................................................................... 77 Determining Metapopulation Viability ............................................................................. 78 Colonization Rate: ............................................................................................................ 78 Extinction Rate: ................................................................................................................ 79 Incorporating Dispersal Directly into the Metapopulation Model:................................... 80 Regional Stochasticity ...................................................................................................... 81 Comparing Viability of Metapopulations Using Direct versus Indirect Parameterization ................................................................................................... 81 vii
Results ...................................................................................................................................... 82 Metapopulation Dynamics ................................................................................................ 82 Metapopulation and Patch Viability ................................................................................. 83 Discussion ................................................................................................................................ 87 Acknowledgements .................................................................................................................. 92 References ................................................................................................................................ 92 Appendices ............................................................................................................................. 100
Section II: ................................................................................................................................... 102 Reducing uncertainty from methodological choices using consensus methods ......................... 102 Chapter 5: Reducing uncertainties in modeling the distribution of endangered species using habitat-based ensemble models................................................................................................... 103 Abstract .................................................................................................................................. 103 Introduction ............................................................................................................................ 104 Methods .................................................................................................................................. 106 Study Area and Species .................................................................................................. 106 Building Individual Models ............................................................................................ 109 Evaluating Individual Models ......................................................................................... 113 Building Ensemble Models ............................................................................................. 113 Evaluating Ensemble Models ......................................................................................... 114 Results .................................................................................................................................... 115 I) Individual Models ....................................................................................................... 115 II) Ensemble Models ....................................................................................................... 117 Discussion .............................................................................................................................. 120 Conclusion.............................................................................................................................. 123 Acknowledgements ................................................................................................................ 123 viii
References .............................................................................................................................. 124 Appendices ............................................................................................................................. 131 Chapter 6: Using consensus methods to identify (and reduce) sensitivity from methodological choices when measuring functional diversity ............................................................................. 136 Abstract .................................................................................................................................. 136 Introduction ............................................................................................................................ 137 Methods .................................................................................................................................. 139 Using Consensus Methods to Identify Uncertainty when Measuring FD ...................... 140 Results .................................................................................................................................... 141 The Relationship between FD, Distance Measure & Clustering Algorithm .................. 141 Identifying Sensitivity in FD Using Consensus Methods ............................................... 144 Discussion .............................................................................................................................. 144 Acknowledgements ................................................................................................................ 147 References .............................................................................................................................. 147 Appendices ............................................................................................................................. 152 Chapter 7: General Conclusions ................................................................................................. 153 A) Conservation-based Models in General ............................................................................ 153 B) Functional Diversity .......................................................................................................... 154 Recommendations .................................................................................................................. 157 References .............................................................................................................................. 158
ix
List of Tables Table 2.1 – The maximum probability of FD sensitivity for five communities previously used to examine FD (Petchey and Gaston 2007; Podani and Schmera 2006). The number of species and the number of assemblage pairs tested are also shown. Data sources are as follows: A) Holmes et al. (1979), B) Munoz and Ojeda (1997), C) Jaksic and Medel (1990), D) Golluscio and Sala (1993), and E) Chapin et al. (1996). ............................................................................................. 29 Table 3.1 – Summary of ordination techniques, similarity coefficients and exclusion of rarely sampled species being compared. Abbreviations are indicated in parentheses and used in subsequent figures and tables. All four approaches described in the “Exclusion of Rarely Sampled Species” were used in each of the four “Similarity Coefficient” combinations with both PCoA and NMDS. As Correspondence Analysis has the implicit chi-squared distance measure, only the four approaches used in the “Exclusion of Rarely Sampled Species” were included in that set of analyses. ....................................................................................................................... 49 Table 3.2 –Partitioning of variation in sum of squared deviations of Procrustes analyses (m2 statistic) across various choices in multivariate analyses, including: i) removal of rare species; ii) ordination technique; and, iii) choice of similarity measure. Abbreviations are those noted in Table 1. ......................................................................................................................................... 56 Table 4.1 – Summary of mark-recapture information of the endangered fish the reside dace (Clinostomus elongatus) used to directly parameterize stochastic patch occupancy metapopulation models. Shown are two locations in the Greater Toronto Area, Ontario, Canada: A) Leslie Tributary, and B) Berczy Creek. Note: items denoted with a single asterisk (*) represent a significant difference between populations (Mann-Whitney U for average dispersal, G-test for stationary tags, p<0.05). ............................................................................................... 82 Table 4.2 – Intrinsic mean time (in years) to extinction (Grimm and Wissel 2004) of two metapopulations of the endangered fish the redside dace (Clinostomus elongatus) for different levels of regional stochasticity. ..................................................................................................... 84 Table 5.1 – Summary of the hierarchical habitat-based model used for predicting the presence of endangered minnow, redside dace (Clinostomus elongatus). Seven predictive models were used x
including: LR (logistic regression), CT (classification trees), MARS (multivariate adaptative regression splines), ANN (artificial neural networks), DA (discriminant analysis), RF (random forest), and BR (boosted regression trees). Variables were derived using forward-selection procedures on five independent datasets and are shown as a percentage of datasets where each variable was selected (parentheses indicate negative associations). Variables selected from only one dataset were omitted. A prioi predictions (ap) based on habitat predictors thought to influence the decline of the species are shown for reference, where + indicates a positive correlation, - negative correlation, 0 none. ................................................................................. 108 Table 5.2 – A comparison of model sensitivity (correct classification of species presence), specificity (correct classification of species absence), and overall classification (correct classification of both species presence and absence) for redside dace (Clinostomus elongatus). Single models are: LR (logistic regression); CT (classification trees); BR (boosted regression trees); MARS (multivariate adaptive regression splines); ANN (artificial neural networks); DA (discriminant analysis); and, RF (random forest). Ensemble forecasts are: consensus model (CM); principal component analysis (PCA); weighted average using overall classification (WA); mean (Mn); and, median (Md). ................................................................................................... 116 Table 6.1 – Consensus measures of dendrogram group fidelity across distance measures (Euclidean and cosine) for each clustering algorithm: single linkage, unweighted pair group method with arithmetic means (UPGMA), and complete linkage. Group fidelity was determined by majority rules consensus tress using CI(C) consensus index................................................. 144
xi
List of Figures Figure 2.1 –Measuring the sensitivity of FD in a hypothetical eleven species community. The procedure consists of randomly dividing the community into two assemblages, noting how FD orders the two assemblages and assessing how sensitivity effects this order (i.e. did assemblage 1 always have higher FD given differences in methodology?). Each species is represented as a letter and the assemblages are represented as overlapping rectangles that contain the letters associated with their component species. The first set of rectangles represents one random division of the community into two assemblages. A new random division can be obtained by leaving the overlapping rectangles fixed and simply permuting the order of the species. The second set of overlapping rectangles gives an example of such a permutation. ........................... 24 Figure 2.2 - The effect of alpha and beta diversity on the probability of FD sensitivity for three communities crossed with four sets of construction methods. Darker shading represents a higher probability of sensitivity. The communities are A) Insectivorous birds (Holmes et al. 1979), B) Intertidal fish (Munoz and Ojeda 1997) and C) Predatory vertebrates (Jaksic and Medel 1990). Each column is for a different set of construction methods. For the first and second columns, overall PS and PG, all nine methods of dendrogram construction (three distance measures times three clustering algorithms) were used with the PS and PG methods respectively. For the third column, all three clustering algorithms were used with Gower’s distance and the PS method. For the fourth column, all three distance measures were used with UPGMA and the PS method. .... 27 Figure 2.3 –The effect of alpha and beta diversity on the probability of FD sensitivity for two communities crossed with two methods of FD calculation. Darker shading represents a higher probability of sensitivity. The communities are D) Patagonian forbs (Golluscio and Sala 1993) and E) Artic vegetation (Chapin et al. 1996). For these communities, only Gower’s distance could be calculated and so only three construction methods could be compared, corresponding to the three clustering algorithms. Each column is for a different method of FD calculation. The first is for the PS method and the second is for the PG method. .................................................. 28 Figure 3.1 – Example of rank-ordered, site-level vector residuals of Procrustean multivariate comparison. The length of a vector residual indicates an overall lack of fit for a site between two multivariate analyses. Shown is a comparison of full dataset of Principal Coordinates with xii
Jacaard’s distance and the same dataset where species occurring at 5% of sites were removed. Vectors shown in grey indicate those sites where at least one species was removed, whereas vectors in black indicate sites where no species were removed. The ratio of mean vector residuals between sites where species were removed versus those sites that did not have species removed indicates the distribution of impacts of the removal of rare species across multivariate analyses. Where most vector residuals for sites having species removed are largest, they indicate that these observations (sites) have been changed the most in their position between two ordinations..................................................................................................................................... 53 Figure 3.2 –Principal Coordinates Analysis (PCoA) of the sum of squares deviations (m2 statistic) comparing the concordance between solutions based on different ordination techniques, similarity coefficients and treatments of excluding rarely sampled species. A minimum-spanning tree was overlaid on Axes 1 and 2 to highlight connections between groups of points. Dashed lines indicate deviations from group membership in cases where clear groupings do not exist (e.g. M10 for Axes 2 and 3). Short forms are continued from Table 3.2. .................................... 55 Figure 3.3 – Site-level impact of the removal of rare species. Shown are box and whisker plots of the ratios of Procrustes vector residuals between sites for which rare species were removed and those sites that did not have any species removed. All comparisons were done by comparing site-level Procrustes vector residuals from the full datasets and with the removal of rare species across all similarity coefficients and ordination methods. ............................................................ 57 Figure 4.1 – Study sites on Rouge River, Ontario where redside dace (Clinostomus elongatus) dispersal and patch dynamics were monitored. Study locations: A) Leslie Tributary, and B) Berczy Creek, were sub-divided into extensive sites (black), where redside dace were tagged with a color-coded visual implant elastomer tag, and extended sites (grey), which were monitored for tag movement. ........................................................................................................ 76 Figure 4.2 – Metapopulation viability of the endangered species the redside dace (Clinostomus elongatus) in two stream metapopulations: A) Leslie Tributary, and B) Berczy Creek. Shown are the probabilities of extinction (y-axis) in years (x-axis) of a stochastic patch-based metapopulation model. Models were parameterized using: I) indirect parameterization of colonization and dispersal via patch distance, and; II) direct parameterization of colonization and dispersal using empirical estimates from a mark-recapture study. Legend: Vertical hashes xiii
represent a time interval of 100 years, solid lines indicate population trajectories where regional stochasticity was set to 0, dashed lines set to 0.1 and dotted lines set at 0.2. ............................... 85 Figure 4.3 – Differences in patch viability parameterized using indirect (y-axis) and direct (xaxis) patch dynamics of the endangered species the redside dace (Clinostomus elongatus) in two stream metapopulations: A) Leslie Tributary (L6-L15), and B) Berczy Creek (B6-B18). Shown are the mean probabilities of persistence of a given patch across 10,000 simulations. To demonstrate the variability in patch viability, 25% quantiles are overlaid as the negative of both the vertical and horizontal axes, while 75% quantiles are overlaid as the positive vertical and horizontal axes. The dashed line is a 1:1 line. .............................................................................. 86 Figure 5.1 – Distribution of sampling locations between 1997-2007. Closed circles indicate redside dace occurrences, whereas, open circles indicate where redside dace were absent. ...... 107 Figure 5.2 .................................................................................................................................... 110 Figure 5.3- Cluster analysis showing the relationship with the observed distribution (Obs.) of redside dace (Clinostomus elongatus) and: A) the seven individual modeling approaches alone, and B) with ensemble forecasts included. Model short forms are carried over from Table 5.2. 117 Figure 5.4 – Box and whisker plots showing the variability in consensus ensemble forecasts for predicting the presence of an endangered redside dace. Consensus ensemble models were compared across all combinations of one (n =7), three (n =35), five (n =21) and seven (n=1) input models (x-axis). Boxes are 25th and 75th percentiles, horizontal lines indicate the median, vertical lines indicate the upper and lower values, diamonds indicate the mean and are connected by dashed lines. Modeling metrics were: A) model sensitivity (i.e. correct classification of species presence); B) model specificity (i.e. correct classification of species absence); and, C) overall classification. Dashed lines indicate 95% confidence intervals. .................................... 119 Figure 6.1 - The relationship between distance measure (Euclidean or cosine) and clustering algorithm (SL = single linkage / nearest neighbor, UPGMA = unweighted pair group method with arithmetic mean, CL = complete linkage / farthest neighbor) with FD using five community data sets: A) Arctic vegetation (Chapin et al. 1996); B) Insectivorous birds (Holmes et al. 1979); C) Patagonian forbs (Golluscio and Sala 1993); D) Intertidal fish (Munoz and Ojeda 1997); and, E) Predatory vertebrates (Jaksic and Medel 1990). FD values were re-scaled relative to the xiv
Arctic vegetation data, which has the highest FD values. This standardization leads to the appearance of a constant outcome for the Arctic dataset, but this consistency is solely an artifact of using it as the reference point rather than the outcomes not differing depending on the resemblance measure or clustering algorithm............................................................................. 142 Figure 6.2 - The relationship between distance measure (solid lines = Euclidean distance, dashed lines = cosine distance) and building a dendrogram using a clustering algorithm (1 = complete linkage / farthest neighbor, 2 = unweighted pair group method with arithmetic mean / UPGMA, 3 = single linkage / nearest neighbor) where species are individually removed when calculating FD. Five community data sets are shown: A) Arctic vegetation (Chapin et al. 1996), B) Insectivorous birds (Holmes et al. 1979), C) Patagonian forbs (Golluscio and Sala 1993), D) Intertidal fishes (Munoz and Ojeda 1997), and E) Predatory vertebrates (Jaksic and Medel 1990). Shown inset are 50% majority rule consensus trees demonstrating lack of between group fidelity of species where calculating functional diversity using different distance measures, but the same clustering approach. .................................................................................................................... 143
xv
List of Appendices Appendix 2.1 – Derivation of Equation 1. .................................................................................... 36 Appendix 2.2 – MatLAB Code for testing sensitivity of FD ........................................................ 38 Appendix 2.3 – MatLAB Code for calculating FD via Podani and Schmera ............................... 41 Appendix 2.4 – MatLAB Code for calculating probabilities of sensitivity .................................. 42 Appendix 3.1 – Site-level effects of methodological choices in bioassessments. Shown are the ratios between mean site-level vector residuals from Procrustes analyses of sites having species removed and those sites having no species removed. Mean site-level vector residual values were separated by sites which had rare species removed (M1: n=2; M5: n=19; and M10: n=63); and compared with those sites that not. ............................................................................................... 71 Appendix 3.2 – Summary of three-dimensional ordination results. Shown are eigenvalues for Principal Coordinates Analyses (PCA) and Correspondence Analyses (CA), with percent variance explained shown in parentheses. Stress values are shown for Non-Metric Multidimensional Scaling (NMDS). ............................................................................................. 72 Appendix 4.1 –The endangered redside dace (Clinostomus elongatus). Photo credit: Mark Poos. ..................................................................................................................................................... 100 Appendix 4.2 –Visual implant elastomer (VIE) tag inserted subcutaneously on the ventral surface of the endangered redside dace (Clinostomus elongatus). Photo credit: Mark Poos. ................. 100 Appendix 4.3 – Model parameters used in the stochastic patch-occupancy models. Shown are the number of emigrants (mi), the mean probability of detection (PD), the number of emigrants adjusted for potentially missed tags (Mi), the number of immigrants needed to start a new subpopulation, and the rate of extinction (Ei). Other parameters include: the incidence function for Leslie Tributary (dI = 210, x = 0.4926, e = 3.685, y = 6.12), and Berczy Creek (dI = 150, x = 0.5652, e = 4.187, y = 7.01). Not shown: dij given it is a pairwise estimate rather than unique for each pool. .................................................................................................................................... 101
xvi
Appendix 5.1 – Model metrics for all combinations of consensus ensemble models. Models are: LR (logistic regression), CT (classification trees), MARS (multivariate adaptive regression splines), RF (random forest), ANN (artificial neural networks), BR (boosted regression trees, and DA (discriminant analysis). ........................................................................................................ 131 Appendix 5.2 – R Code for testing configurations of 1,3, 5 and 7 prediction consensus models ..................................................................................................................................................... 133 Appendix 6.1 – MatLAB Code for calculating total branch lengths of dendrograms from various species combinations .................................................................................................................. 152
xvii
1
Chapter 1
General Introduction “All models are wrong, some models are useful” George Box, Statistics Professor The overall aim of ecology is to broaden the understanding of the relationship between species and their environment (Krebs1998). This objective can be achieved through a myriad of ways from broadening the understanding of the structural components of an organism (e.g. molecular or cellular biology, physiology, evolutionary ecology), to understanding specific aspects of the organism itself (e.g. behavioral or evolutionary ecology, genetics), to understanding the interface between the organism and its biotic and abiotic environment (e.g. community or population ecology). Due to the complexity involved in each of the disciplines, and regardless of how ecologists view the world, simplification of natural processes are needed and a statistical model is often required. Statistical models provide the framework for an interpretation of nature. A statistical model is, by definition, merely a mathematical representation of some aspect of nature (Quinn and Keough 2002) and can take many forms. Models can be predictive in that they can be used to forecast, among other things, species distributions (e.g. Araújo and Guisan 2006; Elith et al. 2006), habitat importance (e.g. Olden and Jackson 2001; Guissan and Thuiller 2005), or potential changes in the environment and species response (e.g. climate change; Thuiller 2004; Araújo et al. 2005). Models can also be used to test hypotheses, e.g. is environment A better than B for a given species (e.g. Matthews et al. 1992; Grossman et al. 2002; Skyfield and Grossman 2008)? Lastly, models can be strictly informative, such as what is the current condition of environment A or species B (e.g. Grossman 1984; Mathews et al. 1982; 1986)? In all cases, models are assumed as useful if they can provide a realistic representation of the underlying natural world (Hilborn and Mangel 1997). However, many factors may influence models and understanding methodological choices may provide improved understanding of how ecologists interpret nature. The history of statistical methods has demonstrated that understanding the impact of methodological choices in today’s modeling approaches is sorely needed. Prior to the development of modern computers, ecologists opted for models where the structure was arbitrarily simplified (Quinn and Keough 2002). As such, there were relatively few
2
methodological choices, as data were often ‘massaged’ to fit the statistical approach or, alternatively, practitioners used their favorite method ad-hoc (Jackson 1993). Now, with complex statistical software and powerful computers readily available, ecologists have a plethora of approaches to choose from and many methodological choices to make. Understanding how methodological choices may impact results is an area of active and ongoing research (Jackson et al. 1989; Hilborn and Mangel 1997; Grossman et al. 2006). Over the last decade alone, methods for the analysis of species distribution have diversified and, at this point, dozens of different approaches are available (Guisan and Zimmermann 2000; Guisan and Thuiller 2005; Elith et al. 2006). Previous studies have demonstrated that methodological choices made at one level of analysis may have cascading impacts on subsequent choices (Jackson et al. 1989; Grossman et al. 1991; Dormann et al. 2008). For example, Dormann et al. (2008) demonstrated that changes in both data quality and choice of model type had large impacts for modeling the distribution of the critically endangered bird the Great Grey Shrike (Lanius excubitor). Further, in a review of modeling species distributions with presence only data, Elith et al. (2006) showed large differences in predictive performance among modeling methods at both regional and species levels. These examples highlight that, despite substantial effort in improving statistical methods for conservation-based species, there remains great uncertainty regarding the impacts of methodological choices. There are numerous debates in ecological literature regarding methodological choices in diverse areas such as how to measure functional traits to account for ecosystem function (Petchey and Gaston 2002; Podani and Schmera 2006), the impact of including rare species in biological assessments (Cao et al. 1998; Marchant 1999), and whether species-specific dispersal can be measured using distance-based functions (Heinz et al 2005; 2006). Given the choices that ecologists must make, knowing the pitfalls and assumptions of conservation-based models become increasingly important. For example, ecologists should know the impact of methodological choices given the sampling design and kind of data that they wish to collect (Quinn and Keough 2002). My thesis attempts to address issues in methodological choices in conservation-based models in two ways. In the first section of the thesis, the impacts of methodological choices inherent in conservation-based models will be examined across a broad selection of available approaches
3
from measuring functional diversity (Chapter 2), to developing bio-assessments in community ecology (Chapter 3), to assessing dispersal in metapopulation analyses (Chapter 4). It is the goal of this section to illustrate the potential for methodological choices to impact conservation-based models, regardless of the scale, study system or species involved. For clarity, a conservation-based model refers to a model used to interpret/advance knowledge of species with conservation concern (e.g. conservation designation, rare, or declining). These models can include population viability analyses, models of habitat suitability or models of ecological function. A methodological impact refers to a difference in conclusions (e.g. population viability, habitat suitability, and ecological function) derived from methodological choices that alter ecological interpretation and/or lead to differences in recovery/management. Species with conservation concern refers to those which occur infrequently, have a conservation status (e.g. endangered or threatened) or are declining (i.e. undergoing reductions in population size but do not have a conservation designation). As there are many forms of rarity (Rabinowitz et al. 1986), the definition of conservation-based models and species with conservation concern used in this thesis is purposefully inclusive to allow for the greatest possible impact on conservation applications. In the second section of this thesis, the use of consensus methods will be demonstrated as a potential tool for reducing uncertainty with methodological choices in conservation-based models. Two separate applications of consensus methods will be highlighted, including how consensus methods can reduce uncertainty from choosing a modeling type (Chapter 5) and can be used to identify when methodological choices may be a problem (Chapter 6).
I) Impacts of Methodological Choices in Conservation-Based Models To some researchers, there is perhaps no greater goal of ecology than to preserve the biological diversity of species (Ehlrich and Wilson 1991). The loss of biological diversity through human impact is of concern for the structure of ecological communities and its affects may, in turn, affect the structure of ecosystems (e.g. keystone species). Current rates of decline of species worldwide are thought to be hundreds to thousands of times higher than pre-human levels (Jelks et al. 2008). In particular, freshwater systems may have rates of decline many times higher than terrestrial systems (Ricciardi and Rasmussen 1999). With the rate of species decline increasing, the need for conservation-based models that reflect biological phenomena is becoming timely.
4
There are many methodological difficulties with developing conservation-based models and understanding these challenges may help improve species management. First, species that are of conservation concern suffer from low sample sizes and spatially fragmented distributions, thereby making accurate enumeration difficult (Green and Young 1993) and traditional statistical approaches often invalid (Ellison and Agrawal 2005). As a result, models fitted to a response of species with conservation concern usually have low power or high uncertainty (Cunningham and Lindenmayer 2005). Further, in many multivariate methodologies, researchers down weight or exclude rare species altogether (Gauch 1982), thereby reducing or eliminating their usefulness for conservation applications. The role of methodological choices has been debated in several types of conservation-based models. In the following section a few ecological examples will be highlighted for further analysis in this thesis.
Functional Diversity Given the current loss of biodiversity, it is important to be able to model how functional diversity will respond to species loss. Quantification of the relationship between species loss and functional diversity is necessary to highlight unique species traits that may be lost and their implication to ecosystem function (Srivastava and Vellend 2005). Extinction should not affect overall ecosystem function if all species have similar functions, but it will have a major effect if each species has a different function (Fonseca and Ganade 2001). Simulations using natural communities have shown that the loss of species reduces functional diversity disproportionately relative to the number of species (Petchey and Gaston 2002; Larsen et al. 2005). These findings are in agreement with empirical studies noting the loss of rare species disproportionately impacts ecosystem function (Hooper et al. 2005). One of the current metrics of functional diversity (sensu Petchey and Gaston 2002) uses the total branch length of a dendrogram (tree) constructed from a matrix of ecological traits. This method requires three steps before the total branch lengths can be calculated, including: 1) obtaining a trait matrix; 2) converting the trait matrix into a distance matrix; and, 3) clustering of the distance matrix to produce a dendrogram. Petchey and Gaston (2002) have suggested that the particular analytical methods used to produce the functional dendrogram are of limited relevance
5
to the resultant metric as the relationship between functional diversity and species richness is robust to changes in distance metric and clustering algorithm. However, previous research has indicated that changes in distance metric have dramatic effects on clustering dendrograms (Jackson et al. 1989; Jackson 1993) and, similarly, choice of clustering algorithm (Legendre and Legendre 1998) can alter biological conclusions. The role of methodological choices for altering functional diversity (FD; Petchey and Gaston 2002) has only been recently debated (Petchey and Gaston 2006, 2007; Podani and Schmera 2006, 2007; Poos et al. 2009). Functional diversity is an important component of conservationbased models as researchers are often interested in the amount of functional variation with an ecosystem (Loreau et al. 2001). Comparisons of the amount of functional variation that is lost when species have been removed from an ecosystem has demonstrated that rare species play an important role for maintaining ecosystem function (Fonseca and Ganade 2001; Petchey and Gaston 2002b; Schmera et al. 2009). In Chapter 2, the importance of methodological choices on estimating functional diversity is examined. In particular, Chapter 2 highlights whether how functional diversity is measured can alter ecological insight.
Bioassessments Similar to functional diversity, there has been a dramatic increase in the amount of bioassessment literature in the past few decades (Bailey et al. 2004). Bioassessments have been used by managers for decades to evaluate communities undergoing anthropogenic impact and species decline (Green 1979; Barbour et al. 1999). Debates regarding the impact of methodological choices on the assessment of biological communities (i.e. bioassessment) are varied and include criticisms of multi-metric approaches (Hannaford and Resh 1995; Wallace et al. 1996; Bowman and Somers 2006); the development of multivariate and predictive models (Bailey et al. 2004), taxonomic resolution (Somers et al. 1998; Hewlett 2000), improvements to rapid methods (Hannaford and Resh 1995; Somers et al. 1998), redundancy in metrics (Barbour et al. 1992), and improvements to the study design of reference conditions (Norris and Thoms 1999; Bowman and Somers 2005; Bailey et al. 2009). One ongoing debate in bioassessments is the importance of rare species (Cao et al. 1998; 2001; Marchant 1999; Marchant 2002). Methods for identifying trends in multivariate analyses which
6
include rare species are often limited because they typically represent points of unusually high leverage in the analysis. When rare species (i.e. species that occur infrequently) are included in multivariate analyses, they often alter the analysis because, relative to more common species, they over-fit an association of relatively few occurrences to a specific habitat type that, in turn, unduly influences species-habitat differences (Legendre and Gallagher 2001). To compensate researchers either down weight or exclude rare species from multivariate analyses without consideration of the degree of influence they may have on the underlying biological data (Gauch 1982). Obviously, the benefits of this decision are negated when the species or relationship of interest relates to a species that is rare, as is the case with many conservation applications. Imperiled species (i.e. species with a conservation status), like rare species, are scarce and often found in small numbers in specific habitat types (Mace 1994). In Chapter 3, the impact of removing rare species is assessed relative to other multivariate methodological choices, such as the choice of distance measure and multivariate method, which are all common to multivariate bioassessment approaches.
Metapopulation Viability Models Trajectories of populations through time are often needed to extrapolate population viability. Approaches such age or stage based matrix population models require parameterization of several species life history characteristics (e.g. fecundity, age/size of maturation), which may be limited or unknown for species with conservation concern. As such, simplifications in how many life history characteristics are needed to model population trajectories are needed. Metapopulation viability analyses offer simplification in over traditional population viability analyses as they only require parameterization of patch occupancy and estimates of rates of colonization and extinction (Hanski 1999). One popular method of metapopulation viability analysis is stochastic patch occupancy models (SPOM, Hanksi 1999; Moilanen 1999; 2004). Stochastic patch occupancy models provide a spatially realistic model for patch dynamics as they incorporate estimates of patch location and species dispersal into the modeling procedure. For example, estimates of patch colonization (i.e. “reachability”) are quantified using the distance between a starting patch and a target patch, and the ability of species to disperse (Hanksi 1994; Hanski et al. 1996; Heniz et al. 2005). Most often, this measure is estimated by assuming that colonization potential declines exponentially with distance (i.e. exponential decay; Hanski 1994;
7
Vos et al. 2001; Frank and Wissel 2002). However the ability of how well even simple formulae are able to model patch dynamics remains an open question (Heinz et al. 2006; Marsh 2008). Chapter 4 aims to better understand the impact of estimating population viability using direct versus indirect parameterization of dispersal in metapopulation models.
II) Using Consensus Methods to Reduce Uncertainty from Methodological Choices Conservation-based models provide a useful tool for understanding consequences of species loss (e.g. functional diversity), anthropogenic impacts (e.g. bioassessments) and the predictions for population viability and habitat suitability. However, the difficulty with using the standard sets of conservation-based models is that they are often inappropriate when analyzing data limited to a few sites and/or species (Ellison and Agrawal 2005). This condition of sparse data, in turn, contributes to many statistical limitations, including zero-inflated bias, increased collinearity between variables and inflated coefficient of variation (Dixon et al. 2005; Edwards et al, 2006; Dormann et al. 2008). However, uncertainties may arise during all stages of modeling such as obtaining species level data obtaining accurate species counts, and linking species to landscapes/habitats (Burgman et al. 2005; Cunningham and Lindenmayer 2005; Edwards et al. 2005; Elith et al. 2006). Therefore, developing methods that may reduce impacts from methodological choices may enhance the utility of conservation-based models. The use of consensus methods has been highlighted as one way to address problems related to using conservation-based models (Araújo and New 2007; Marrimon et al. 2009). Laplace (1820) demonstrated that the probability of error will decline as more models are included (i.e. an ensemble or consensus approach). This old idea has only recently been applied to problems of uncertainty with the influential work of Bates and Granger (1969) that in part, contributed to the award of the Nobel Memorial Prize in Economics to Clive Granger in 2003. Recently, ecologists have applied consensus models to issues with uncertainty in climate-based predictions (Thuiller 2004; Araújo et al. 2005; Thuiller et al 2005; Buisson and Grenouillet 2009) and, to some extent, to functional diversity (Mouchet et al. 2008). Yet, the application of consensus methods seems due to other problems in ecology.
8
The application of consensus methods to resolve issues of methodological choices is in its infancy. Issues such as how to build the best ensemble model have yet to be quantitatively evaluated. The second section of the thesis applies consensus models to impacts from methodological choices including whether consensus methods can improve prediction over several singular methodological approaches (Chapter 5) and whether consensus methods can identify instances where methodological choices may be an issue as in the measure of functional diversity (Chapter 6). In Chapter 5, a quantitative evaluation of consensus models is applied to the habitat modeling of the endangered fish, Redside Dace (Clinostomus elongatus). Chapter 5 aims to evaluate the ability of consensus models to improve model performance (e.g. model specificity, sensitivity and overall classification). In Chapter 6, consensus methods are used to evaluate whether or not methodological choices are an issue for measuring functional diversity.
Scope The difficulty with developing any model for conservation-based purposes are the linkages needed across many facets of ecological knowledge, from how species traits relate to the environment, to the types of environments undergoing impact. These relationships require knowledge not only of the species, but each system, the suite of species found there, the relationships between those species and their related habitats, and the relationship between each habitat and large scale (e.g. landscape) and small scale (e.g. in-stream habitat) ecological processes (Smith and Powell 1971; Matthews 1998; Jackson and Harvey 1989; Jackson et al. 2001). Therefore, there is no easy starting point (e.g. scale, method or species) for determining modeling approaches appropriate for conservation-based issues. The scope for identifying the influence of methodological choices in conservation-based models is broad by its very nature. As such, this thesis is not restricted to a given study organism, study system, data series, or methodology. Instead, data is used that best addresses the specific questions and objectives of the various chapters. For example, in Chapters 2, 3 and 6 data is from published studies to demonstrate the impact of methodological issues. In Chapter 4, empirical information was gathered over the course of a one-year field study. In Chapter 5, over 420,000 fish records were gathered from various not-for-profit and government agencies, universities and other academic institutions and studies.
9
In addition, this thesis examines different spatial and temporal scales for identifying methodological choices. In Chapters 4 and 5, data was collected on dispersal and habitat suitability of Redside Dace (Clinostomus elongatus). In Chapter 4, the spatial scale is the province of Ontario and temporal scale are data points within the last decade. In Chapter 5, the spatial scale is one watershed and the temporal scale is one year. These spatial and temporal scales were chosen to highlight the approaches being compared (e.g. predictive models in Chapter 4, and metapopulation models in Chapter 5). The goal of this thesis is that through the consistency of the results within this thesis, across study organisms, study systems, data series, methodologies, and spatial and temporal scales, it will provide ecological insight into making methodological decisions in ecology. Each chapter is linked based on the influence of methodological choices and an overall interest on how methodological choices can alter conclusions in ecological studies.
Statement of Contribution Parts of this thesis would not have been possible without the contribution of many collaborators. Chapters 2 and 6 were completed in collaboration with Steve Walker, who helped in the derivation of Equation 1 (Appendix 2.1), the development of ideas, and in writing code for assessing the robustness of the functional diversity measure. Chapter 4 was completed with the help of staff from the Toronto Region Conservation Authority, which provided valuable field and logistical assistance. Chapter 5 was completed with the help of various governmental, nongovernmental and not-for profit agencies who contributed data.
Publication of Thesis Material Chapter 2 has been published and can be found as: Poos, M.S., S.C. Walker and D.A Jackson. 2009. Functional diversity indices can be driven by methodological choices and species richness. Ecology 90(2): 341-347. (doi: 10.1890/08-1638.1) Chapters 3 and 5 are currently submitted and undergoing peer review. All published material is provided under approved copyright from the related journal.
10
References Araújo, M. B., and A. Guisán. 2006. Five (or so) challenges for species distribution modelling. Journal of Biogeography 33:1677-1688. Araújo, M. B., and M. New. 2007. Ensemble forecasting of species distributions. Trends in Ecology and Evolution 22:42-47. Araújo, M. B., R. J. Whittaker, R. J. Ladle, and M. Erhard. 2005. Reducing uncertainty in projections of extinction risk from climate change. Global Ecology and Biogeography 14:529-538. Bailey, R. C., R. H. Norris, and T. B. Reynoldson. 2004. Bioassessment of freshwater ecosystems: Using the reference condition approach. Springer, New York. Bailey, R. C., M. G. Kennedy, M. Z. Dervish, and R. M. Taylor. 2008. Biological assessment of freshwater ecosystems using a reference condition approach: comparing predicted and actual benthic invertebrate communities in Yukon streams. Freshwater Biology 39:765774. Barbour, M.T., J. Gerritsen, B.D. Snyder, and J.B. Stribling. 1999. Rapid bioassessment protocols for use in streams and wadeable rivers: Periphyton, benthic macroinvertebrates and fish. 2nd edition. EPA-841-B-99-002. Environmental Protection Agency. Bates, J. M., and C. W. Granger. 1969. The combination of forecasts. Operations Research Quarterly 20:451-468. Bowman, M. F., and K. M. Somers. 2005. Considerations when using the reference condition approach for bioassessment of freshwater ecosystems. Water Quality Research Journal of Canada 40:347-360. Bowman, M. F., and K. M. Somers. 2006. Evaluating a novel Test Site Analysis (TSA) bioassessment approach. Journal of the North American Benthological Society 25:712727.
11
Burgman, M., D. B. Lindenmayer, and J. Eltih. 2005. Managing landscapes for conservation under uncertainty. Ecology 86: 2007-2017. Buisson, L., and G. Grenouillet. 2009. Contrasted impacts of climate change on stream fish assemblages along an environmental gradient. Diversity and Distributions 15:613-626. Cao, Y., D.D. Williams, and N.E. Williams. 1998. How important are rare species in community ecology and bioassessment? Limnology and Oceanography 43: 1403–1409. Cao, Y., D.P. Larsen, and R.S. Thorne. 2001. Rare species in multivariate analysis for bioassessment: some consideration. Journal of the North American Benthological Society 20: 144–153. Cunningham, R.B., and D.B. Lindenmayer. 2005. Modeling count data of rare species: Some statistical issues. Ecology 86: 1135-1142. Dixon, P.M., A.M. Ellison, and N.J. Gotelli. 2005. Improving the precision of estimates of the frequency of rare events. Ecology 86: 1114-1123. Dormann, C. F., O. Purschke, J. R. G. Márquez, S. Lautenbach, and B. Schröder. 2008. Components of uncertainty in species distribution analysis: A case study of the great grey shrike. Ecology 89:3371-3386. Edwards, T.C., D.R. Cutler, N.E. Zimmerman, L. Geiser, and J. Alegriae. 2005. Model-based stratifications for enhancing the detection of rare ecological events. Ecology 86: 10811090. Ehlrich, P.R. and E.O. Wilson. 1991. Biodiversity studies: Science and policy. Science 253: 758763. Elith, J., C. H. Graham, R. P. Anderson, M. Dudík, S. Ferrier, A. Guisan, R. J. Hijmans, F. Huettmann, J. R. Leathwick, A. Lehmann, J. Li, L. G. Lohmann, B. A. Loiselle, G. Manion, C. Moritz, M. Nakamura, Y. Nakazawa, J. M. Overton, A. T. Peterson, S. J. Phillips, K. Richardson, R. Scachetti-Pereira, R. E. Schapire, J. Soberón, S. Williams, M. S. Wisz, and N. E. Zimmermann. 2006. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29:129-151.
12
Ellison, A.M., and A.A. Agrawal. 2005. The statistics of rarity. Ecology 86: 1079-1080. Frank, K., and C. Wissel. 2002. A formula for the mean lifetime of metapopulations in heterogeneous landscapes. The American Naturalist 159:530-552. Fonseca, C. R., and G. Ganade. 2001. Species functional redundancy, random extinctions and the stability of ecosystems. Ecology 89:118-125. Gauch, H. G. 1982. Multivariate analysis in community ecology. Cambridge University Press, London, England. Green, R.H. and R.C. Young. 1993. Sampling to detect rare species. Ecological Application 3: 351-366. Grossman, G. D., and V. Devlaming. 1984. Reproductive ecology of female Oligocottus snyderi greeley - A North American intertidal sculpin. Journal of Fish Biology 25:231-240. Grossman, G. D., D. M. Nickerson, and M. C. Freeman. 1991. Principal component analyses of assemblage structure data - Utility of tests based on eigenvalues. Ecology 72:341-347. Grossman, G. D., P. A. Rincon, M. D. Farr, and R. E. Ratajczak. 2002. A new optimal foraging model predicts habitat use by drift-feeding stream minnows. Ecology of Freshwater Fish 11:2-10. Grossman, G. D., R. E. Ratajczak, J. T. Petty, M. D. Hunter, J. T. Peterson, and G. Grenouillet. 2006. Population dynamics of mottled sculpin (Pisces) in a variable environment: information theoretic approaches. Ecological Monographs 76:217-234. Guisan, A., and W. Thuiller. 2005. Predicting species distribution: offering more than simple habitat models. Ecology Letters 8: 993-1009. Guisan, A., and N. E. Zimmermann. 2000. Predictive habitat distribution models in ecology. Ecological Modelling 135:147-186.
13
Hannaford, M. J., and V. H. Resh. 1995. Variability in macroinvertebrate rapid-bioassessment surveys and habitat assessments in a Northern California stream. Journal of the North American Benthological Society 14:430-439. Hanski, I. 1999. Metapopulation ecology. Oxford University Press, New York, New York. Hanski, I. 1994. A practical model of metapopulation dynamics. The Journal of Animal Ecology 63:151-162. Hanski, I., A. Moilanen, T. Pakkala, and M. Kuussaari. 1996. The quantitative incidence function model and persistence of an endangered butterfly metapopulation. Conservation Biology 10:578-590. Heinz, S. K., L. Conradt, C. Wissel, and K. Frank. 2005. Dispersal in fragmented landscapes: Deriving a practical formula for patch accessibility. Landscape Ecology 20:83-99. Heinz, S. K., C. Wissel, and K. Frank. 2006. The viability of metapopulations: individual dispersal behaviour matters. Landscape Ecology 21:77-89. Hewlett, R. 2000. Implications of taxonomic resolution and sample habitat for stream classification at a broad geographic scale. Journal of the North American Benthological Society 19:352-361. Hilborn, R. and M. Mangel. 1997. The ecological detective: confronting models with data. Princeton University Press, Princeton, New Jersey, 330 pp. Hooper, D. U., F. S. Chapin, J. J. Ewel, A. Hector, P. Inchausti, S. Lavorel, J. H. Lawton, D. M. Lodge, M. Loreau, S. Naeem, B. Schmid, H. Setala, A. J. Symstad, J. Vandermeer, and D. A. Wardle. 2005. Effects of biodiversity on ecosystem functioning: A consensus of current knowledge. Ecological Monographs 75:3-35. Jackson, D.A. 1993. Multivariate analysis of benthic invertebrate communities: The implication of choosing particular data standardizations, measures of association, and ordination methods. Canadian Journal of Fisheries and Aquatic Sciences 50:2641-2651.
14
Jackson, D.A. and H.H. Harvey. 1993. Fish and benthic invertebrates: Community concordance and community-environment relationships.Canadian Journal of Fisheries and Aquatic Sciences 50: 2641-2651. Jackson D. A., K.M. Somers, and H.H. Harvey. 1989. Similarity coefficients: measures of cooccurrence and association or simply measures of occurrence? American Naturalist 133: 436453. Jackson, D.A., P.R. Peres-Neto, and J.D. Olden. 2001. What controls who is where in freshwater fish communities – the roles of biotic, abiotic, and spatial factors. Canadian Journal of Fisheries and Aquatic Sciences 58: 157-170. Jelks, H. L., S. J. Walsh, N. M. Burkhead, S. Contreras-Balderas, E. Diaz-Pardo, D. A. Hendrickson, J. Lyons, N. E. Mandrak, F. McCormick, J. S. Nelson, S. P. Platania, B. A. Porter, C. B. Renaud, J. J. Schmitter-Soto, E. B. Taylor, and M. L. Warren. 2008. Conservation status of imperiled North American freshwater and diadromous fishes Fisheries 33:372-386. Krebs, C.J.1998. Ecological Methodology (2nd ed.). Addison-Welsey Educational Publishers Inc., New York, New York. Laplace, P. S. 1820. Théorie analytique des probabilités. Courcier, Paris. Larsen, T. H., N. M. Williams, and C. Kremen. 2005. Extinction order and altered community structure rapidly disrupt ecosystem functioning. Ecology Letters 8:538-547. Legendre, P. and E.D. Gallagher. 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129: 271-280 Legendre, P. and L. Legendre. 1998. Numerical Ecology (2nd edition). London, Elsevier. Loiselle, B. A., C. A. Howell, C. H. Graham, J. M. Goerck, T. Brooks, K. G. Smith, and P. H. Williams. 2003. Avoiding pitfalls of using species distribution models in conservation planning. Conservation Biology 17:1591-1600.
15
Loreau, M., S. Naeem, P. Inchausti, J. Bengtsson, J. P. Grime, A. Hector, D. U. Hooper, M. A. Huston, D. Raffaelli, B. Schmid, D. Tilman, and D. A. Wardle. 2001. Biodiversity and ecosystem functioning: current knowledge and future challenges Science 294: 804-808. Mace, G.M. 1994. Classifying threatened species: means and ends. Philosophical Proceedings of the Royal Society of London. Series B. 344: 91-97. Matthews, W. J., E. Surat, and L. G. Hill. 1982. Heat death of the orangethroat darter Etheostoma spectabile (Percidae) in a natural environment. Southwestern Naturalist 27:216-217. Matthews, W. J., M. E. Power, and A. J. Stewart. 1986. Depth distribution of Campostoma grazing scars in an Ozark stream. Environmental Biology of Fishes 17:291-297. Matthews, W. J., F. P. Gelwick, and J. J. Hoover. 1992. Food and habitat use by juveniles of species of Micropterus and Morone in a southwestern reservoir. Transactions of the American Fisheries Society 121:54-66. Matthews, W. J. 1998. Patterns in freshwater fish ecology. Chapman and Hall, New York, New York. Marchant, R. 1999. How important are rare species in aquatic ecology and bioassessment? A comment on the conclusions of Cao et al. Limnology and Oceanography 44: 1840–1841. Marchant, R. 2002. Do rare species have any place in multivariate analysis for bioassessment? Journal of the North American Benthological Society 21: 311-313 Marmion, M., M. Parviainen, M. Luoto, R. K. Heikkinen, and W. Thuiller. 2009. Evaluation of consensus methods in predictive species distribution modelling. Diversity and Distributions 15: 59-69. Marsh, D. 2008. Metapopulation viability analysis for amphibians. Animal Conservation 11:463465. Moilanen, A. 1999. Patch occupancy models of metapopulation dynamics: Efficient parameter estimation using implicit statistical inference. Ecology 80:1031-1043.
16
Moilanen, A. 2004. SPOMSIM: software for stochastic patch occupancy models of metapopulation dynamics. Ecological Modelling 179:533-550. Mouchet, M., F. Guilhaumon, S. Villéger, N. W. H. Mason, J.-A. Tomasini, and D. Mouillot. 2008. Towards a consensus for calculating dendrogram-based functional diversity indices Oikos 117:794-800. Norris, R. H., and M. C. Thoms. 1999. What is river health? Freshwater Biology 41:197-209. Olden, J. D., and D.A. Jackson. 2001. Fish habitat relationships in lakes: gaining predictive and explanatory insight by using artificial neural networks. Transactions of the American Fisheries Society 130:878-897. Petchey, O. L. and K. J. Gaston. 2002. Functional diversity (FD), species richness and community composition. Ecology Letters 5: 402-411. Petchey, O. L., and K. J. Gaston. 2002b. Functional diversity (FD), species richness and community composition. Ecology Letters 5:402-411. Petchey, O. L., and K. J. Gaston. 2006. Functional diversity: back to basics and looking forward. Ecology Letters 9:741-758. Petchey, O. L., and K. J. Gaston. 2007. Dendrograms and measuring functional diversity. Oikos 116:1422-1426. Podani, J. and D. Schmera. 2006. On dendrogram-based measures of functional diversity. Oikos 115: 179-185. Podani, J., and D. Schmera. 2007. How should a dendrogram-based measure of functional diversity function? A rejoinder to Petchey and Gaston. Oikos 116:1427-1430. Quinn, G.P.and M.J. Keough 2002. Experimental design and data analysis for biologists. Cambridge University Press: London, 537 pp.
17
Rabinowitz, D., S. Cairns and T. Dillon. 1986. Seven forms of rarity, and their frequency in the flora of the British Isles. In Conservation biology: the science of scarcity, and diversity (ed. M.E. Soule), pp. 184-204. Sunderland, Massachusetts, Sinauer Associates. Ricciardi, A. and J.B. Rasmussen. 1999. Extinction rates of North American freshwater fauna. Conservation Biology 13: 1220-1222. Schmera, D., J. Podani, and T. Eros. 2009. Measuring the contribution of community members to functional diversity. Oikos 118: 961-971. Sharma, S. and D.A. Jackson. 2008. Predicting smallmouth bass incidence across North America: Comparison of statistical approaches. Canadian Journal of Fisheries and Aquatic Sciences 65: 471-481. Skyfield, J. P., and G. D. Grossman. 2008. Microhabitat use, movements and abundance of gilt darters (Percina evides) in southern Appalachian (USA) streams. Ecology of Freshwater Fish 17:219-230. Somers, K. M., R. A. Reid, and S. M. David. 1998. Rapid biological assessments: how many animals are enough? Journal of the North American Benthological Society 17:348-358. Srivastava, D. S., and M. Vellend. 2005. Biodiversity-ecosystem function research: Is it relevant to conservation? Annual Review of Ecology, Evolution, and Systematics 36: 267-294. Thuiller, W. 2004. Patterns and uncertainties of species’ range shifts under climate change. Global Change Biology 10: 2020-2027. Thuiller, W., S. Lavorel, M. B. Araújo, M. T. Sykes, and I. C. Prentice. 2005. Climate change threats to plant diversity in Europe. Proceedings of the National Academy of Sciences 102:8245-8250. Wallace, J. B., J. W. Grubaugh, and M. R. Whiles. 1996. Biotic indices and stream ecosystem processes: results from an experimental study. Ecological Applications 6:140-151. Vos, C. C., J. Verboom, P. F. M. Opdam, and C. J. F. ter Braak. 2001. Toward ecologically scaled landscape indices. American Naturalist 157: 24-41.
18
Section I:
The importance of methodological choices in ecological models
19
Chapter 2: The importance of methodological choices in influencing the measure of functional diversity across ecological communities
Abstract Functional diversity is an important concept in community ecology because it captures information on functional traits absent in measures of species diversity. One popular method of measuring functional diversity is the dendrogram-based method, FD. To calculate FD, a variety of methodological choices are required and it has been debated about whether biological conclusions are sensitive to such choices. We studied the probability that conclusions regarding FD were sensitive, and that patterns in sensitivity were related to alpha and beta components of species richness. We developed a randomization procedure which iteratively calculated FD by assigning species into two assemblages and calculating the probability that the community with higher FD varied across methods. We found evidence of sensitivity in all five communities we examined, ranging from a probability of sensitivity of 0 (no sensitivity) to 0.976 (almost completely sensitive). Variations in these probabilities were driven by differences in alpha diversity between assemblages and not by beta diversity. Importantly, FD was most sensitive when it was most useful (i.e. when differences in alpha diversity were low). We demonstrate that trends in functional diversity analyses can be largely driven by methodological choices or species richness, rather than functional trait information alone.
Keywords: multivariate statistics, functional diversity, community ecology, species richness, biodiversity, ecological organization, dendrogram.
20
Introduction Functional diversity is the amount of inter-specific variation in functional traits in an ecological community. The concept of functional diversity has received considerable attention recently, largely because of the following intuitive argument. Species diversity indices treat all species identically, whereas functional diversity indices do not. Therefore, it is reasonable to expect that functional diversity is likely to be more ecologically relevant because species differ from one another in functionally important ways (Petchey and Gaston 2002). For example, several studies have concluded that measures of ecosystem function tend to correlate more strongly with functional diversity indices than with species diversity indices (Loreau et al. 2001). These studies have spurred continued interest in developing new and improved functional diversity indices (Mouchet et al. 2008; Villeger et al. 2008). Despite the conceptual simplicity of functional diversity, ecologists wishing to measure it must choose from a number of approaches. Mason (2005) developed a typology of functional diversity indices with three types: functional richness; functional evenness; and, functional divergence. This typology is similar in spirit to the distinction between species richness and evenness in species diversity studies. For example, functional richness indices measure the amount of trait space filled by the species in a community, whereas functional evenness indices measure the evenness in the distribution of abundance in trait space (Mason 2005; Villeger et al. 2008). Using rarefaction techniques, functional richness and evenness can also be thought of as extremes along a gradient of functional diversity indices (Walker et al. 2008). Rarefaction also makes clear the close relationship between species and functional richness. It is therefore important to ensure that accepted indices of functional richness provide information beyond that of species richness, as data on functional traits can be costly to obtain. One approach to measuring functional richness, which has shown promise as a predictor of ecosystem function, is the dendrogram-based approach known as FD (Petchey and Gaston 2002). This approach consists of measuring functional richness as the total branch-length of a dendrogram that clusters species based on the similarity of their functional-trait characteristics. There are numerous methods for constructing a dendrogram; in particular, both a resemblance measure, which measures the difference between species in their functional-trait characteristics,
21
and a clustering algorithm, which specifies the manner in which similar species are grouped together, must be chosen. There is the possibility that ecological conclusions drawn from an analysis of dendrogram-based functional diversity may be sensitive to the methodological choices that are required for producing a dendrogram. This may be a serious issue given that dendrogram topology may change considerably with changes in the methods used (e.g. Sneath and Sokal 1973, Jackson et al. 1989). There has been considerable recent debate about the importance of the method of dendrogram construction for the measurement of dendrogram-based functional diversity (Petchey and Gaston 2006, 2007, Podani and Schmera 2006, 2007; Mouchet et al. 2008). To facilitate resolution, we conducted a detailed analysis of the sensitivity of dendrogram-based functional diversity measures to differences in species richness and methodological choices. For this resolution, we need a quantitative definition of sensitivity. Given a pair of species assemblages and set of dendrogram construction methods, we make the following definition: conclusions are insensitive if all construction methods result in the same assemblage being identified as having higher functional diversity. Conclusions are sensitive if at least one construction method identifies a different assemblage as having higher functional diversity. With this definition, we seek answers to the following questions. First, through a systematic study of previously analyzed data from ecological communities (Petchey and Gaston 2002, Podani and Schmera 2006, Petchey and Gaston 2007), what is the probability that conclusions regarding FD are sensitive to methodological choices? Second, if sensitivity is found to be likely in many communities analyzed, is the probability of sensitivity related to the difference in local species richness (i.e. alpha diversity) between the two assemblages? We hypothesize that the probability of sensitivity should be low when differences in alpha diversity are very high. Intuitively, we expect species richness to drive functional richness patterns in these cases, no matter how it is measured. This is a null hypothesis; it assumes that functional richness (as measured by FD) does not provide information beyond that provided by species richness. Failure to reject this hypothesis would suggest that FD and alpha diversity are largely redundant, provided that the most species rich assemblage also tends to have the highest FD. Third, is the probability of sensitivity related to the amount of species turnover (i.e. beta diversity) between assemblages? As species turnover can be measured in numerous ways, hereafter we use the term beta diversity to refer to Lande’s species turnover (Lande 1996). We hypothesize that the
22
probability of sensitivity should be high when beta diversity is high. High beta diversity will tend to lead to lower redundancy across traits between assemblages, in comparison to low beta diversity. Therefore, high beta diversity produces conditions under which we intuitively expect small differences in functional diversity. Small differences will presumably be more sensitive to methodological choices. Fourth, is the probability of sensitivity related to certain types of methodological choices? We hypothesize that conclusions will be more sensitive to the choice of distance measure than to the choice of clustering algorithm, because the distance measure can completely change the order of functional similarity amongst the species whereas the clustering algorithm is more limited in that it can only alter how groups of species relate to one another in multivariate space. We note that there are reasons to believe that FD will also be quite sensitive to the choice of a clustering algorithm. Indeed, different clustering algorithms can generate quite different tree topologies, which may translate into FD sensitivity. We address these questions by assessing the probability of sensitivity of pairs of randomly drawn sub-assemblages from five ecological communities.
Methods All of these analyses were based on data from ecological communities obtained from the literature. We used the same five datasets used in previous studies of FD (Petchey and Gaston 2002, 2007, Podani and Schmera 2006). These datasets represent variation in the number (from 13 to 37) and type of species, and the number and type of functional traits (from 6 to 27). For example, the three vertebrate datasets use characteristics ranging from foraging behavior to the consumption of prey species as their functional traits (Holmes et al. 1979, Jaksic and Medel 1990, Munoz and Ojeda 1997), whereas the remaining two datasets rely on vegetative characteristics, such as rooting depth and herbivore palatability, of the plants being studied (Golluscio and Sala 1993, Chapin et al. 1996). The general approach to assessing the sensitivity of FD to methodological choices was as follows (see Fig. 2.1 for an example). For each community (i.e. dataset), we organized all of the species, γ, into two groups, hereafter referred to as assemblages. Let the average species richness over the two assemblages be α . Each species in the community was included in either one of the assemblages or in both. For a given level of beta diversity, β = γ − α , and difference in alpha diversity between the assemblages, Δα, the total number of unique pairs of assemblages is
23
. (1)
The numerator is the total number of ways that one can order γ species. The three factorials in the denominator are, respectively, the total number of ways that (i) the number of shared species can be ordered, (ii) the species that are unique to assemblage 1 can be ordered, (iii) the species that are unique to assemblage 2 can be ordered (see Appendix 2.1). In Fig. 2.1, we give two examples of such orderings when γ = 11, Δα = 1 and β = 4.5. Note however, that assemblage pairs for which (2β – Δα) is an odd number are not possible given the inter-dependencies of these parameters. For each possible combination of β and Δα, we randomly selected 1000 pairs of assemblages using code programmed in MATLAB version 7.1 (see Appendices 2.1 – 2.4). For each of these randomly selected assemblages, we calculated FD based on several different dendrogram construction methods. FD was considered insensitive to methodological choices for a particular pair of assemblages if the assemblage with the higher FD was the same for all construction methods; FD was otherwise considered sensitive. We then calculated the proportion of the 1000 random iterations that were sensitive. We refer to this proportion as the probability of sensitivity. When the probability of sensitivity is high for a particular combination of Δα and β, it is very likely that the conclusions drawn from an FD analysis in this context will be dependent on methodological choices, rather than on the data alone. In order to calculate FD, two methodological choices must be made. First, a distance (or resemblance) measure must be chosen. Distance measures quantify the difference between two entities based on their characteristics (e.g. species based on their functional traits). There are a large number of resemblance measures from which to choose (Jackson et al. 1989; Legendre and Legendre 1998). We used three distance measures: Euclidean distance as suggested by Holmes et al. (1979); cosine distance; and, Gower’s distance as it allows mixed and missing data types
24
Figure 2.1 –Measuring the sensitivity of FD in a hypothetical eleven species community. The procedure consists of randomly dividing the community into two assemblages, noting how FD orders the two assemblages and assessing how sensitivity effects this order (i.e. did assemblage 1 always have higher FD given differences in methodology?). Each species is represented as a letter and the assemblages are represented as overlapping rectangles that contain the letters associated with their component species. The first set of rectangles represents one random division of the community into two assemblages. A new random division can be obtained by leaving the overlapping rectangles fixed and simply permuting the order of the species. The second set of overlapping rectangles gives an example of such a permutation.
(Gower 1971; Podani 1999; Podani and Schmera 2006, 2007). For Euclidean distance, we standardized all trait matrices so that all traits have a mean = 0 and variance =1 (i.e. z-scores; Holmes et al. 1979, Gaston and Petchey 2002). We used cosine distance because it more accurately reflects proportional changes in traits whereas the Euclidean distance emphasizes absolute differences. For the Patagonian forb and Arctic vegetation datasets we used only Gower’s distance because these datasets contained missing values and mixed data types; the Euclidean and cosine distances were not appropriate for such datasets (e.g. Podani and Schmera 2006). Second, a clustering algorithm must be chosen. We used three clustering algorithms in this analysis: 1) unweighted pair group method with arithmetic mean (UPGMA); 2) single linkage (i.e. nearest neighbor); and, 3) complete linkage (i.e. maximum or farthest neighbor). These algorithms represent natural endpoints across a methodological continuum of dendrogram construction methods, where single linkage lies on one end, complete linkage on the other and UPGMA lies somewhere in the middle (Podani and Schmera 2006). We considered several different collections of construction methods because the sensitivity of FD is defined in terms of a particular set of construction methods. For cases where multiple comparisons could be made (e.g. several distance measures), we calculated four separate
25
probabilities of sensitivity: 1) sensitivity with respect to all nine construction methods; 2) sensitivity with respect to the three distance measures with UPGMA clustering (i.e. clustering algorithm is held constant); and, 3) sensitivity with respect to the three clustering algorithms with Gower’s distance measure (i.e. distance measure is held constant). In cases where data were deficient and only Gower’s distance could be used, only overall probabilities of sensitivity were calculated. The sensitivity when the clustering algorithm is held constant could be calculated, but these results would be identical to the overall values. Finally, we also calculated probabilities of sensitivity holding other distance measures and clustering algorithms constant and consider pairs of assemblages that do not contain all of the species in the complete datasets. However, we do not present these additional results because they do not alter any of the conclusions. There exists an ongoing debate regarding a standard procedure for calculating FD. Petchey and Gaston (2002) based their measure of FD on a dendrogram derived from a dataset that included all species that were of interest (i.e. the entire community). For an assemblage that does not contain all of the species in the entire community, FD is measured as the total branch length of the dendrogram minus the branch lengths of the species that are not included in the assemblage (see Petchey and Gaston 2002, 2007 for more details). We refer to this approach as the PetcheyGaston (PG) method. Alternatively, Podani and Schmera (2006) suggested that FD should be calculated as the total branch length of a dendrogram that is unique to each assemblage, i.e. recalculated from the reduced dataset. We refer to this measure as the PS (Podani-Schmera) method. As this debate remains unresolved, we tested whether FD was sensitive using both methods. To calculate FD using the PG method, we calculated a species-by-branch matrix and a vector of branch lengths for the complete community using the code of Petchey and Gaston (2002) for the R programming language. We then used this code to calculate FD using the PG method for each assemblage (see Petchey and Gaston 2002 for more details). We repeated this approach for each of the nine construction methods (i.e. three distance measures for each of the three clustering algorthms). To calculate dendrograms using the PS method, we calculated unique dendrograms for all assemblages and construction methods. We used MatLAB (v.7.1) to calculate the sum of dendrogram lengths for each assemblage and construction method. To display all of these results, we constructed image plots with the R programming language. Image plots can be used to show how a variable changes over a two-dimensional grid. The
26
shading of each square on the grid represents the value of the variable at that grid location. In this case, the variable of interest is the probability of sensitivity and the grid is defined by beta diversity, β, and the difference in alpha diversity, Δα, between the two assemblages. However, only certain combinations of β and Δα are possible. For example, for an eleven species community, it is not possible to create two assemblages such that β = 6, Δα = 3 and all of the species are in at least one of the two assemblages. Therefore, for identification purposes, these impossible grid locations are plotted in white whereas all other levels of sensitivity are some shade of grey. Higher levels of sensitivity are represented by darker shades of grey. This results in a checkerboard pattern. However, it is important to keep in mind that the checkerboard pattern is solely an artifact of the impossibility of certain combinations of β and Δα.
Results We identified numerous cases for which FD had a high probability of sensitivity across all communities; that is, it is easy to find cases for which conclusions derived from FD analyses will be driven primarily by methodological choices. In the worst-case scenario, FD sensitivity reached probabilities of 0.976 using the PS method and 0.594 using the PG method (Table 1). Variation in the probabilities of sensitivity was largely driven by variation in alpha diversity, with the highest probabilities of sensitivity found when assemblages were similar in alpha diversity (Figures 2.2 & 2.3). In every case where the probability of sensitivity was zero, FD was larger for the assemblage with more species; this result indicates that FD and alpha diversity lead to identical conclusions about the diversity of assemblages in these cases. Therefore, the hypothesis presented here concerning the relationship between alpha diversity and probability of sensitivity is consistent with these results. Contrary to hypotheses presented here, there were no consistent patterns in the relationship between beta diversity and probability of sensitivity (Figures 2.2 & 2.3). Decisions about distance measures were more important than decisions about clustering algorithms. For example, when UPGMA clustering was kept constant and only distance measures were compared, FD was more sensitive than when Gower’s distance was held constant and clustering methods were compared (Figures 2.2 & 2.3). These results were not altered by the distance measure held constant (e.g. Euclidean, cosine or Gower’s) or by the clustering
27
algorithm held constant (e.g. UPGMA, single linkage and complete linkage), and so we only present the results for holding constant Gower’s distance and UPGMA respectively (Figure 2.2).
Overall PS
Overall PG
Gower PS
UPGMA PS 15
1.0
Difference in alpha diversity
10 5
A)
0.5 3.5
6
8.5 11
0.5 3.5
6
8.5 11
0.5 3.5
6
8.5 11
0.5 3.5
6
8.5 11
0.8
0
0.6 8 6 4 2
B)
0.5
2.5
4.5
6.5
0.5
2.5
4.5
6.5
0.5
2.5
4.5
6.5
0.5
2.5
4.5
0
0.2
6
0.0
6.5
4 2
C)
0.4
Probability of Sensitivity
0 1.5
3.5
5.5
1.5
3.5
5.5
1.5
3.5
5.5
1.5
3.5
5.5
Beta diversity Figure 2.2 - The effect of alpha and beta diversity on the probability of FD sensitivity for three communities crossed with four sets of construction methods. Darker shading represents a higher probability of sensitivity. The communities are A) Insectivorous birds (Holmes et al. 1979), B) Intertidal fish (Munoz and Ojeda 1997) and C) Predatory vertebrates (Jaksic and Medel 1990). Each column is for a different set of construction methods. For the first and second columns, overall PS and PG, all nine methods of dendrogram construction (three distance measures times three clustering algorithms) were used with the PS and PG methods respectively. For the third column, all three clustering algorithms were used with Gower’s distance and the PS method. For the fourth column, all three distance measures were used with UPGMA and the PS method.
28
Difference in alpha diversity
Overall PS
Overall PG 20
1.0
10
0.8 0.6
0
D)
2
7
12
2
7
12
0.4 30
0.2 20 10 0
E)
1
6
11
16
1
6
11
16
0.0 Probability of Sensitivity
Beta diversity Figure 2.3 –The effect of alpha and beta diversity on the probability of FD sensitivity for two communities crossed with two methods of FD calculation. Darker shading represents a higher probability of sensitivity. The communities are D) Patagonian forbs (Golluscio and Sala 1993) and E) Artic vegetation (Chapin et al. 1996). For these communities, only Gower’s distance could be calculated and so only three construction methods could be compared, corresponding to the three clustering algorithms. Each column is for a different method of FD calculation. The first is for the PS method and the second is for the PG method.
29
There are some additional trends worth mentioning. The PG method of FD calculation led to lower probabilities of sensitivity than the PS method in all cases (Figures 2.2 & 2.3, Table 2.1). Also, where greater numbers of dendrogram construction methods are compared, the probabilities of sensitivity increase. For example, compare the overall probabilities (nine construction methods) with the probabilities obtained by holding the clustering method at UPGMA (three construction methods) (Figure 2.2). This difference makes intuitive sense because as one considers more construction methods, it becomes more likely to find a method that leads to different conclusions regarding the ranking of the assemblages in terms of FD. Table 2.1 – The maximum probability of FD sensitivity for five communities previously used to examine FD (Petchey and Gaston 2007; Podani and Schmera 2006). The number of species and the number of assemblage pairs tested are also shown. Data sources are as follows: A) Holmes et al. (1979), B) Munoz and Ojeda (1997), C) Jaksic and Medel (1990), D) Golluscio and Sala (1993), and E) Chapin et al. (1996).
Maximum
Maximum
No. of
Probability of
Probability of
Assemblage
species
Sensitivity: PS
Sensitivity: PG
Pairs Tested
A) Insectivorous birds
22
0.818
0.497
134
B) Intertidal fish
13
0.976
0.366
46
C) Predatory vertebrates
11
0.610
0.594
32
D) Patagonian forbs
24
0.364
0.196
159
E) Arctic vegetation
37
0.244
0.142
370
No. Community
Discussion These results demonstrate that FD is sensitive to choices of distance measure and clustering algorithm in many cases. The major factor contributing to a high probability of sensitivity is low variation in alpha diversity between the assemblages being compared. By contrast, beta diversity between assemblages was a very poor predictor of sensitivity. This did not support the initial hypothesis that lower beta diversity (i.e higher redundancy between traits across assemblages) would lead to a higher probability of sensitivity. The consistency and severity of the results
30
suggest that this sensitivity is not likely to be unique to the examples we present. Indeed, we did not actively search for atypical data to support this position; we merely used the same data that have been used consistently by investigators when evaluating FD (Petchey and Gaston 2002, 2007, Podani and Schmera 2006). If our results are so clear, why did others (e.g. Petchey and Gaston 2007) conclude that decisions regarding methodological choices have only a minor affect on FD, especially given that they used the same data that we use here? There are two possible reasons for this discrepancy. First, to evaluate sensitivity, previous studies have shown that FD calculated using Gower’s distance was strongly collinear with FD calculated using the Euclidean distance across many functional trait matrices (Petchey and Gaston 2007). However, these trait matrices differed widely in number of species. In this analysis, we demonstrate that FD becomes more sensitive as variation in alpha diversity becomes small. Therefore, in the light of this new work, it is not surprising that others have found low sensitivity to methodological choices; in their case, the results strongly suggest that variation in FD was being driven largely by differences in alpha diversity, no matter what methodological choices were made. Second, we compared more distance measures than previously investigated (Podani and Schmera 2007; Petchey and Gaston 2007). We feel this is a more appropriate comparison as there are a large number of distance measures in the multivariate literature deemed to be appropriate. Additionally, when we restricted the analysis to comparing only Gower’s distance and Euclidean distance (with PG dendrogram construction and UPGMA held constant), we found that rates of sensitivity remained high when differences in alpha diversity were low (maximum probability of sensitivity: 0.260 for the bird data, 0.162 for the fish data, and 0.319 for the mammal data). Thus, FD did not provide much additional information in this case, beyond that provided by alpha diversity. The preceding discussion leads to the following important conclusion regarding FD. FD is most sensitive to methodological choices when it genuinely provides new information beyond that provided by alpha diversity. This is because conditions under which FD is sensitive coincide with relatively little variation in alpha diversity between assemblages. Thus in these cases, FD could potentially provide useful information about the differences between the assemblages and ecosystem function. Unfortunately it is precisely in these cases, where FD would genuinely be useful, that it is expected to be highly sensitive to the choice of a distance measure or clustering algorithm. On the other hand, FD is not sensitive to methodological choices, in those cases when
31
it provides very little information beyond that already provided by species richness (alpha diversity). This is because, when FD is insensitive, the results show that alpha diversity is largely redundant with FD no matter what methodological choices are made. Newer approaches to measuring functional richness (e.g. convex hull volume or consensus dendrograms) have been proposed that may reduce the subjectivity of multivariate decisions (Cornwell et al. 2006; Mouchet et al. 2008; Villeger et al. 2008); however, decisions are still required that may alter results (e.g. trait scaling and transformations or what to include in the consensus). Further research into understanding these methodological choices will likely enhance the ability to measure functional richness. Here we wish to raise awareness about the importance of species richness and methodological choices for calculating functional richness and identify cases for which sensitivity is likely to be an issue. What can be done to minimize the impact of sensitivity? One simple approach could be to analyze data from ecological communities using several different construction methods to ensure that sensitivity is not an issue. However, if sensitivity is an issue, a decision must be made. The results suggest potential approaches for reducing the probability of sensitivity. First, we found that probabilities of sensitivity were systematically lower for the PG method of FD calculation than for the PS method. Therefore, one might be tempted to recommend the PG method for general use. There is an important issue with this recommendation however. The PG method assumes that the entire community is known whereas the PS method does not. In a recent paper (Walker et al. 2008), they emphasized the importance of assuming that there may be species in the community that are undiscovered or undetected in the study area when estimating FD from field data. In some cases, this might not be a problem. For example, Barnett et al. (2007) have recently published a list of species to be used in studies of FD in zooplankton communities. However, in the vast majority of cases, there will typically be a high degree of uncertainty about the composition of the entire community. The PG method does not provide the same estimate as the PS method for a subset of the community. Given that the PS method provides the correct dendrogram length for that particular subset, as it is based on a distance matrix constructed from this subset, such differences between the methods remain a concern. Therefore, even though the PS method is more sensitive than the PG method, the PS method is recommended for general use and the PG method when the species list for the entire community is known. Second, we found that FD is much more sensitive to the choice of a distance measure than to the choice of a
32
clustering algorithm. Therefore, one might be tempted to simply adopt a particular distance measure as a standard. However, FD is not completely sensitive to the choice of clustering algorithm (e.g. range in maximum probability of sensitivity across communities: 0.137 to 0.260 for PG method and 0.248 to 0.364 for PS Method). Furthermore, the choice of a distance measure must be made very carefully. It is unlikely that a single distance measure can be found that is justifiable in all situations; indeed, the history of multivariate statistics teaches us that there is no distance measure that can be uniformly recommended in all cases (Sneath and Sokal 1973; Legendre and Legendre 1998). To calculate functional richness, a method for quantifying inter-specific differences in functional traits is required. However, the flexibility to use more than one trait is often required to understand even simple natural systems (Villeger et al. 2008). Unfortunately in these multivariate situations, complications arise as researchers have to make several key decisions during data analysis (e.g. choice of a distance measure, clustering algorithm, data transformations, scaling). Ideally, these decisions should have minimal impact on scientific conclusions. Here we demonstrate that in the case of the popular index of functional richness, FD, decisions inherent in multivariate analyses can drastically alter conclusions of functional diversity and that sensitivity in FD is highest when alpha diversity is low. These results suggest that in cases where information captured by dendrogram-based functional diversity would be most useful, it is redundant with alpha diversity.
Acknowledgements Funding was provided by NSERC Canada and OGS Scholarships to M.S.P and S.C.W., an NSERC Discovery Grant to D.A.J., and the University of Toronto. We thank D.A.R. Drake, J. Podani, D. Schmera, O. Petchey, and anonymous reviewers for comments on early drafts of this paper.
33
References Barnett, A.J., K. Finlay and B.E. Beisner. 2007. Functional diversity of crustacean zooplankton communities: towards a trait-based classification. Freshwater Biology 52: 796-813. Chapin, F. S. I., M.S. Bret-Harte, S.E. Hobbie, and Z. Hailan. 1996. Plant functional types as predictors of transient responses of arctic vegetation to global change. Journal of Vegetation Science 7: 347-358. Cornwell, W. K., D. W. Schwilk, and D. D. Ackerly. 2006. A trait-based test for habitat filtering: convex hull volume. Ecology 87:1465–1471. Golluscio, R. A. and O. E. Sala. 1993. Plant functional types and ecological strategies in Patagonian forbs. Journal of Vegetation Science 4: 839-846. Gower, J. C. 1971. A general coefficient of similarity and some of its properties. Biometrics 27: 857–874. Holmes, R. T., R. E. J. Bonney and S. W. Pacala. 1979. Guild structure of the Hubbard Brook bird community: a multivariate approach. Ecology 60: 512-520. Jackson, D. A., K. M. Somers and H. H. Harvey. 1989. Similarity coefficients: measures of cooccurrence and association or simply measures of occurrence? American Naturalist 133:436-453. Jaksic, F. M., and R.G. Medel. 1990. Objective recognition of guilds: testing for statistically significant species clusters. Oecologia 82: 87-92. Legendre, P., and L. Legendre. 1998. Numerical Ecology, 2nd ed. Elsevier B.V. Amsterdam, The Netherlands. Lande, R. 1996. Statistics and partitioning of species diversity, and similarity among multiple communities. Oikos 76: 5-13. Loreau, M., S. Naeem, P. Inchausti, J. Bengtsson, J.P. Grime, A. Hector, D. U. Hooper, M. A. Huston, D. Raffaelli, B. Schmid, D. Tilman, and D. A. Wardle. 2001. Biodiversity and ecosystem functioning: Current knowledge and future challenges. Science 294: 804-808.
34
Mason, N.W.H., D. Mouillot, W. G. Lee, and J. B. Wilson. 2005. Functional richness, functional evenness and functional divergence: the primary components of functional diversity. Oikos 111: 112-118. Mouchet, M., F. Guilhaumon, S. Villeger, N. W. H. Mason, J. A. Tomasini, and D. Mouillot. 2008. Towards a consensus for calculating dendrogram-based functional diversity indices. Oikos 117: 794-800. Munoz, A. A., and F. P. Ojeda. 1997. Feeding guild structure of a rocky intertidal fish assemblage in central Chile. Environmental Biology of Fishes 49: 471-479. Petchey, O. L., and K. J. Gaston. 2002. Functional diversity (FD), species richness and community composition. Ecology Letters 5:402-411. Petchey, O. L., and K. J. Gaston. 2006. Functional diversity: back to basics and looking forward. Ecology Letters 9: 741-758. Petchey O. L., and K. J. Gaston. 2007. Dendrograms and measuring functional diversity. Oikos 161: 1422-1426. Petchey, O.L., A. Hector, and K. J. Gaston. 2004. How do different measures of functional diversity perform. Ecology 85: 847-857. Podani J., and D. Schmera. 2006. On dendrogram-based measures of functional diversity. Oikos 115: 179-185. Podani J., and D. Schmera. 2007. How should a dendrogram-based measure of functional diversity function? A rejoinder to Petchey and Gaston. Oikos 116: 1427-1430. Sneath, P.H.A., and R. R. Sokal. 1973. Numerical Taxonomy. The principles and practice of numerical classification. W.H. Freeman and Company, San Francisco, United States. Walker, S.C, M. S. Poos, and D. A. Jackson. 2008. Functional rarefaction: estimating functional diversity from field data. Oikos 117: 286-296.
35
Villeger S., N.W.H. Mason, and D. Mouillot. 2008. New multidimensional functional diversity indices for a multifaceted framework in functional ecology. Ecology 89: 2290-2301.
36
Appendices Appendix 2.1 – Derivation of Equation 1. Definitions Let γ be the total number of species in both assemblages. Let α1 be the number of species in assemblage one. Let α2 be the number of species in assemblage two. Let Δα = α1 – α2. Let <α> be the average of α1 and α2. Let β = γ – <α>. Let ν be the number of shared species (i.e., the number of species in both assemblages). Assumption All of the γ species are in assemblage 1, assemblage 2, or both. Main result (Eq. 1) For a given γ, Δα and β, the total number of unique assemblage pairs is [γ!] / [(γ – 2β)!(0.5(2β + Δα))!(0.5(2β – Δα))!]. Deriving this result is much easier once there are three simpler results. Simpler result 1 X! is the total number of ways to order X objects. Simpler result 2 (γ – 2β) is the number of shared species, ν.
37
Derivation of simpler result 2 Given that all species are in at least one of the two assemblages (see assumption above), it follows from the inclusion-exclusion principle that γ = α1 + α2 – ν. It must also be that β = γ – 0.5(α1 + α2) from the definitions of β and <α>. Eliminating the alphas and solving for ν leads to simpler result 2. Corollary of simpler result 2 2β = α1 + α2 – 2ν is the total number of unshared species (i.e., species that are only in either assemblage 1 or assemblage 2 but not in both). Simpler result 3 0.5(2β + Δα) and 0.5(2β – Δα) are the numbers of unshared species in assemblages 1 and 2, respectively. Derivation of simpler result 3 It is only necessary to derive the first claim; the second follows immediately after the first. Note that one can write Δα = (α1 – ν) – (α2 – ν), from the definition of Δα, and 2β = (α1 – ν) + (α2 – ν), from the corollary of simpler result 2. Eliminating α2 – ν and solving for α1 – ν one can obtain α1 – ν = 0.5 (Δα + 2β). Simpler result 3 follows once it is recognized that α1 – ν is the number of unshared species in assemblage 1. Derivation of main result It follows from simpler result 1 that the numerator in Eq. 1 is the number of ways that γ species can be ordered. From Fig. 1, one can see that the numerator gives the total number of ways that iterations from the simulations could occur. However, many of these ways result in identical assemblage pairs. This is because the order that species are listed in Fig. 1 determines a particular assemblage pair, only insofar as it determines whether each species is in assemblage 1 only, assemblage 2 only or in both assemblages 1 and 2. Therefore one must divide by the total number of ways to order the species within each of these three groups. It follows from simpler results 1 through 3 that the denominator in Eq. 1 gives this number.
38
Appendix 2.2 – MatLAB Code for testing sensitivity of FD % X is the trait matrix (S by T) % S is species richness % T is # of traits % robustness measures % 1. overall % 2. euc % 3. cos % 4. gower % 5. UPGMA % 6. single % 7. complete S = length(X(:,1)); gamma = 10; iters = 10; n_measures = 9; % number of different ways to calculate FD G = zeros(iters,n_measures); output = -999.*ones((gamma+1),(gamma+1),iters,n_measures); % sign matrix nonrobust_probs = -999.*ones((gamma+1),(gamma+1),7); % use when all dist measures %nonrobust_probs = -999.*ones((gamma+1),(gamma+1),1); % use when only gowers (remove % and put one on the above) for B = 4:gamma for A = max(B,(gamma-B)):gamma delta = A - B overlap = A + B - gamma unsharedA = A - overlap unsharedB = B - overlap startA = 1; endA = unsharedA; startB = unsharedA + 1; endB = unsharedA + unsharedB; startShare = endB + 1; endShare = startShare + overlap - 1; Alist = zeros(iters,A); Blist = zeros(iters,B); Slist = zeros(iters,S);
39
FDA = zeros(iters,n_measures); FDB = zeros(iters,n_measures); if overlap == 0 for i = 1:iters currlist = randperm(S); Slist(i,:) = currlist; Alist(i,:) = [currlist(startA:endA)]; Blist(i,:) = [currlist(startB:endB)]; FDA(i,1) = FD(X(Alist(i,:),:),'euclidean','average'); FDB(i,1) = FD(X(Blist(i,:),:),'euclidean','average'); FDA(i,2) = FD(X(Alist(i,:),:),'euclidean','single'); FDB(i,2) = FD(X(Blist(i,:),:),'euclidean','single'); FDA(i,3) = FD(X(Alist(i,:),:),'euclidean','complete'); FDB(i,3) = FD(X(Blist(i,:),:),'euclidean','complete'); FDA(i,4) = FD(X(Alist(i,:),:),'cosine','average'); FDB(i,4) = FD(X(Blist(i,:),:),'cosine','average'); FDA(i,5) = FD(X(Alist(i,:),:),'cosine','single'); FDB(i,5) = FD(X(Blist(i,:),:),'cosine','single'); FDA(i,6) = FD(X(Alist(i,:),:),'cosine','complete'); FDB(i,6) = FD(X(Blist(i,:),:),'cosine','complete'); FDA(i,7) = FD(X(Alist(i,:),:),'gowers','average'); % change needed FDB(i,7) = FD(X(Blist(i,:),:),'gowers','average'); % change needed FDA(i,8) = FD(X(Alist(i,:),:),'gowers','single'); % change needed FDB(i,8) = FD(X(Blist(i,:),:),'gowers','single'); % change needed FDA(i,9) = FD(X(Alist(i,:),:),'gowers','complete'); % change needed FDB(i,9) = FD(X(Blist(i,:),:),'gowers','complete'); % change needed end else for i = 1:iters currlist = randperm(S); Slist(i,:) = currlist; Alist(i,:) = [currlist(startA:endA),currlist(startShare:endShare)]; Blist(i,:) = [currlist(startB:endB),currlist(startShare:endShare)];
40
FDA(i,1) = FD(X(Alist(i,:),:),'euclidean','average'); FDB(i,1) = FD(X(Blist(i,:),:),'euclidean','average'); FDA(i,2) = FD(X(Alist(i,:),:),'euclidean','single'); FDB(i,2) = FD(X(Blist(i,:),:),'euclidean','single'); FDA(i,3) = FD(X(Alist(i,:),:),'euclidean','complete'); FDB(i,3) = FD(X(Blist(i,:),:),'euclidean','complete'); FDA(i,4) = FD(X(Alist(i,:),:),'cosine','average'); FDB(i,4) = FD(X(Blist(i,:),:),'cosine','average'); FDA(i,5) = FD(X(Alist(i,:),:),'cosine','single'); FDB(i,5) = FD(X(Blist(i,:),:),'cosine','single'); FDA(i,6) = FD(X(Alist(i,:),:),'cosine','complete'); % change needed FDB(i,6) = FD(X(Blist(i,:),:),'cosine','complete'); % change needed FDA(i,7) = FD(X(Alist(i,:),:),'gowers','average'); % change needed FDB(i,7) = FD(X(Blist(i,:),:),'gowers','average'); % change needed FDA(i,8) = FD(X(Alist(i,:),:),'gowers','single'); % change needed FDB(i,8) = FD(X(Blist(i,:),:),'gowers','single'); % change needed FDA(i,9) = FD(X(Alist(i,:),:),'gowers','complete'); % change needed FDB(i,9) = FD(X(Blist(i,:),:),'gowers','complete'); % change needed end end output(overlap+1,delta+1,:,:) = sign(FDA-FDB); G(:,:) = output(overlap+1,delta+1,:,:); nonrobust_probs(overlap+1,delta+1,1) = sum(abs(sum((G(1:iters,:))'))
41
Appendix 2.3 – MatLAB Code for calculating FD via Podani and Schmera function output = FD(X,distance,cluster) if strcmp(distance,'gowers') Y = gowers(X); else Y = pdist(X,distance); end Z = linkage(Y,cluster); output = sum(sum(branch_lengths(Z)));
42
Appendix 2.4 – MatLAB Code for calculating probabilities of sensitivity S = length(X(:,1)); iters = 1000; n_measures = 9; % number of different ways to calculate FD Amax = (ceil(S/2)-3); Bmax = (floor(S/2)-3); Omax = (S-8); dist.probs = zeros(Bmax,Amax,Omax,7); G = zeros(iters,n_measures); for BBB = 1:Bmax for AAA = BBB:Amax for OOO = 1:(S-5-AAA-BBB) G(:,:) = output(BBB,AAA,OOO,:,:); dist.probs(BBB,AAA,OOO,1) = sum(abs(sum((G(1:iters,:))'))
43
Chapter 3: Addressing the removal of rare species in bioassessments with other choices in multivariate analyses
Abstract The inclusion or exclusion of rare species in the bioassessment of aquatic communities has been greatly debated. Researchers may include rare species in bioassessments as they are likely better indicators of ecosystem stress than more common species (i.e. the biological argument). Alternatively, researchers may exclude rare species due to the potential influence on statistical analyses (i.e. the statistical argument). As this debate remains unresolved, the objective of this study was to determine the impacts of removing rare species in multivariate bioassessments. These approaches were tested independently using multivariate comparisons of fishes from a thoroughly sampled system. The biological argument was assessed using sites-level vector residuals across treatments of species removal and demonstrated that the removal of rare species had important site-level implications relative to full dataset, including up to a nine-fold change in multivariate vector residuals at sites where single species were removed. The statistical argument was assessed using variation partitioning of multivariate decisions such as ordination method, distance measure and removal of rare species, and found that the removal of rare species demonstrated similar levels of multivariate variation (e.g. 24.8% variation) as other choices inherent in multivariate bioassessments, such as the choice of ordination technique (26% variation) and similarity measure (11%). This study demonstrates that contrary to the common held practice of removing rare species in multivariate bioassessments, that the removal of rare species may be less important than previously thought, whereas other multivariate decisions may be at least as equally important. Better justification for the removal of rare species, along with all decisions in multivariate analyses, is needed to ensure bioassessments are developed in a rigorous manner. Keywords: Multivariate analysis; ordination; similarity measures; rare species; community ecology; bioassessment; procrustes analysis.
44
Introduction The use of multivariate analyses has become an important tool in the biological assessment of aquatic communities (Norris 1995; Wright et al. 2000). Several national bioassessment programs are based on multivariate measures, notably those in the UK (e.g. RIPVACS; Wright et al. 2000) and Australia (AUSRIVAS; Simpson and Norris 2000; Metzeling et al. 2006); and the use multivariate analyses has become widespread elsewhere (Reynoldson et al. 2001; Joy and Death 2002; Bailey et al. 2004). In multivariate analyses, researchers order entities in the data (e.g. species or observations) on the basis of the similarity of their characteristics (e.g. observations or species) (Wartenberg et al. 1987). The goal of such analyses is to determine the basis for the order of entities; for example, differences in species abundance or occurrence between impacted sites versus sites with little or no impact (Barbour et al. 1999; Wright et al. 2000). From these orderings, one may be able to infer causative relationships between species and their environment so that site-level impacts can be identified and mitigated (Hawkins et al. 2000). The application of multivariate analyses to bioassessments of aquatic communities has been riddled with controversy. Debates in bioassessment literature include the use of multi-metric versus multivariate approaches (Hannaford and Resh 1995; Wallace et al. 1996; Bowman and Somers 2006); the amount of taxonomic resolution needed to determine site level impacts (Somers et al. 1998; Hewlett 2000); the use of rapid assessment methods (Hannaford and Resh 1995; Somers et al. 1998); and issues with using reference conditions (Norris and Thoms 1999; Bowman and Somers 2005; Bailey et al. 2008). In particular, the debate regarding the treatment of rare species has received much attention (e.g. Faith and Norris 1989; Norris 1995; Cao et al. 1998; Cao and Williams 1999; Marchant 1999; 2002). On one hand, researchers remove rare species with the perceived notion that they may add noise to multivariate analyses and provide little additional information beyond more common species (Gauch 1982; Clarke and Green 1988; Marchant 1990; 2006; McCune and Grace 2002; Paukert and Wittig 2002). On the other hand, researchers retain rare species in multivariate analyses because they may be better indicators of ecosystem stress than common species (Faith and Norris 1988; Cao et al. 1998; 2001), as they may be more sensitive to the stressor(s). In either case, the debate regarding the treatment of rare
45
species has remained unresolved and researchers need to be aware of the impact of their decision of how to treat rare species (among others). There are many difficulties in attempting to resolve the debate regarding the treatment of rare species in bioassessments. For example, most multivariate approaches require several decisions beyond the removal of rare species, and these decisions may reduce insight into the effect of rare species on resultant analyses. Researchers using multivariate methods typically must choose a type of similarity coefficient and ordination technique, where such choices have been shown to significantly alter results (Podani 2000; Podani and Schmera 2006; Hirst and Jackson 2007; Poos et al. 2009). As such, the resolution to the debate regarding the impact of rare species cannot proceed until the effect of removing rare species is placed into a context comparable to other decisions inherent in multivariate bioassessments. Unfortunately, relatively little effort has gone into comparing these methods of bioassessment (Norris 1995; Marchant et al. 2006), and few studies have viewed decisions in analyses in a holistic manner (e.g. how do all of the decisions inherent in multivariate bioassessments affect results?). There are two general arguments for the inclusion or exclusion of rare species in multivariate bioassessments. The first argument for removing rare species from bioassessments is that rare species provide limited interpretative value (Marchant 1999). Proponents of this argument suggest that rare species may simply reflect stochastic sampling effects and, as such, add noise to the statistical solution (Gauch 1982; Clarke and Green 1988; Bailey et al. 2004). This argument is referred to as the “statistical argument” for excluding rare species. Support for this argument has come from results from multivariate methods could be driven by the inclusion of rare species alone (Webb et al. 1967, Austin and Greig-Smith 1968, Day et al. 1971, Orloci and Mukkattu 1973). To some degree this argument has been examined in the literature with analyses of data standardizations (Jackson 1993; Cao et. al 1999), similarity coefficients (Jackson et al. 1989), ordination method (Marchant 1990), or their combinations (e.g. data standardization and similarity coefficients; Jackson 1993; Hirst and Jackson 2007; taxonomic resolution and rarity, Arscott et al. 2006). Unfortunately, a quantitative evaluation of the role of rare species in community assessments is largely absent, including a more holistic evaluation that answers the practical question of how important rare species are relative to other decisions in multivariate analyses (but see Faith and Norris 1988). In this context, the statistical argument can be tested as a hypothesis, with the prediction that multivariate analyses should show greater variation among
46
one another where rare species are removed than multivariate analyses where other decisions have been altered (e.g. similarity measure, ordination method) – specifically, does the inclusion or exclusion of rare species lead to greater changes in multivariate analyses than those arising due to other decisions? If not, then the inclusion of rare species may be warranted as this decision contributes similar or lower amounts of variation versus other decisions typical in a bioassessment. The second argument for the inclusion or exclusion of rare species is that such species should be included in multivariate analyses because they are better indicators of ecosystem stress than are common species (Cao et al. 1999; 2001a), i.e. common species tend to have broad ranges of tolerance to many conditions and, therefore, are not good indicators. This argument is referred to as the “biological argument” for the inclusion of rare species. Support for this argument has come from empirical studies that note the importance of including rare species for conservation issues (Margules 1986; DeVelice et al. 1988; Norris and Hawkins 2000). Proponents of this argument suggest that exclusion of rare species may lead to an underestimation of differences between impacted and un-impacted sites (Cao et al. 1999; 2001a). As one of the main goals of a bioassessment is to determine site-level impacts (Barbour et al. 1999; Wright et al. 2000), this argument assumes that sites with rare species represent the strongest signals, such as decreases in species diversity or changes in community composition (Cao et al. 1998). Therefore, the biological argument can be tested also as a hypothesis, with the prediction that once rare species are removed from multivariate analyses, sites which were chosen for the removal of rare species should be more affected across multivariate analyses than sites that were not (i.e. when rare species are removed from the analysis, the site-level signal will change in greater proportion at sites with rare species than those sites without). If this result is not found, it would indicate that the exclusion of rare species is warranted as they do not provide information beyond that captured by more common species. Given these hypotheses, the objective of this study is to determine the biological and statistical effect of removing rare species relative to other methodological decisions inherent in multivariate analyses (e.g. choice of ordination method and distance measure). For this evaluation, the Sydenham River is used as a model system because it has the highest diversity of aquatic fauna in Canada, as well as the highest number of species at risk in Canada (Staton et al. 2003). Also, the Sydenham River has undergone detailed sampling (Poos et al. 2007; 2008),
47
which provides a high-quality dataset. As there are many types of rarity (Raboniwitz et al. 1986; Gaston 1994), for clarity rare species are defined as those which occur infrequently (i.e. at few locations or low prevalence, e.g. <1%, 5% and 10% occurrence). As well, species that are rare due to declines in population sizes or number of locations and have a conservation designation (e.g. endangered, threatened, special concern) are referred to as species at risk.
Methods Fishes were collected from the Sydenham River using the Ontario Stream Assessment Protocol (OMNR 2007). The Ontario Stream Assessment Protocol is a typical bioassessment protocol for monitoring impacts to aquatic systems (e.g. Barbour et al. 1999; Wright et al. 2000) and has been used to monitor the changes in riverine communities (Stanfield and Jones 1998; Poos et al. 2008). Fishes were used as model organisms rather than benthic macroinvertebrates as fish are relatively easy to identify, and enumerate and have been used extensively in bioassessments (e.g. Fausch et al. 1990; Joy and Death 2002; Boys and Thoms 2006; Kennard et al. 2006; Mugodo et al. 2006); however, this approach is equally well suited to macroinvertebrates or any other group of species. Further, most bioassessments using benthic macroinvertebrates (Marchant 1990; Marchant et al. 1997; Cao et al. 1998) are sensitive to sampling method and taxonomic resolution (Arscott et al. 2006; Nichols and Norris 2006). Fishes were sampled using a variety of approaches (see Poos et al. 2007 for sampling protocol); however, for this analysis only electrofishing data were used as it is the most commonly used bio-monitoring protocol and regarded as the most effective gear type for sampling stream-fish assemblages (Bohlin et al. 1989; Reynolds 1996). As sample representativeness may be an issue (Cao et al. 2001a; 2002), electrofishing was assessed relative to other methods and determined to be the most effective method for evaluating fish species at risk (Poos et al. 2007). Species were collected at 50 sites in 2002 and 25 additional sites in 2003. Sampling sites were chosen at random across the entire watershed, with the exception of non-wadeable sections of the river in the lower portions of the watershed which were not sampled due to constraints with using wadeable sampling gear.
48
Evaluating Decisions in Multivariate Bioassessments Prior to analysis, four treatments of the removal of species were applied to the species matrix (Table 3.1). Data transformation and standardization have been previously shown to influence multivariate analyses (Austin and Greg-Smith 1968; Jackson et al. 1993; Cao et al. 1999) and, not wanting to provide potential bias, data were reduced to presence/absence as it removes one additional source of variation from the comparisons and provides focus on specific comparisons of multivariate methods. Traditionally, researchers decide what characteristics define a rare species within a sample (Faith and Norris 1989; McCune and Grace 2002). Some researchers suggest eliminating species that occur at single sites because of the inflated correlations created by attempting to relate potentially random features at that site to its lone occurrence (Legendre and Legendre 1998). Others suggest removing species that occur at less than five percent (McGardial et al. 2000) or ten percent of sites (Marchant 1990; McCune and Grace 2002) or at even higher thresholds (Boulton et al. 1992; Marchant et al. 1997). The following treatments of removing rare species were used: all species included, single-occurrence species removed; species found at less than five percent of sites removed; and, species found in less than 10% of sites removed. These criteria removed 0, 2, 8, and 21 species respectively of the 67 species dataset.
49 Table 3.1 – Summary of ordination techniques, similarity coefficients and exclusion of rarely sampled species being compared. Abbreviations are indicated in parentheses and used in subsequent figures and tables. All four approaches described in the “Exclusion of Rarely Sampled Species” were used in each of the four “Similarity Coefficient” combinations with both PCoA and NMDS. As Correspondence Analysis has the implicit chi-squared distance measure, only the four approaches used in the “Exclusion of Rarely Sampled Species” were included in that set of analyses.
Ordination Technique
Similarity Coefficient
Principal Co-ordinates • Jaccard’s (J) Analysis (PCoA) • Phi (Φ) • Russell and Rao (RR) • Simple Matching (SM) * Non-metric multidimensional scaling (NMDS) Correspondence analysis (CA)
• • • •
Jaccard’s (J) Phi (Φ) Russell and Rao (RR) Simple Matching (SM)
• X2 distance (no other choice of similarity coefficient)
Exclusion of Rarely Sampled Species • No species removed (All) • Single occurrences removed (M1) • < 5% occurrences removed (M5) • < 10% occurrences removed (M10) • Same as above.
• Same as above.
Note: * Principal Co-ordinates Analysis (PCoA) using a simple matching similarity is identical to a Principle Components Analysis (PCA) using a correlation matrix (Gower 1966).
50
Similarity coefficients were calculated from each of the four matrices of rarely sampled species (Table 3.1). Several dozen similarity coefficients have been developed for use with presence/absence data (Gower 1966; Gower and Legendre 1986; Legendre and Legendre 1998; Podani 2000). The choice of similarity coefficient has been largely subjective and is often based on tradition or on a posteriori criteria without clear justification (Jackson et al. 1989; Krebs 1998). As different similarity coefficients emphasize different aspects of the relation between observations, the exclusion of rare species may alter species relationships and subsequent results of analyses such as ordinations (McGargal et al. 2000). Treatments of similarity coefficient included: Jaccard’s, phi (Φ), Russell and Rao, and simple matching coefficients, and were chosen because they represent standard examples amongst the continuum of similarity coefficients. Jaccard’s similarity does not consider joint absences; the phi coefficient is the correlation coefficient for binary data; and, Russell and Rao and simple matching are variations that consider joint absences (Jackson et al. 1989; Legendre and Legendre 1998; Podani 2000). All similarity coefficients were transformed into metric distances having Euclidean properties for subsequent analysis (Jackson et al. 1989; Legendre and Legendre 1998). Three types of ordination technique were compared for each combination of treatments excluding rare species and using different similarity coefficients: principal co-ordinates analysis (PCoA), non-metric multidimensional scaling (NMDS), and correspondence analysis (CA). As PCoA measured from simple matching similarity is identical to PCA measured using a correlation coefficient (Gower 1966), the results can be used to interpret both ordination types. These ordination techniques were chosen because they represent typical multivariate methods used by the majority of biologists (Legendre and Legendre 1998; Podani 2000; McCune and Grace 2002). In addition, more current approaches (e.g. NMDS, CA) may provide robust alternatives to previous methods (e.g. PCA) where non-linear relationships occur between variables (Cao et al. 2001a). Whereas both NMDS and PCoA allow the user to choose a similarity coefficient, CA do not provide the same option given its inherent resemblance measure (chi-square distance). For NMDS, a random set of 20 starting configurations were used as input configurations, and the solution having the lowest stress was retained. In NMDS, stress was measured as an objective function of a regression analysis where the goodness of between the fitted values and forecasted values was fit using a least square criterion (Legendre and Legendre 1998). A broken-stick model was used to compare the eigenvalues from PCoA and CA to those
51
expected from random relationships. This method has been shown useful in identifying nonrandom patterns of association in multivariate analyses (Jackson 1993; King and Jackson 1999; Peres-Neto et al. 2003). Axes not representing a meaningful contribution of the variation were removed from resultant analyses. All analyses were completed using the R programming language v2.70 plus statistics libraries simba (Jurasinski 2007), vegan (Oksanen et al. 2008) and ecodist (Goslee and Urban 2007).
Assessing the Statistical Argument The statistical argument in multivariate bioassessments was assessed in several ways. First, all variants in ordination method, similarity coefficients and exclusion of rare species were compared using Procrustes analysis (Gower 1971; Jackson 1995). Procrustes analysis is appropriate for comparing separate ordination results because it is an orthogonal rotation that best matches two or more ordinations (Olden et al. 2001; Peres-Neto and Jackson 2001; Paavola et al. 2006). The first three dimensions from each ordination solution were retained for comparisons as they represented the greatest portion of variance explained using a broken-stick model (Legendre and Legendre 1998) and the majority of the ordination methods were best represented by three-dimensional solutions, including the NMDS results. The sum-of-squareddeviations (i.e. m2 statistic) was used as a metric of association, with lower sum-of-squareddeviations representing greater similarity of multivariate configurations (Gower 1971; Jackson 1993; Peres-Neto and Jackson 2000) and was calculated between each pair of three-dimensional ordination solutions to produce a matrix of m2 distances between all 36 exclusion-distanceordination combinations. The resultant 36-by-36 matrix of m2 distances was analyzed using a PCoA to determine the relative effect of each methodological choice. This type of “ordination of ordinations” (see Digby and Kempton 1987; Jackson 1993; Hirst and Jackson 2007) provides a useful characterization of methodological decisions, where larger distances between objects represent more dissimilar associations. A minimum spanning tree was calculated to determine the most similar groups and super-imposed onto the first two axes of the ordination diagram. Partitioning of variation of multivariate data can provide quantitative and objective determination of the influence of methodological choices. For multivariate datasets, partitioning of variation is often thought of in a spatial or temporal context, where the influence of variables can be partitioned across various spatial or temporal scales (e.g. Borcard et al. 1992; Borcard et
52
al. 2004; Dray et al. 2006). Yet, partitioning of variation in multivariate data is also possible through the analysis of residuals across various matrix comparisons (e.g. Rundle and Jackson 1996; Olden et al. 2001; Paavola et al. 2006). For example, the total of among-group variation of removing rare species can be summarized relative to all treatments using the sum of squared deviations from a Procruste’s analysis, and represents the amount of multivariate variation explained. In this study, the variation of all treatments was quantified using a partitioning method of multivariate matrices (see Rundle and Jackson 1996), which separated the sum-of-squareddeviations for within- and among-treatment groups (e.g. removal of rare species, ordination technique, and similarity measure).
Assessing the Biological Argument To assess the biological argument of removing rare species from bioassessments, changes at sites where rare species were removed were evaluated across the various multivariate analyses. Sitelevel differences were calculated for each pair-wise Procrustes analysis using vector residuals from PROtest (Jackson 1995). Vector residuals provide a means of investigating deviations in position of individual samples between two superimposed ordinations (Olden et al. 2001; Paavola et al. 2006), i.e. the degree to which any given observation changes from one ordination to another ordination. The length of the vector residual represents a lack of fit of scores for an individual sample between two ordinations, with low values indicating close agreement between multivariate methods. Vector residuals were separated between sites where rare species were removed with sites where rare species where no species were removed. For example in Figure 3.1, a typical example of site level vector residuals is shown across a comparison of multivariate analyses (e.g. PCoA with Jaccard’s distance) between all data and with species occurring at 5% of sites removed. From this comparison, the effect of the removal of rare species can be assessed as the ratio of site-level vector residuals for sites where rare species were removed (i.e. grey bars; Fig. 3.1) versus the sites where species were not removed (i.e. black bars; Fig. 3.1). Ratios over 1 indicate situations where site-specific differences are more related to the removal of rare species than those that not (i.e. bioassessments may be affected by species removal as sites where rare species occur, change in greater magnitude than sites where rare species do not occur when compared to a full dataset).
53
Vector Residual
0.012
0.006
0 1
75
Ordered Sites Figure 3.1 – Example of rank-ordered, site-level vector residuals of Procrustean multivariate comparison. The length of a vector residual indicates an overall lack of fit for a site between two multivariate analyses. Shown is a comparison of full dataset of Principal Coordinates with Jacaard’s distance and the same dataset where species occurring at 5% of sites were removed. Vectors shown in grey indicate those sites where at least one species was removed, whereas vectors in black indicate sites where no species were removed. The ratio of mean vector residuals between sites where species were removed versus those sites that did not have species removed indicates the distribution of impacts of the removal of rare species across multivariate analyses. Where most vector residuals for sites having species removed are largest, they indicate that these observations (sites) have been changed the most in their position between two ordinations.
54
Results The multivariate analyses used in this study provided generally good representation of the data. Variance explained by the first three axes from all combinations of multivariate analyses ranged from 24.6% (PCoA with Jaccard’s distance and all species included) to 38.4% (PCoA with simple matching similarity and species occurring at less than 10% sites removed). In all cases, the variance explained by each multivariate method increased with the exclusion of more species. For fish species in the Sydenham River, the removal of rare species had similar effects on resultant analyses to those arising from other decisions in multivariate analyses, such as the choice of ordination type or similarity coefficient. Recall that in considering the relative role of the various decisions to be made, the partitioning of variation provides a measure summarizing the relative importance. Variation across multivariate analyses was highest for ordination method (26.15%) across all comparisons, followed by the removal of rare species (24.81%) and similarity measure (10.99%; Table 3.2). These results were also evident by the well-defined clustering of treatments of rarely sampled species (All, M1, M5) in close proximity to one another relative to the clustering of the differences in ordination technique (PCoA, NMDS, CA, PCA) or in similarity coefficients (J, Φ RR, SM) in the ordination of m2 distances, i.e. the comparison of the various ordination results (Fig. 3.2). One clear exception to this result was the removal of 10% of the least occurring species, which showed deviations from the general multivariate groupings (Fig. 3.2a,b), and variation that exceeded most other choices (6.82%; Table 3.2). There was large variation between individual choices across multivariate methods. Whereas CA and PCA demonstrated overall low amounts of variation among analyses (0.19%; 0.039% respectively), there was large variation among analyses based on NMDS (21.60%; Table 3.2). These differences may be influenced by the smaller number of comparisons for correspondence analysis (as a choice of similarity measure is implicit and not selected); however, NMDS also showed an almost seven-fold increase in variation over PCoA, which involved same number of
55
0.3
Axis 2 (20.7%)
PCoA-J M10
All, M1, M5
CA
M10
All, M1, M5 PCoA-RR M10 All, M1, M5 All, M1, M5
M10
0
All M1 M5
M10
PCoA-Φ PCA
All, M1, M5
PCoA-SM
All, M1, M5
M10
NMDS-RR
M10
All, M1, M5
NMDS-J
M10
NMDS-Φ
All M1 M5
NMDS-SM
M10
-0.3 -0.45
0.45
0
Axis 1 (47.3%)
a) 0.15
M10
M10 M5
Axis 3 (15.6%)
M10
0
M5 M1 All
NMDS-SM M10
M10
M1 All
NMDS-RR
PCA PCoA-RR
M5 M5 M1 M1 All All
PCoA-SM
M10 M5 M1 All
M10
PCoA-J
M5 All M1
NMDS-J
PCoA-Φ
M1 All M5
M10
M5,M1,All
NMDS-Φ
M10
-0.3
b)
All,M1,M5
-0.3
0
CA
0.30
Axis 2 (20.7%) Figure 3.2 –Principal Coordinates Analysis (PCoA) of the sum of squares deviations (m2 statistic) comparing the concordance between solutions based on different ordination techniques, similarity coefficients and treatments of excluding rarely sampled species. A minimum-spanning tree was overlaid on Axes 1 and 2 to highlight connections between groups of points. Dashed lines indicate deviations from group membership in cases where clear groupings do not exist (e.g. M10 for Axes 2 and 3). Short forms are continued from Table 3.2.
56
comparisons (Table 3. 2). Decisions, such as, which ordination method to choose, may be as important (or more important) than other choices like the removal of rare species. In fact, the inclusion or exclusion of rare species did not impact the resulting multivariate analyses any more than the choice of ordination method (i.e. subtotals of 24.8% vs 26.2%, respectively) and likely less given that four comparisons regarding species inclusion are used and only three comparisons for ordination method. Levels of variation were similar between the removal of single occurrence species (6.03%) and species occurring at 5% of the sites (5.93%) as they were for using the entire dataset (6.03%; Table 3.2). Finally, the choice of similarity measure showed lower levels of variation in general. Variation in X2 values was lowest (0.19%), followed by simple matching (1.53%), Jaccard (2.03%) and the Phi (2.16%) coefficients (Table 3.2). Table 3.2 –Partitioning of variation in sum of squared deviations of Procrustes analyses (m2 statistic) across various choices in multivariate analyses, including: i) removal of rare species; ii) ordination technique; and, iii) choice of similarity measure. Abbreviations are those noted in Table 1.
Variation component
%
I) Removal of Rare Species All 6.03 M1 6.03 M5 5.93 M10 6.82 Subtotal 24.81 II) Ordination Technique PCoA 4.37 NMDS 21.60 CA 0.19 26.15 Subtotal III) Dissimilarity Measure J 2.03 Phi 2.16 RR 5.08 SM 1.53 2 0.19 X 10.99 Subtotal
57
There was high site-level effect of removing rare species. In most cases, when rare species were removed from the analysis, the effects were driven by differences in sites that contained species that were removed. For example, the site-level residuals were much higher in sites impacted by the removal of species than the full data (Figure 3.3; Appendix 3.1). Recall that site-level vector residuals represent the degree to which any given observation changes from one ordination to another (Olden et al. 2001; Paavola et al. 2006). Therefore, the ratio of site-level residuals between sites impacted by species removal and those sites that did not have species removed provides an indication as to magnitude that rare species may alter site-level assessments. Overall, sites impacted by the removal of rare species changed in multivariate space over nine fold when single-occurrence species were removed from the analysis, relative to the full data set, and over two fold when species having prevalence less than 5% were removed from the analysis. Interestingly, once species that occurred at less than 10% of sites were removed from the analyses, there was virtually no difference between the two categories of sites (1.13 difference; Figure 3.3), and in some cases represented less of an impact (e.g. NMDS-J, NMDS-RR; Appendix 3.1). 15
Procruste's vector residual
12 9 6 3 0 M1
M5
M10
Removal of Missing Species Figure 3.3 – Site-level impact of the removal of rare species. Shown are box and whisker plots of the ratios of Procrustes vector residuals between sites for which rare species were removed and those sites that did not have any species removed. All comparisons were done by comparing site-level Procrustes vector residuals from the full datasets and with the removal of rare species across all similarity coefficients and ordination methods.
58
Discussion The treatment of rare species in multivariate bioassessments has been debated widely. Advocates for the inclusion of rare species argue that rare species are useful indicators of environmental stressors and their removal may result in the unnecessary loss of ecological information (i.e. biological argument; Cao et al. 1999; 2001). On the other hand, advocates for the exclusion of rare species argue that when rare species are removed from analysis, the resultant analyses (and conclusions) are not altered (i.e. statistical argument; Gauch 1982; Marchant 1999; 2002). To date, the debate surrounding the removal of rare species in bioassessments has remained controversial and resolution is needed. One of the difficulties with assessing the importance of removing rare species in bioassessments is the lack of context from which to judge the consequences of the decision. For example, how can one evaluate whether the inclusion of rare species provides redundant information with more common species or provides undue influence (Gauch 1982; Marchant 1999; 2002; Bailey et al. 2004)? Alternatively, how can one determine whether rare species are more sensitive to ecosystem stress than more common species (Faith and Norris 1989; Cao et al. 1999; 2001)? Here, both the statistical and biological arguments were tested as separate hypotheses using data collected from fish species in a well-sampled aquatic system. This study demonstrates contrary to the previously held notions that rare species provide redundant information as compared to more common species, or unduly influence multivariate analyses (Marchant 1999; Marchant et al. 2006); neither was supported. In the case of fish species in the Sydenham River, the hypothesis for the biological argument for the inclusion of rare species in bioassessments was supported, while the hypothesis for the statistical argument was not. The removal of rare species may have large biological consequences for bioassessments. First, as rare species may not be as rare as perceived simply as a result of sampling bias (Preston 1948, Resh et al. 2005; Arscott et al. 2006), the removal of rare species may limit the number of species from which to assess the biological community. Second, and perhaps more importantly, the removal of rare species may fundamentally change conclusions of multivariate bioassessments. In this study, when rare species were removed from the analysis, sites impacted by this removal shifted in multivariate space to a greater degree than those not directly changed by this decision (e.g. nine fold change between the full data set and when single occurring
59
species were removed; Fig. 3.3). Interestingly, these results were reduced as more species were removed (Fig. 3.3), demonstrating that as species are removed from the analyses, site-level species assemblages become more homogenized and differences across multivariate analyses are minimized. Therefore, the removal of rare species may also remove important indicators of ecosystem stress. In the case of the Sydenham River, large-scale agricultural activity and increases in turbidity has led to declines in several species (Staton et al. 2003; Poos et al. 2007), with issues of turbidity shown to be more related to rare species (e.g. species at risk; Poos et al. 2008). Despite claims to the contrary (e.g. Marchant et al. 2006) this study demonstrates a scenario where the choice of removing rare species can have limited statistical effect on the overall multivariate analyses but, at the same time, have large biological effect on particular observations (i.e. sampling sites). Therefore, the assumption that more common species can sufficiently define impact or reflect the response of the whole community may not be justified (Marchant 2002; Marchant et al. 2006). The importance of choices in multivariate analyses need to be better justified for bioassessments. For example, this study demonstrates that the removal of rare species had similar (and often less) influence in multivariate analyses as other choices inherent in its calculation. Previous research has noted that differences in similarity coefficients (Jackson et al. 1989; Jackson 1997; Cao et al. 1998; Legendre and Gallagher 2001; Podani 2005; Podani and Schmera 2007; Poos et al. 2009) and ordination techniques (Jackson 1993; Podani 2000; Heino et al. 2003) can lead to divergent results. Researchers often select methods based on past experience and assume that the resultant summary adequately models the underlying data, or they choose solutions that are most interpretable with regard to a priori hypotheses (Jackson et al. 1989; Jackson 1997; Podani 2000). This approach may have severe consequences for the ultimate goal of inferring community responses for bioassessment. Here, comparisons of choices inherent in multivariate analyses demonstrated that choices, such as ordination methods (e.g. NMDS), can provide largely divergent results (Fig. 3.2; Table 3.2). As a result, the removal of rare species may be less of a concern than previously noted (e.g. Marchant 1999; Marchant et al. 2006), whereas other choices (e.g. type of ordination) may be more important. These results indicate that researchers must be mindful of the statistical decisions they make including ordination technique, similarity coefficient and the exclusion of rarely sampled species, as each choice may have potential to influence community responses and meaningful conclusions. Other issues,
60
such as sample size (Cao et al. 2001b; 2002), seasonal effects (Furse et al. 1984), experience (Metzeling et al. 2003), data standardization (Jackson 1993; Cao et al. 1999), taxonomic resolution (Arscott et al. 2006) and data quality (Cao et al. 2003; Nichols and Norris 2006), also require adequate justification. There are issues with the inclusion of rare species in multivariate analyses that should be taken into account for the bioassessment of aquatic communities that this study could not assess. Species found infrequently, but with varying abundances have been shown strongly to influence multivariate analyses (Legendre and Legendre 1998). Researchers who wish to minimize the impacts of rare species can choose from a variety of data transformations that can downweight the influence of rare species, although caution is warranted when using small sample sizes (Jackson et al. 1989; Jackson 1993; Cao et al. 1999). For example, Legendre and Gallagher (2001) have suggested the use of Helinger transformation for reducing the impacts of rare species. In addition, authors can choose to reduce their data from abundance to presenceabsence (as done here), which represents the strongest form of standardization and is less likely to influence analyses (Cao et al. 1999) than many other decisions. The decision to include or remove rare species will be context dependant. This study demonstrates that, contrary to the common practice of removing rare species in multivariate analysis, that the impact of leaving rare species in may be minimal relative to other methodological choices while maintaining important site-level information on species, including some with species at risk. Ultimately, the decision to include or remove rare species should be justified by the goals of the bioassessment. In cases, such as the Sydenham River, where rare species are used as targets for evaluating ecosystem recovery, rare species should remain in the analyses as their inclusion did not alter the results (e.g. levels of variation were similar among All, M1 and M5 treatment groups; Table 2) and they also represent components of the community being assessed. Naturally, researchers wish to limit their data to reflect the most appropriate number of species, the most practical similarity coefficient and the most useful ordination technique; however, no such set of criteria exists. One alternative is to use a consensus approach where several methods (and choices within methods) are used and compared (Green 1979; Jackson et al. 1989; Jackson 1993). If the methods produce similar results then one can have greater confidence that the results are more robust and representative rather than being dominated by the set of choices used
61
in the analysis (Jackson 1993). In summary, better justifications of all of decisions in analyses are needed to ensure bioassessments are rigorous.
Acknowledgements Funding was provided by NSERC and OGS Scholarships to M.S.P, an NSERC Discovery Grant to D.A.J, and funding from the Ontario Ministry of Natural Resources and University of Toronto. This manuscript was greatly enhanced by conversations with C. Harpur.
References Arscott, D. B., J. K. Jackson, and E. B. Kratzer. 2006. Role of rarity and taxonomic resolution in a regional and spatial analysis of stream macroinvertebrates. Journal of the North American Benthological Society 25:977-997. Austin, M.P. and P. Greig-Smith. 1968. The application of quantitative methods to vegetation survey. II. Some methodological problems of data from rain forests. Journal of Ecology 56:827–844 Bailey, R. C., R. H. Norris, and T. B. Reynoldson. 2004. Bioassessment of freshwater ecosystems: Using the reference condition approach. Springer, New York. Bailey, R. C., M. G. Kennedy, M. Z. Dervish, and R. M. Taylor. 2008. Biological assessment of freshwater ecosystems using a reference condition approach: comparing predicted and actual benthic invertebrate communities in Yukon streams. Freshwater Biology 39:765774. Barbour, M.T., J. Gerritsen, B.D. Snyder, and J.B. Stribling. 1999. Rapid bioassessment protocols for use in streams and wadeable rivers: Periphyton, benthic macroinvertebrates and fish. 2nd edition. EPA-841-B-99-002. Environmental Protection Agency. Borcard, D., P. Legendre, and P. Drapeau. 1992. Partialling out the Spatial Component of Ecological Variation. Ecology 73:1045-1055.
62
Borcard, D., P. Legendre, C. Avois-Jacquet, and H. Tuomisto. 2004. Dissecting the spatial structure of ecological data at multiple scales. Ecology 85:1826-1832. Bohlin, T., S. Hamrin, T.G. Heggberget, G. Rasmussen, and S.J. Saltveit. 1989. Electrofishing theory and practice with special emphasis on salmonids. Hydrobiologia 173: 9-43. Boulton, A. J., 1999. An overview of river health assessment: philosophies, practice, problems and prognosis. Freshwater Biology 41: 469–479. Bowman, M. F., and K. M. Somers. 2005. Considerations when using the reference condition approach for bioassessment of freshwater ecosystems. Water Quality Research Journal of Canada 40:347-360. Bowman, M. F., and K. M. Somers. 2006. Evaluating a novel Test Site Analysis (TSA) bioassessment approach. Journal of the North American Benthological Society 25:712727. Boys, C. A. and M. C. Thoms, 2006. A large-scale, hierarchical approach for assessing habitat associations of fish assemblages in large dryland rivers. Hydrobiologia 572: 11–31. Cao, Y., and D.D. Williams. 1999. Rare species are important in bioassessment (Reply to the comment by Marchant). Limnology and Oceanography 44:1841-1842. Cao Y., D.D. Williams and N.E. Williams. 1998. How important are rare species in community ecology and bioassessment. Limnology and Oceanography 43: 1403–1409. Cao, Y., D. D. Williams, and N. E. Wiliams. 1999. Data transformation and standardization in the multivariate analysis of river water quality. Ecological Applications 9:669-677. Cao Y., D.P. Larsen, and R.S. Thorne. 2001a. Rare species in multivariate analysis for bioassessment: some consideration. Journal of the North American Benthological Society 20: 144–153. Cao Y., D.P. Larsen, and R.M. Hughes. 2001b. Evaluating sampling sufficiency in fish assemblage survey: a similarity-based approach. Canadian Journal of Fisheries and Aquatic Sciences 58: 1782–1793.
63
Cao Y., D.D. Williams, and D.P. Larsen. 2002. Comparison of ecological communities-the problem of sample representativeness. Ecological Monographs 72: 41–56. Cao, Y., C.P. Hawkins, and M.R. Vinson. 2003. Measuring and controlling data quality in biological assemblage surveys with special reference to stream benthic macroinvertebrates. Freshwater Biology 48: 1898–1911 Clarke, K. R., and R. H. Green. 1988. Statistical design and analysis for a 'biological effects' study. Marine Ecology Progress Series 46: 213-226. COSEWIC. 2007. Canadian Species at Risk. Committee on the Status of Endangered Wildlife in Canada, Ottawa, Ontario. Day, J.H., J.G. Field, and M.P. Montgomery. 1971. The use of numerical methods to determine the distribution of the benthic fauna across the continental shelf of North Carolina. Journal of Animal Ecology 40: 93–125. DeVelice, R. L., J. W. DeVelice, and G. N. Park. 1988. Gradient analysis in nature reserve design: a New Zealand example. Conservation Biology 2: 206-217. Digby, P. G. N. and R.A. Kempton. 1987. Multivariate analysis of ecological communities. Chapman & Hall, London. Dray, S., P. Legendre, and P. R. Peres-Neto. 2006. Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecological Modelling 196:483-493. Faith, D. P., and R. H. Norris. 1989. Correlation of environmental variables with patterns of distribution and abundance of common and rare freshwater macroinvertebrates. Biological Conservation 50:77-98. Fausch, K. D., J. Lyons, J. R. Karr, and P. L. Angermeir. 1990. Fish communities as indicators of environmental degradation. American Fisheries Society Symposium 8:123-144. Furse, M. T., D. Moss, J.F. Wright, and P.D. Armitage. 1984. The influence of seasonal and taxonomic factors on the ordination and classification of running-water sites in Great
64
Britain and on the prediction of their macroinvertebrate communities. Freshwater Biology 14: 257-80. Gaston, K.J. 1994. Rarity. Chapman and Hall Press, London. Gauch, H. G. 1982. Multivariate analysis in community ecology. Cambridge University Press, Cambridge. Goslee, S. and D. Urban. 2007. Ecodist: Dissimilarity-based functions for ecological analysis. R package version 1.1.3. Gower, J.C. 1966. Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika 53 : 325-338. Gower J. C. 1971. A general coefficient of similarity and some of its properties. Biometrics 27: 857–874. Gower, J. C., and P. Legendre. 1986. Metric and Euclidean properties of dissimilarity coefficients Journal of Classification 3: 5-48. Green, R. H. 1979. Sampling design and statistical methods for environmental biologists. Wiley, New York. Hannaford, M. J., and V. H. Resh. 1995. Variability in macroinvertebrate rapid-bioassessment surveys and habitat assessments in a Northern California stream. Journal of the North American Benthological Society 14:430-439. Hawkins, C. P., R. H. Norris, J. N. Hogue, and J. W. Feminella. 2000. Development and evaluation of predictive models for measuring biological integrity of streams. Ecological Applications 10:1456-1477. Heino, J., T. Muotka, H. Mykra, R. Paavola, H. Hamalainen, and E. Koskenniemi. 2003. Defining macroinvertebrate assemblage types of headwater streams: implications for bioassessment and conservation. Ecological Applications 13:842-852.
65
Hewlett, R. 2000. Implications of taxonomic resolution and sample habitat for stream classification at a broad geographic scale. Journal of the North American Benthological Society 19:352-361. Hirst, C. N., and D. A. Jackson. 2007. Reconstructing community relationships: the impact of sampling error, ordination approach, and gradient length. Diversity and Distributions 13:361-371. Jackson, D.A. 1993. Multivariate analysis of benthic invertebrate communities: the implication of choosing particular data standardizations, measures of association, and ordination methods. Jackson, D.A. 1995. PROTEST: A PROcrustean randomization test of community environment concordance. Ecoscience 2: 297-303. Jackson, D.A. 1997. Compositional data in community ecology: the paragidm or peril of proportions? Ecology 78: 929-940. Jackson D. A., K.M. Somers, and H.H. Harvey. 1989. Similarity coefficients: measures of cooccurrence and association or simply measures of occurrence? American Naturalist 133: 436-453. Joy, M. K., and R. G. Death. 2000. Development and application of a predictive model of riverine fish community assemblages in the Taranaki region of the North Island, New Zealand. New Zealand Journal of Marine and Freshwater Research 34:241-252. Jurasinski, G. 2007. Simba: A collection of functions for similarity calculation of binary data. R package version 0.2-5. Kennard, M. J., B. J. Pusey, A. H. Arthington, B. D. Harch and S. J. Mackay, 2006. Development and application of a predictive model of freshwater fish assemblage composition to evaluate river health in eastern Australia. Hydrobiologia 572: 33–57. King, J.R., and D.A. Jackson. 1999. Variable selection in large environmental data sets using principal components analysis. Environmetrics 10: 67-77.
66
Krebs, C.J. 1998. Similarity coefficients and cluster analysis. Pages 293-323 in Ecological methodology. Harper and Row, New York. Legendre, P., and E.D. Gallagher. 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129: 271-280 Legendre, P. and L. Legendre. 1998. Numerical ecology. Elsevier Science B.V., Amsterdam. 853 pp. Mace, G.M. 1994. Classifying threatened species: means and ends. Philosophical Proceedings of the Royal Society of London. Series B. 344: 91-97. Marchant, R. 1990. Robustness of classification and ordination techniques applied to macroinvertebarte communities from the La Trobe River, Victoria. Australian Jouranl of Marine and Freshwater Research 41: 493-504 Marchant, R., A. Hirst, R.H. Norris, R. Butcher, L. Mixzeling, and D. Tiller. 1997. Classification and ordination of macroinvertebrate assemblages from running waters in Victoria, Australia. Journal of North American Benthological Society 16: 664-681. Marchant, R. 1999. How important are rare species in aquatic ecology and bioassessment? A comment on the conclusions of Cao et al. Limnology and Oceanography 44: 1840–1841. Marchant, R. 2002. Do rare species have any place in multivariate analysis for bioassessment? Journal of the North American Benthological Society 21: 311-313 Marchant, R, R.H. Norris, and A. Milligan. 2006. Evaluation and application of methods for biological assessment of streams: summary of papers. Hydrobiologia 572: 1-7. Margules, C. R. 1986. Conservation evaluation in practice. Pages 297-314 in Wildlife Conservation Evaluation. Chapman and Hall, London. McCune, B. and J.B. Grace. 2002. Analysis of ecological communities MjM Software Design, Gleneden Beach, Oregon.
67
McGarigal, K., S. Cushman, and S. Stafford. 2000. Multivariate statistics for wildlife and ecology research. New York: Springer-Verlag. Metzeling, L., D. Tiller, P. Newall, F. Wells, and J. Ree. 2006. Biological objectives for the protection of rivers and streams in Victoria, Australia. Hydrobiologia 572:287-299. Mugodo, J., M. Kennard, P. Liston, S. Nichols, S. Linke, R. H. Norris and M. Lintermans. 2006. Local stream habitat variables predicted from catchment scale characteristics are useful for predicting fish distribution. Hydrobiologia 572: 59–70. Nichols, S. and R. H. Norris, 2006. River condition assessment may depend on the sub-sampling method: field live-sort versus laboratory sub-sampling of invertebrates for bioassessment. Hydrobiologia 572: 195–213. Norris, R. H. 1995. Biological monitoring: the dilemma of data analysis. Journal of the North American Benthological Society 14:440-450. Norris, R. H., and C. P. Hawkins. 2000. Monitoring river health. Hydrobiologia 435:5-17. Norris, R. H., and M. C. Thoms. 1999. What is river health? Freshwater Biology 41:197-209. Paavola, R., T, R. Muotka, R. Virtanen, J. Heino, D.A. Jackson, and A. Maki-Petays. 2006. Spatial scale affects community concordance among fishes, benthic macroinvertebrates and bryophytes in boreal streams. Ecological Applications 16:368-379. Paukert, C.P., and T.A. Wittig. 2002. Applications of multivariate statistical methods in fisheries. Fisheries 27: 16-22. Peres-Neto, P. R., and D. A. Jackson. 2001. How well do multivariate data sets match? Evaluating the association of multivariate biological data sets: comparing the robustness of Mantel test and a Procrustean superimposition approach. Oecologia 129:169-178. Peres-Neto, P. R., D. A. Jackson, and K. M. Somers. 2003. Giving meaningful interpretation to ordination axes: assessing loading significance in principal component analysis. Ecology 84: 2347-2363.
68
Podani, J. 2000. Introduction to the exploration of multivariate biological data. Backhuys Publishers, Leiden, Netherlands. Podani, J. 2005. Multivariate exploratory analysis of ordinal data in ecology: pitfalls, problems and solutions. Journal of Vegetation Science 16: 497-510. Podani, J. and D. Schmera. 2006. On dendrogram-based measures of functional diversity. Oikos 115: 179-185. Poos, M.S., N.E. Mandrak, and R.L. McLaughlin. 2008. A practical framework for selecting among single species, multi-species and ecosystem-based recovery plans. Canadian Journal for Fisheries and Aquatic Science 65: 2656-2666. Poos, M.S., N.E. Mandrak, and R.L. McLaughlin. 2007. The effectiveness of two common sampling methods for sampling imperiled freshwater fishes. Journal of Fish Biology 70: 691-708. Poos, M.S., S.C. Walker, and D.A. Jackson. 2009. Functional diversity indices can be driven by methodological choices and species richness. Ecology 90: 341-347. Preston, F. W. 1948. The commonness, and rarity, of species. Ecology 29:254-283. Oksanen, J., R. Kindt, P. Legendre, B. O'hara, G.L. Simpson, M. Henry, and H. Stevens. 2008. Vegan: Community ecology package. R package version 1.11-0. Olden, J. D., D. A. Jackson, and P. R. Peres-Neto. 2001. Spatial isolation and fish communities in drainage lakes. Oecologia 127:572-585. OMNR 2007.Stream assessment protocol for southern Ontario. Picton, Ontario: Ontario Ministry of Natural Resources. Orloci, L. and M.M. Mukkattu. 1973. The effect of species number and type of data on the resemblance structure of a phytosociological collection. Journal of Ecology 61:37–46.
69
Rabinowitz, D., S. Cairns, and T. Dillon. 1986. Seven forms of rarity and their frequency in the flora of the British Isles. Pages 182–204 in M. E. Soule (ed.). Conservation biology: the science of scarcity and diversity. Sinauer Associates, Sunderland,Massachusetts Resh, V. H., L. A. Beche, and E. P. McElravy. 2005. How common are rare taxa in long-term benthic macroinvertebrate surveys? Journal of the North American Benthological Society 24:976-989. Reynolds, J. B. 1996. Electrofishing. Pages 221-253 in B. R. Murphy and D. W. Willis (ed.). Fisheries techniques. American Fisheries Society, Bethesda, Maryland. Reynoldson, T. B., D. M. Rosenberg, and V. H. Resh. 2001. Comparison of models predicting invertebrate assemblages for biomonitoring in the Fraser River catchment, British Columbia. Canadian Journal of Fisheries and Aquatic Sciences 58:1395-1410. Roberts, D.W. 2008. Ordist: Ordination and Multivariate Analysis for Ecology. R package version 1.3. Rundle, H. D., and D. A. Jackson. 1996. Spatial and temporal variation in littoral-zone fish communities. Canadian Journal of Fisheries and Aquatic Sciences 53:2167-2176. Simpson, J. C., and R. H. Norris. 2000. Biological assessment of river quality: development of AUSRIVAS models and outputs. Pages 125-142 in J. F. Wright, D. W. Sutcliffe, and M. T. Furse (ed.). RIVPACS and Similar Techniques for Assessing the Biological Quality of Freshwaters. Biological Association and Environment Agency, Ambleside, Cumbria, U.K. Sloane, P. I. W. and R. H. Norris, 2003. Relationship of AUSRIVAS-based macroinvertebrate predictive model outputs to a metal pollution gradient. Journal of the North American Benthological Society 22: 457–471. Somers, K. M., R. A. Reid, and S. M. David. 1998. Rapid biological assessments: how many animals are enough? Journal of the North American Benthological Society 17:348-358. Stanfield, L.W. and M.L. Jones. 1998. A comparison of full-station visual and transect-based methods of conducting habitat surveys in support of habitat suitability index models for southern Ontario. North American Journal of Fisheries Management, 18: 657-675.
70
Staton, S.K., A. Dextrase, J.L. Metcalfe-Smith, J. DiMaio, M. Nelson, J. Parish, and E. Holm. 2003. Status and trends of Ontario's Sydenham River ecosystem in relation to aquatic species at risk. Environmental Monitoring and Assessment 88: 283-310. Wallace, J. B., J. W. Grubaugh, and M. R. Whiles. 1996. Biotic indices and stream ecosystem processes: results from an experimental study Ecological Applications 6:140-151. Wartenberg, D., S. Ferson, and F.J. Rohlf. 1987. Putting things in order: A critique of detrended correspondence analysis. American Naturalist 129: 434-448. Webb, L. J., J.G. Tracey, W.T. Williams, and G.N. Lance. 1967. Studies in the numerical analysis of complex rain forest communities. II. The problem of species sampling. Journal of Ecology 55: 525– 538. Wright, J. F., D. W. Sutcliffe, and M. T. Furse. 2000. Assessing the biological quality of fresh water: RIVPACS and other techniques. Freshwater Biological Association, Cumbria, UK.
71
Appendices Appendix 3.1 – Site-level effects of methodological choices in bioassessments. Shown are the ratios between mean site-level vector residuals from Procrustes analyses of sites having species removed and those sites having no species removed. Mean site-level vector residual values were separated by sites which had rare species removed (M1: n=2; M5: n=19; and M10: n=63); and compared with those sites that not.
PCoA-J PCoA- Φ PCoA-RR PCoA-SM NMDS-J NMDS- Φ NMDS-RR NMDS-SM CA PCA* Overall Average
M1 13.83 2.40 9.64 12.61 8.66 9.25 7.91 9.58 13.97 12.61 9.76
M5 M10 2.58 1.10 1.45 1.03 1.70 1.26 3.47 1.11 2.41 0.93 2.00 1.03 1.91 1.27 1.21 0.79 2.46 1.61 3.47 1.11 2.13 1.13
Note: Abbreviations of various treatments are carried forward from Table 1. *PCA values are shown for comparison, but are not included in overall averages.
72 Appendix 3.2 – Summary of three-dimensional ordination results. Shown are eigenvalues for Principal Coordinates Analyses (PCA) and Correspondence Analyses (CA), with percent variance explained shown in parentheses. Stress values are shown for Non-Metric Multidimensional Scaling (NMDS).
Ordination Technique A-PCoA-J 1-PCoA-J 5-PCoA-J 10-PCoA-J A-PCoA- Φ 1-PCoA- Φ 5-PCoA- Φ 10-PCoA- Φ A-PCoA-RR 1-PCoA-RR 5-PCoA-RR 10-PCoA-RR A-PCoA-SM 1-PCoA-SM 5-PCoA-SM 10-PCoA-SM A-NMDS-J 1-NMDS-J 5-NMDS-J 10-NMDS-J A-NMDS- Φ 1-NMDS- Φ 5-NMDS- Φ 10-NMDS- Φ A-NMDS-RR 1-NMDS-RR 5-NMDS-RR 10-NMDS-RR A-NMDS-SM 1-NMDS-SM 5-NMDS-SM 10-NMDS-SM ACA--1CA--5-CA--10-CA---
Axis 1
Axis 2
Axis 3
3.43 (10.67) 3.43 (10.67) 3.49 (10.99) 3.69 (12.01) 4.87 (15.05) 4.93 (15.09) 5.23 (15.51) 6.19 (17.25) 2.21 (13.00) 2.23 (13.12) 2.32 (13.65) 2.58 (15.18) 2.03 (15.63) 2.10 (15.67) 2.36 (16.04) 3.06 (17.56) 0.4158 0.4183 0.4179 0.4235 0.4186 0.4219 0.4245 0.4492 0.3994 0.3990 0.4012 0.3950 0.3830 0.3870 0.3882 0.4120 0.32 (12.17) 0.31 (12.24) 0.30 (14.03) 0.23 (14.88)
2.79 (19.37) 2.83 (19.48) 2.89 (20.05) 3.11 (22.12) 3.97 (27.31) 4.01 (27.37) 4.21 (28.00) 4.40 (29.51) 2.06 (25.12) 2.07 (25.29) 2.15 (26.29) 2.30 (28.71) 1.65 (28.36) 1.71 (28.41) 1.92 (29.07) 2.37 (31.16) --------------------------------0.22 (20.62) 0.22 (21.10) 0.22 (24.31) 0.21 (28.98)
1.67 (24.58) 1.68 (24.71) 1.70 (25.40) 1.70 (27.62) 2.14 (33.93) 2.17 (34.00) 2.30 (34.82) 2.61 (36.79) 1.62 (27.02) 1.63 (27.21) 1.65 (28.23) 1.68 (30.68) 0.86 (34.95) 0.89 (35.02) 1.00 (35.85) 1.25 (38.36) --------------------------------0.15 (26.17 0.14 (26.83) 0.14 (30.99) 0.12 (36.62)
73
Chapter 4: Contrasting direct versus indirect dispersal in metapopulation viability analyses
Abstract Species dispersal is a central component of metapopulation models. Spatially realistic metapopulation models, such as stochastic patch-occupancy models (SPOMs), quantify species dispersal using indirect estimates of colonization potential based on inter-patch distance. In this study, indirect parameterization of SPOMs was compared with dispersal and patch dynamics quantified directly from empirical data. For this purpose two metapopulations of an endangered minnow, redside dace (Clinostomus elongatus), were monitored using mark-recapture techniques across 43 patches, re-sampled across a one year period. More than 2,000 fish were marked with visible implant elastomer tags coded for patch location and dispersal and patch dynamics were monitored. Direct and indirect parameterization of SPOMs provided qualitatively similar rankings of viable patches; however, there were differences of several orders of magnitude in the estimated intrinsic mean times to extinction, from 24 and 148 years to 362 and >100,000 years, depending on the population. In several cases, patches that were in close proximity (high colonization potential) that were not used by redside dace. This study demonstrates the importance of incorporating species and patch-specific data directly into metapopulation models, especially given heterogeneous landscapes.
Keywords: metapopulations, dispersal, population viability analysis, stochastic patch-occupancy models, parameterization.
74
Introduction Species dispersal is a central component in the study of spatially structured populations. At a landscape scale, population viability strongly depends on individual dispersal allowing recolonisation of empty habitats or patches (Hanski 1999a). For this reason, species dispersal is considered the ‘glue’ for maintaining local populations within a network of suitable habitats (Hansson 1991; Vandewoestijne et al. 2004). The degree of dispersal has an impact on local population dynamics, on gene flow, and on adaptation to local conditions. For example, low dispersal can foster isolation and local adaptations (Thomas and Hanski 1997; Resetarits et al. 2005). Alternatively, high species dispersal can have a stabilizing effect on metapopulation dynamics (Hanski 1999a). Many species with spatially structured populations are in decline, and population viability models provide a statistical evaluation of species viability to facilitate informed management decisions (Ackakaya 2000; Nowicki et al. 2007). Metapopulation viability analyses provide a spatially realistic evaluation of the local population structure (Hanski 2001; March 2008). A metapopulation is defined as a system of local populations (patches) connected by dispersing individuals (Hanski and Gilpin 1991). By quantifying patch dynamics, metapopulation viability analyses can be used to understand better the importance of ecological processes such as speciesspecific dispersal, patch quality and landscape influences (Moilanen and Hanski 1998), and to enhance management through evaluation of minimum amount of habitat or population size needed to maintain viability (Hanski 1999b; Robert 2009). Understanding how the parameterization of metapopulation viability analyses may impact results can inform managers as to the potential areas of concern when developing management decisions. One popular type of metapopulation viability analysis is stochastic patch-occupancy models (i.e. SPOMs), which have been used extensively to model the viability of spatially structured populations (Hanksi 1999; Moilanen 1999). For example, SPOMs were used in studies of species with conservation concern, such as capercaillie (Grimm and Storch 2000), American pika (Moilanen et al. 1998) and Glanville fritillary, and silver spotted skipper butterflies (Hanski et al. 1994). As SPOMs provide a simplification over traditional populationviability analyses (Akayaka and Sjögren-Gulve 2000), they do not require demographic or stage
75
data, but only occupancy, colonization and extinction rates, which can be estimated readily from empirical data (Hanski 1994,1999; Moilanen 1999, 2004; Grimm et al. 2004). The influence of dispersal on metapopulation viability is often evaluated using some approximation of colonization potential (Verbroom et al. 1993; Hanski and Gilpin 1997; Frank and Wissel 2002; Heinz et al. 2005). The easiest approach to describe colonization potential (i.e. patch accessibility) is as a function of distance between a starting patch to a target patch and the ability of species to disperse (Hanksi 1994; Hanski et al. 1996; Heniz et al. 2005). This relationship can be quantified in several ways; however, most often this estimation is done by assuming that colonization potential declines exponentially with distance (i.e. exponential decay; Hanski 1994; Vos et al. 2001; Frank and Wissel 2002). It is uncertain how well the assumption of exponential decay can model species-specific dispersal (Hill et al. 1996; Baguette et al.2000; Heinz et al. 2005). Whether simple formulae are adequate in describing species- and patchspecific movement in metapopulation models remains an open question (Heinz et al. 2005; Marsh 2008). The overall aim of this study is to assess whether direct parameterization of species dispersal can impact estimates of metapopulation viability. For this assessment, a detailed mark-recapture survey of the endangered fish, redside dace (Clinostomus elongatus), was conducted in the Greater Toronto Area, Ontario, Canada. The redside dace is a spatially structured, pool-dwelling species that is undergoing declines in the majority of its range due to impacts from urbanization (COSEWIC 2007; Poos and Jackson, submitted). Two locations on the Rouge River were monitored, including one location on Leslie Tributary, and the other location on Berczy Creek (Figure 4.1). These locations were shown previously to have among the highest abundances of redside dace recently sampled across its entire Canadian range (Reid et al. 2008).
76
A) Leslie Tributary
B) Berczy Creek
L20 L19
L18 L17
B17 B16 B15 B14
L16 L13 L12
L15 L14 L11
B11 B8
L10
L7
B7 L9 L8
B6 B5 B2
L6 L5 L3
N 100m
B4
B22 B20
B23 B21
B18 B19 B13 B12 B10 B9
B3 B1
100m
N
L4 L2 L1
N 0 2 4 8 Figure 4.1 – Study sites on Rouge River, Ontario where redside dace (Clinostomus elongatus) dispersal and patch dynamics were monitored. Study locations: A) Leslie Tributary, and B) Berczy Creek, were subdivided into extensive sites (black), where redside dace were tagged with a color-coded visual implant elastomer tag, and extended sites (grey), which were monitored for tag movement.
77
Methods The metapopulation dynamics of redside dace were monitored by enumerating the dispersal of tagged individuals on monthly intervals during a one-year period. For this study, each location was sub-divided into two areas: intensively monitored sites where individuals were tagged; and, extended sites, that were beyond those areas where fish were not tagged but where tagged fish could have potentially moved (Figure 4.1). As meta-populations can be defined in a number of ways (Hanski 1999a), a metapopulation was defined as an assemblage of local populations inhabiting spatially distinct habitat patches (Moilanen and Hanski 1998). Redside dace live primarily in clear water with well-defined pools (COSEWIC 2007); therefore, each spatially distinct pool segregated by a well-defined riffle (e.g. a passable, but natural migratory barrier) was selected as a habitat patch. Leslie Tributary was sub-divided into 20, connected and distinct patches, with 10 intensive sites and five extended sites on both upstream and downstream ends. Similarly, Berczy Tributary was sub-divided into 13 intensive sites with five extended sites on each of the upstream and downstream ends (Figure 4.1). Sampling was conducted using multiple-pass depletion surveys at each pool. Using a twentyfoot bag seine (1/4” mesh), each site was surveyed until depletion of redside dace, with a minimum of three sampling events conducted at each site per time period. At each pool, redside dace were implanted with visual implant elastomer (VIE) tags colour-coded for their location (Plates 1 and 2). VIE tags were chosen because they had good tag retention and negligible effects on survival, growth and behavior when used on other species (Dewey and Ziegler 1996; Goldsmith et al. 2003; Walsh and Winkelman 2004). Tags were injected subcutaneously near the anal fin on the ventral surface (See Plates 1 and 2). All redside dace were held in welloxygenated flow-through bins for 2-4 hours to monitor for potential physiological stress, and then returned to the river at the site of capture. Both intensive and extended sites were resampled for redside dace at monthly intervals, except under winter-ice conditions (NovemberMarch) and when redside dace were spawning (June) to not disrupt this important life stage for an endangered species. All redside dace sampled in a recapture event were examined for the presence of a VIE tag. Redside dace dispersal and metapopulation dynamics were tracked and mark-recapture data were recorded. If redside dace were re-captured at a new location, they were subsequently tagged posterior to the existing tag, with a new colour code for the recapture
78
location. The dispersal patterns, such as average distance dispersed, and proportion of stationary tags, of the metapopulations were compared using non-parametric Mann-Whitney U-tests and log-linear models (G-tests with Yates continuity correction; Zar 1999) respectively, in the R language v2.80 (R Development Team 2008).
Determining Metapopulation Viability The viability of the two stream metapopulations was quantified using direct and indirect parameterization of stochastic patch-occupancy model (SPOM). SPOMs use a time-continuous Markov-chain model (Hanski 1999; Grimm et al. 2004). Each patch (i) is assumed to be in one of two states, vacant (xi =0) or occupied (xi=1). Changes in these states can occur from a patch becoming vacant due to local extinction (xi: 1Æ 0) or correlated extinction (i.e. regional stochasticity) from another patch (xj, xi: 1Æ0). Alternatively, a vacant patch can become occupied (xi: 0Æ1) via colonization from another patch (j). The state of the whole metapopulation (xi, … xn) is given by a vector of states xi of these individual patches. The metapopulation models were parameterized using a combination of indirect (i.e. dispersal ability from empirical data and incidence functions; Hanski 1999) and direct parameterization of SPOMs (using rates of actual patch colonization from empirical data alone). The models were quantified as follows.
Colonization Rate: Colonization between two patches i and j (bij) was defined using an incidence-function model (Hanski 1994):
·
· exp
/
I
, where y is a parameter, and Mi is the number of
emigrants from pool i. The mean number of emigrants leaving a pool was estimated using data from the tagging study. To account for the potential uncertainty with missing emigrants leaving a patch, the probability of detection at each pool (PDi) was quantified using maximum likelihood from the n-pass depletion surveys (Zippin 1956, 1958) with the Bayesian modification from Carle and Strub (1978). These were coded in the R v2.80 (R Development Team 2008) using the fisheries-assessment package FAS (Ogle 2009). The total number of emigrants leaving each patch per year was measured as
∑
…
1
1
, where mi is the
uncorrected number of emigrants. Similar to most metapopulation models, a distance-based dispersal kernel using a negative exponential decay was used, where exp
/
, and dij is the
79
distance from patch i to patch j and d0 is the mean dispersal ability of redside dace. This type of dispersal kernel has been used extensively in metapopulation models and assumes that patch accessibility is dependent on distance (Hansson 1991; Hanski et al. 1996; Hokit et al. 1999; Moilanen 2004). As dispersal data are quantified indirectly using patch area or depth, hereafter, this approach is referred to as indirect parameterization of the metapopulation model. This type of indirect parameterization allows researchers to extrapolate relationships in patch occupancy, often by using species life-history characteristics, without the need of labor-intensive field studies (Hanski, 1994; Moilanen 2004; Heinz et al. 2005). A non-linear (polynomial) regression of mean distance of dispersal of recaptured fish through time was used to identify potential dispersal across the patches. The fits of these non-linear regressions were highly significant (Leslie Tributary; r2 = 0.92, p < 0.01, Berczy Creek; r2 = 0.88, p < 0.01) and the average distance dispersed (dI) for a one-year period (one time step in final SPOM model) was 210m for redside dace in Berczy Creek and 150m for Leslie Tributary (Table 1). As there were several consecutive surveys, it was possible to estimate y from the number of transitions (i.e. an empty patch becoming occupied and vice versa; Hanski 1999a). For this, a GLM procedure was used which considered multiple snapshots of the sampling events using a binomial distribution and logistic-link function developed in the R programming language v2.80 (R Development Team 2008) using the incidence function (see Oksanen 2004 for details). The value of the y parameter for Leslie Tributary was 0.0816 and Berczy Creek was 0.0713. Finally, for dij, a distance matrix was measured using the river distance between patches as calculated by a geographic information system.
Extinction Rate: Extinction rates can be quantified in many ways (Hanksi 1999a). The simplest form of determining extinction rate (Ei) is using the area of the patch (Ai), and given by
,
where e defines the extinction probability of a patch of unit size, and x defines the scaling of the extinction risk with patch area (Hanski 1998; Moilanen 2004). This model assumes that probability of extinction generally depends on population size, which, in turn, is usually given by a simple linear or power function of patch area. This relationship has been demonstrated on both empirical and theoretical grounds (Lande 1993; Foley 1994; Hanski 1994; 1999a; Hanski et al. 1996). Here extinction rates were calculated using patch (i.e pool) depth (d). Patch depth was
80
calculated by taking the average of 60 equidistant point measurements as suggested by the Ontario Stream Assessment Protocol (OMNR 2007). Patch depth was used because redside dace are known to be pool-dwelling species (COSEWIC 2007) and, therefore, depth may be more relevant to model patch dynamics. Indeed, redside dace abundances were more correlated with patch depth (r=0.44, p=0.0018) than patch area (r = 0.39, p = 0.0048; Poos and Jackson, unpublished data). The extinction rate was fitted using an incidence function relating species presence in the patches over time and depth. For this modeling, a GLM procedure was used that considered multiple snapshots of sampling events using a binomial distribution developed in the R programming language v2.80 (R Development Team 2008) using the incidence function (Oksanen 2004). The parameters of the incidence function for Leslie Tributary and Berczy Creek were x = 0.4926, 0.5652, and e = 3.685, 4.187, respectively.
Incorporating Dispersal Directly into the Metapopulation Model: Recent theoretical studies on the impact of species movement have found that it can alter metapopulation viability (Heinz et al. 2005; 2006; Revilla and Wiegand 2008). Therefore, the incidence function models were extended by incorporating species dispersal directly into the metapopulation model using empirical data of patch-specific movement. As dispersal data was directly used to quantify the metapopulation model, hereafter, this is referred to as direct parameterization of the metapopulation model. For colonization rate, a model developed by Frank and Wissel (1998; 2002) was used, which (in this case) is identical to the incidence function model and allows the incorporation of patch dynamics (Grimm et al. 2004; Heinz et al. 2006). This model took into account three processes; emigration of individuals from occupied patches; dispersal to a target patch; and, the establishment of a new subpopulation on the target patch. The rate of colonization, bij, was defined as ·
· 0
.
where, Mi was the number of emigrants leaving the occupied patch i per year
(previously defined), ni was the number of connections from patch i to other patches, rij was the probability of an individual started at patch i successfully dispersing to patch j, and Ij was the number of immigrants needed to establish a new subpopulation (Frank and Wissel 1998; 2002). For the probability of dispersal between patches (rij), a patch-colonization matrix was developed using the empirical tagging results for each time period. As the tags were colour coded for patch
81
location (at each time period), rij was defined empirically as the ratio of fish that started at patch i that dispersed to location j, across all recaptured fish. In this instance, tag loss was not accounted for as no tag-related behavioral response or tagging related mortality was assumed (i.e. this ratio was adequate given equal likelihood of mortality of a tagged fish versus non-tagged fish). Given the monitoring data (>75 hours), only minor amounts of tag-related mortality occurred (<0.01%), and all occurred on the first sampling day (likely due to experience bias). Finally, to quantify Ij, an incidence-function model was developing using the probability a patch persistence (across all time periods) given the starting population size (StPopn) of each patch at the start of the study (t1). For this approach a GLM function was used with binomial distribution and logit link in the R programming language v2.80 (R Development Team 2008; Appendix 4.1).
Regional Stochasticity Regional stochasticity refers to the level of correlated extinctions caused by factors influencing a shared geographic location, such as weather or disease (Lande et al. 1988; Lande 1993; Foley 1994). Regional stochasticity has the ability to impact metapopulation viability by incorporating the influence of the fate of proximal patches. The influence of regional stochasticity was quantified at three levels: 0 (no influence of regional stochasticity); 0.1 (a moderate level of regional stochasticity); and, 0.2 (more severe regional stochasticity).
Comparing Viability of Metapopulations Using Direct versus Indirect Parameterization The ultimate viability of patch (i) was defined using the intrinsic mean time to extinction Tm = 1/λ, determined using the reciprocal value of the overall extinction rate λ calculated using a plot of −ln(1 − P0(t)), where P0 is the probability of extinction at a given time (t) (Verbroom et al. 1991; Grimm et al. 2004; Grimm and Wissel 2004). Intrinsic mean time to extinction has been previously shown to be an adequate currency in assessing the viability of metapopulations and can be easily extracted from simulation data (Frank and Wissel 1998; Grimm and Wissel 2004; Heinz et al. 2006). Transitions in metapopulations were simulated 10,000 times using ‘stochastic time steps’ (Frank et al. 2004; Grimm et al. 2004) of transition probabilities of extinction and colonization rates. For this estimation a manually created sub-routines was created in the software program Meta-X (Frank et al. 2003; Grimm et al. 2004), a metapopulation program flexible for incorporating behavior into metapopulation-viability analysis (Heinz et al. 2006).
82
Results In total, 2,141 redside dace were tagged and monitored across 43 patches during in a one-year period from 2007-2008. Due to logistical issues, the stream systems were not sampled during winter-ice conditions (November-March) or when redside dace were spawning in June given the potential to disrupt spawning activities. Recapture rates for redside dace - calculated as the proportion of fish marked during the preceding marking period that were recaptured - were generally high (>25%) during the initial four monitoring events, ending in October 2007. These numbers were greatly reduced by the following spring, with recapture rates < 10% likely due to high over-winter mortality or the re-distribution of tagged fishes to areas beyond the study location (Table 4.1). In addition, the capture efficiency – as determined by probability of detection using n-pass depletion surveys (Zippin et al. 1954; 1956; Carle and Strub 1978) – was also very high for both study systems: Leslie Tributary (mean 71%); and, Berczy Creek (mean 65.6%, Appendix 4.1). Table 4.1 – Summary of mark-recapture information of the endangered fish the reside dace (Clinostomus elongatus) used to directly parameterize stochastic patch occupancy metapopulation models. Shown are two locations in the Greater Toronto Area, Ontario, Canada: A) Leslie Tributary, and B) Berczy Creek. Note: items denoted with a single asterisk (*) represent a significant difference between populations (MannWhitney U for average dispersal, G-test for stationary tags, p<0.05).
A) Leslie Tributary Time T1: July T2: Aug. T3: Sept. T4: Oct. T5: Apr. T6: May T7: July Avg.
Cum. Recap. Avg. Max. Tags (%) (m) (m) 133 0.46 21* 175 305 0.26 91 227 404 0.28 71 547 483 0.26 139 680 * 503 0.04 196 680 542 0.09 182 649 662 0.10 180* 547 0.14 105
B) Berczy Creek Stationary Tags (%) 0.67 0.19 0.52 0.03* 0.00 0.14 0.17 0.31
Cum. Recap. Avg. Tags (%) (m) 342 0.61 9* 511 0.26 84 770 0.26 66 1,045 0.25 77 1,137 0.03 129* 1,376 0.03 174 1,479 0.06 82* 0.18 53
Max. Stationary (m) Tags (%) 125 0.74 290 0.23 357 0.44 315 0.41* 411 0.00 411 0.10 275 0.29 0.49
Legend: Cum. Tags refers to the cumulative number of tags released in the population; Recap (%) refers to the recapture rate for cumulative tags, Avg. (m) refers to the average distance the recaptured tags were captured at, Max. (m) refers to the maximum distance recaptured tags were captures at, and Stationary Tags (%) refers to the percentage of tags, at each time interval, that were recaptured in the same location as tagged.
Metapopulation Dynamics In all cases, dispersal was higher in Leslie Tributary as compared with the Berczy Creek; however, this difference was only statistically significant in time periods 1, 5, 7 (Mann-Whitney U-test; p=0.015, 0.004, 0.033, respectively; Table 4.1). The difference in dispersal was not due
83
to more tags dispersing as there was no significant difference in the proportion of stationary tags, except for time period 4 (G-test, p=0.0003; Table 4.1). More likely, the increase in dispersal was due to larger average dispersal per individual, as indicated by mean and maximum dispersal through time (Table 4.1).
Metapopulation and Patch Viability Metapopulation viability, as indicated by both probability of extinction through time and the intrinsic mean time to extinction, were orders of magnitude different depending on whether the patch-occupancy model was parameterized indirectly using colonization potential or directly using observed colonization (Figure 4.2; Table 4.2). For example, when the patch-occupancy models were parameterized directly using observed colonization, the Leslie Tributary population was inviable long-term and the Berczy Creek showed much longer viability. The intrinsic mean time to extinction for Leslie Tributary was 24 years, and occurred in as little as 12 years (regional stochasicity set at 0, 0.2, respectively). The probability of extinction was plotted over time and showed that 95% of the simulations were extinct in less than 100 years. Similarly, the intrinsic mean time to extinction for Berczy Creek was calculated as a maximum of 148 years and occurred in as little as 32 years (regional stochasticity = 0, 0.2 respectively), with 95% of simulations showing metapopulation extinction in under 1000 years (Figure 4.2; Table 4.2).
84 Table 4.2 – Intrinsic mean time (in years) to extinction (Grimm and Wissel 2004) of two metapopulations of the endangered fish the redside dace (Clinostomus elongatus) for different levels of regional stochasticity.
B) Berczy Creek
A) Leslie Tributary Regional Stochasicity
0
0.1
0.2
0
0.1
0.2
Intrinsic mean time to extinction (Indirect Parameterization)
362
95
48
109,594
3,417
764
Intrinsic mean time to extinction (Direct Parameterization)
24
17
12
148
54
32
When the patch-occupancy models were parameterized using indirect colonization, based on estimated dispersal ability, one population (Berczy Creek) was deemed as viable and quasistationary (regional stochasticity = 0; Figure 4.2; Table 4.2), with 95% of simulations showing viability beyond 250,000 years (Figure 4.2). The remaining population estimates varied considerably in their viability, ranging in intrinsic mean times to extinction from 48 to 348 years in Leslie Tributary and from 764 to >109,000 years in Berczy Creek (regional stochasticity 0.2, 0, respectively).
85
I) Direct parameterization of dispersal
II) Indirect parameterization of dispersal
A)
B) YEARS Figure 4.2 – Metapopulation viability of the endangered species the redside dace (Clinostomus elongatus) in two stream metapopulations: A) Leslie Tributary, and B) Berczy Creek. Shown are the probabilities of extinction (y-axis) in years (x-axis) of a stochastic patch-based metapopulation model. Models were parameterized using: I) indirect parameterization of colonization and dispersal via patch distance, and; II) direct parameterization of colonization and dispersal using empirical estimates from a mark-recapture study. Legend: Vertical hashes represent a time interval of 100 years, solid lines indicate population trajectories where regional stochasticity was set to 0, dashed lines set to 0.1 and dotted lines set at 0.2.
86
Specific-patch viability mirrored overall metapopulation viability, with all patches showing reduced viability when parameterized directly from empirical data (Figure 4.3). In all cases (except L12) patch viability was over-estimated with the indirect parameterization relative to direct parameterization of the SPOM (Figure 4.3). Mean patch viability was significantly higher with indirect parameterization of the SPOM for both Leslie Tributary (mean indirect patch viability = 0.51 ± 0.10, mean direct patch viability = 0.35 ± 0.08; Welch’s t-test; t= 26.85, pvalue << 0.0001) and Berczy Creek (mean indirect patch viability = 0.69 ± 0.09, mean direct patch viability = 0.39 ± 0.08; Welch’s t-test; t = 11.204, p-value << 0.0001).
Patch Viability (Indirect Parameterization)
1 B18 L11 B10
B8
B9 B12
L8 L6
L10
B16 L9
L7
0.5
B11 B6
B7
B13
B15 L13 B14 L14
0
L15 B17
0
L12
0.5
1
Patch Viability (Direct Parameterization) Figure 4.3 – Differences in patch viability parameterized using indirect (y-axis) and direct (x-axis) patch dynamics of the endangered species the redside dace (Clinostomus elongatus) in two stream metapopulations: A) Leslie Tributary (L6-L15), and B) Berczy Creek (B6-B18). Shown are the mean probabilities of persistence of a given patch across 10,000 simulations. To demonstrate the variability in patch viability, 25% quantiles are overlaid as the negative of both the vertical and horizontal axes, while 75% quantiles are overlaid as the positive vertical and horizontal axes. The dashed line is a 1:1 line.
87
Interestingly, the rankings of patch viability did not markedly differ based on the parameterization of the SPOM. For example, both indirect and direct parameterization of the SPOM identified the same five most-viable patches per population (overall) as: L6, L8, L10, L11, L9; and, B7, B13, B6, B11, B18. One notable difference in patch viability was that several patches that were in close proximity to good-quality patches had significantly lower viability when parameterized using direct parameterization (Figure 4.3). Patches such as L9, L11, B6 and B9 had reduced viabilities when directly parameterized (mean indirect viability = 0.57, 0.84, 0.86, 0.90, respectively; mean direct viability = 0.28, 0.46, 0.63, 0.17; Figure 4.3).
Discussion Metapopulation-viability models have a long history of use (e.g. Levins 1969, 1970) and provide advantages over traditional population-viability analyses (PVAs). One clear advantage of using stochastic patch-occupancy metapopulation models (SPOMs) over traditional PVAs (e.g. structured models, demographic models; (Akcakaya and Sjorgen-Gulve 2000; Morris and Doak 2002) is that they require the parameterization of fewer variables (Ovaskainen and Hanski 2004). This reduction is especially advantageous for modeling endangered species, where enumeration is complicated by rarity and where greater uncertainty exists. Simplification is often needed as incorporating species-specific or demographic data into ecological studies can be difficult, time consuming, or not economically possible. Stochastic patch-occupancy models allow for simplification of metapopulations, as only patch occupancy, colonization and extinction rates are needed, even within a single snapshot (Moilanen et al. 2004; Marsh 2008). There is a tradeoff between simplification of PVAs by using less parameters and with the added value and information that those parameters may have (Shreeve et al. 2004). This study demonstrates that differences in species dispersal patterns, a key component of metapopulation models such as SPOMs, have the ability to dramatically impact estimates of metapopulation viability. Further, this study demonstrates that direct parameterization of species-specific dispersal can reduce the overall estimates of viability of redside dace metapopulations by several orders of magnitude over estimates using indirect estimates (e.g. exponential-decay kernels). In several cases (e.g. L9, L11, B6 and B9), patches that were in close proximity to good-quality patches (i.e. had high dispersal potential and high viability with indirect parameterization) but had significantly lower viability when parameterized using direct parameterization (Figure 4.3).
88
From these results it may be inferred that the study patches are likely distributed across a heterogeneous landscape where occupancy rates or habitats may not be equal, as assumed by SPOMs. Landscape and patch heterogeneity has been shown to impact colonization potential, which alters metapopulation dynamics (Gustafson and Gardner 1996; Heinz et al. 2005). Even if habitats were equal, studies of the rosyside dace (Clinostomus funduloides), a sister species, demonstrated that dispersal could not be predicted by suitable habitats (Freeman and Grossman 1993). These conclusions are in agreement with others which indicate that better integration of species-specific behaviour is needed into the analyses of metapopulations (Tischendorf 2001; Tischendorf and Fahrig 2001; Vos et al. 2001; Heinz et al. 2005, 2006; Baguette and vanDyck 2007; Marsh 2008). There is an ongoing debate whether SPOMs are useful in cases where the assumptions of classic (i.e. Levin’s type) metapopulations are not met (see Levins 1969; Harrison 1994; Baguette 2004; Hanski 2004; Shreeve et al. 2004). For example, empirical studies have demonstrated large temporal variation in patch dynamics (i.e. colonization and extinction rates) can lead to sensitivity in SPOMs (Crone et al. 2001; Thomas et al. 2002). In addition, there are few empirical examples of metapopulations that meet the assumption of a constant pulse of extinction-colonization (e.g. pool frog, Sjogren-Gulve 1991; Glanville fritillary butterfly, Hanski et al. 1994; but see Nowicki et al. 2007); whereas, the vast majority do not (Harrison 1994; Baguette et al 2004). However, recently SPOMs have been shown to be appropriate for use over a range of spatially structured populations from classic metapopulations to species found in fragmented landscapes with patchy distributions (Ovaskainen and Hanski 2004). For example, using a SPOM with individual-based background, Ovaskainen and Hanski (2004) demonstrated a unifying framework for incorporating metapopulation dynamics into SPOMs. Their study suggested that instead of attempting to identify metapopulation types, research should focus on relevant processes, such as dispersal. By understanding processes behind the variability shown between (and within) patch dynamics, the conservation and management of endangered species can be improved (Revilla and Wiegand 2008). Comparative simulations of how metapopulation analyses perform when altered are important tools for ensuring appropriate management action. Interestingly, despite finding large differences in the prediction of intrinsic mean time to extinction (Table 4.2), this study demonstrates that regardless of how the SPOMs were parameterized (directly or indirectly), they
89
identified qualitatively similar patches as having the highest viability (Figure 4.2). Such results are reassuring given that SPOMs rarely are used for exact quantitative analyses, but rather used to compare among several scenarios to develop decision support (Lindenmayer and Possingham 1996; Heinz et al. 2006). These results suggest that qualitative comparisons of SPOMs may still be a fruitful management option; however, care should be taken with estimates of probability of extinction through time, which may be over-estimated (as shown here). Indeed, others have shown SPOMs to be comparable to other spatially realistic models (Kindvall et al. 2000; Keeling 2002; Ovaskainen and Hanski 2004). Modifications to SPOMs, such as consideration of patch quality (Hanski and Moilanen 1998; Ovaskainen and Hanski 2002), improved dispersal metrics (Ovaskainen 2004; Heinz et al. 2005), incorporation of transition states (Thomas and Hanski 2004), rigorous parameter-estimation techniques (Moilanen 1999; Dreschler et al. 2003), and/or direct parameterization of dispersal data, will only ensure better integration of biological data into SPOMs. Studying metapopulations in stream settings is challenging and has limitations which should be noted. Fagan (2002) demonstrated how a dendritic network can provide additional isolation of patches not encountered in terrestrial landscapes. Further, Gotelli and Taylor (1999) demonstrated how stream fishes may not fit Levin’s type metapopulation models as migration may be asynchronous in upstream versus downstream movement. Finally, mark-recapture studies have shown that tagged individuals may travel outside recapture territory, thereby, reducing estimates of overall colonization (Ovaskainen 2004). All the above examples may explain the lower viability using direct parameterization of SPOMs in this study. In response, the study design accounted for such potential shortcomings in several ways. First, the sampling was done sequentially, across the entire stream network (using block nets during sampling to eliminate movement induced by the sampling of each pool). Therefore, unlike most studies in terrestrial systems, this study was able to monitor the stream metapopulations within defined boundaries, sampling the entire patch, as well as neighboring patches and connections. Second, this study re-sampled the patches at 7 time intervals to determine intra-and inter- annual asynchrony in movement and rate of patch fidelity. In general asynchrony in movement was identified at a single time period (T4; G-test, p<0.05); however, there was also variation in patch fidelity (Table 4.1). This variation in patch fidelity indicates that redside dace may be more consistent with the patchy population model (Harrison 1991) than the with a Levin’s type
90
metapopulation model (Levins 1969). Temporal variation in patch fidelity was also shown in metpopulations of rosyside dace (Freeman and Grossman 1993). Third, in this study the boundaries of the recapture locations were extended to include patches beyond the marking locations (~ 350 meters on upstream and downstream ends per location; Figure 4.1). These extended patches allowed both average dispersal of tagging individuals and the rates of dispersal out of their marking locations to be determined. Regardless of time period, > 75% of tags were recaptured within 350 m meters of their starting pool for Leslie Tributary and > 90% of tags in Berczy Creek, suggesting that the monitoring was done at the appropriate scale. Finally, this study accounted for the probability of missing tags due to sampling bias by correcting for abundance counts using probability of detection. These approaches suggest that, although asynchrony (and perhaps over-winter mortality) may be an issue, the study design was adequate to monitor species and patch specific movement. Understanding how species- and patch-specific qualities can alter SPOMs is an important area for advancing metapopulation-viability analyses, which are a crucial tool for the management of endangered species (Ackakaya 2000; Morris and Doak 2002; McCoy and Mushinsky 2007). Species- and patch-specific processes can alter metapopulation dynamics in many ways (Roitberg and Mangel 1997; Hanski and Moilanen 1998; Schtickzelle et al. 2006). Patch quality can affect both the probabilities of colonization and extinction of an empty patch (Hanski and Moilanen 1998; Thomas et al. 2002). Changes in landscape structure can alter migration pathways between patches, which may impact dispersal (Roitberg and Mangel 1997; Schtickzelle et al. 2006; North and Ovaskainen 2007). Finally, the configuration of patches may play a role in population viability (Robert 2009). Maintenance not only of high-quality patches, but their connectivity, should be an important aspect of endangered species management, including redside dace. Accounting for regional and stochastic processes are important considerations for the management of endangered species. In this study, when the rate of regional stochasticity was altered, the intrinisic mean time to extinction quickly decreased (Table 4.2). The impact of regional stochasticity was small in cases where populations were already considered to be inviable (e.g. Leslie Tributary with 0 stochasticity), but it had a large impact in cases where a population was viable or quasi-stationary (e.g. Berczy Creek with 0 stochasticity; Figure 4.2). For example, the Berczy Creek population was considered to be viable if regional stochasticity
91
was ignored. Altering rates of regional stochasticity from 0 to 0.1 and 0.2 caused the populations to move towards extinction (intrinsic mean time to extinction from >109,000 to 54 and 32 years, respectively). These results indicate that regional stochastic factors (e.g. weather, drought) that may alter patch dynamics and have undue influence on metapopulation viability, as others also have shown (Grossman et al. 1982; 1985; Lande 1993; Foley 1994; Robert 2009). The conservation applications of species-specific dispersal may help inform management decisions. By incorporating species- and patch-specific data directly into metapopulation models, managers may be better apt at determining the relative importance of spatial and temporal factors, such patch connectivity and seasonal impacts. In this study, empirical estimation of patch viability was shown to be qualitatively similar when parameterizing SPOMs, but that estimates of metapopulation viability were shown to be significantly higher when using indirect parameterization of dispersal. These results indicate that care is needed in ensuring that even simplified metapopulation models, such as SPOMs, are consistent with biological data. Comparisons of how species- and patch-specific data directly impact metapopulation models, as done here, may be one way to accomplish this cautionary aspect (Dreschler et al. 2003; Grimm et al. 2004). Further study into the impact of species-specific behavior and patch dynamics on metapopulation viability will provide additional insight for both management and conservation issues.
92
Acknowledgements Funding was provided by NSERC Canada and OGS Scholarships to M.S.P., an NSERC Discovery Grant to D.A.J., Interdepartmental Recovery Fund #1410 provided by Fisheries and Oceans (DFO), the Ontario Ministry of Natural Resources (OMNR), and the University of Toronto. All field work was conducted under an approved Animal Care Protocol (# 20006805) from the University of Toronto Animal Care Committee; and under the approval of the redside dace recovery team (chair A. Dextrase) and with the guidance of the Ontario Ministry of Natural Resources (J. Pisapio). The Toronto Region Conservation Authority (TRCA) Aquatics Group (Christine Tu, David Lawrie, Tim Rance) provided logistical support and help during this project. Field work was conducted by a dedicated group of volunteers from- the University of Toronto – A. Drake, C. Harpur, B. Edwards, M. St. John, M. Neff, P. Venturelli, A. Manning, M. Granados, J. Ruppert, S. Sharma, N. Puckett, C. Hart, C. Howard, M. Luksenberg, J. Brett, S. Walker – and the Toronto Region Conservation Authority – D. Lawrie, C. Tu, T. Parker, T. Rance, B. Paul, E. Elton, B. Stephens, C. Hart, L. DelGiudice, B. Moyle, and M. Parish – and P. Ng, M. Ken, C. Zehr, B. Edwards, Y. Nozoe, D. Trim, K. Lee, D. Morodvanschi and D. Forder (Ontario Streams). This work was improved by discussion with Marie-Josee Fortin and review of earlier drafts of this manuscript. Finally reviewers were helpful for in providing valuable suggestions.
References Akçakaya, H. R. 2000. Viability analyses with habitat-based metapopulation models Population Ecology 42:45-53. Akçakaya, H. R., and P. Sjögren-Gulve. 2000. Population viability analysis in conservation planning: an overview. Ecological Bulletins 48:9-21. Baguette, M. 2004. The classical metapopulation theory and the real, natural world: a critical appraisal. Basic and Applied Ecology 5:213-224. Baguette, M., S. Petit, and F. Queva. 2000. Population spatial structure and migration of three butterfly species within the same habitat network: consequences for conservation. Journal of Animal Ecology 37:100-108.
93
Baguette, M., and H. VanDyck. 2007. Landscape connectivity and animal behavior: functional grain as a key determinant for dispersal. Landscape Ecology 22:1117-1129. Carle, F. L., and M. R. Strub. 1978. A new method for estimating population size from removal data. Biometrics 34:621-630. COSEWIC. 2007. COSEWIC Assessment and update status report on the redside dace Clinostomus elongatus in Canada. Committee on the Status of Endangered Wildlife in Canada, Ottawa, ON. Crone, E. E., D. Doak, and J. Pokki. 2001. Ecological influences on the dynamics of a field vole metapopulation. Ecology 82:831-843. Dewey, M. R., and S. J. Ziegler. 1996. An evaluation of flourescent elastomer for marking bluegills in experimental studies. Progressive Fish-Culturist 58:219-220. Dreschler, M., K. Frank, I. Hanski, R. B. O’Hara, and C. Wissel. 2003. Ranking metapopulation extinction risk: From patterns in data to conservation management decisions. Ecological Applications 13:990-998. Fagan, W. F. 2002. Connectivity, fragmentation, and extinction risk in dendritic metapopulations. Ecology 83:3243-3249. Fahrig, L. 1992. Relative importance of spatial and temporal scales in a patchy environment. Theoretical Population Biology 41:300-314. Foley, P. 1994. Predicting extinction times from environmental stochasticity and carrying capacity. Conservation Biology 8:124-137. Frank, K., and C. Wissel. 1998. Spatial aspects of metapopulation survival: from model results to rules of thumb for landscape management. Landscape Ecology 13:363-379. Frank, K., and C. Wissel. 2002. A formula for the mean lifetime of metapopulations in heterogeneous landscapes. The American Naturalist 159:530-552.
94
Freeman, M. C., and G. D. Grossman. 1993. Effects of habitat availability on dispersion of a stream cyprinid Environmental Biology of Fishes 37:121-130. Ghimire, S. K., O. Gimenez, R. Pradel, D. McKey, and Y. Aumeeruddy-Thomas. 2008. Demographic variation and population viability in a threatened Himalayan medicinal and aromatic herb Nardostachys grandiflora: matrix modelling of harvesting effects in two contrasting habitats. Journal of Applied Ecology 45:41-51. Goldsmith, R. J., G. P. Closs, and H. Steen. 2003. Evaluation of visible implant elastomer for individual marking of small perch and common bully. Journal of Fish Biology 63:631636. Gotelli, N. J., and C. M. Taylor. 1999. Testing macroecology models with stream-fish assemblages. Evolutionary Ecology Research 1:847-858. Grimm, V., and I. Storch. 2000. Minimum viable population size of capercaillie Terao urogallus: Results from a stochastic model. Wildlife Biology 5:219-225. Grimm, V., and C. Wissel. 2004. The intrinsic mean time to extinction: a unifying approach to analysing persistence and viability of populations. Oikos 105:501-511. Grimm, V., H. Lorek, J. Finke, F. Koester, M. Malachinski, M. Sonnenschein, A. Moilanen, I. Storch, A. Singer, C. Wissel, and K. Frank. 2004. META-X: generic software for metapopulation viability analysis. Biodiversity and Conservation 13:165-188. Grossman, G. D., M. C. Freeman, P. B. Moyle, and J. O. Whittakers. 1985. Stochasticity and assemblage organization in an Indian stream fish assemblage. American Naturalist 126:275-285. Grossman, G. D., P. B. Moyle, and J. O. Whittakers. 1982. Stochasticity in structural and functional characteristics of an Indiana stream fish assemblage: A test of community theory. American Naturalist 120:423-454. Gustafson, E. J., and R. H. Gardner. 1996. The effect of landscape heterogeneity on the probability of patch colonization. Ecology 77:94-107.
95
Hanski, I. 1994. A practical model of metapopulation dynamics. The Journal of Animal Ecology 63:151-162. Hanski, I. 1998. Connecting the parameters of local extinction and metapopulation dynamics. Oikos 83:390-396. Hanski, I. 1999a. Metapopulation Ecology. Oxford University Press, New York. Hanski, I. 1999b. Habitat connectivity, habitat continuity, and metapopulations in dynamic landscapes. Oikos 87:209-219. Hanski, I. 2004. Metapopulation theory, its use and misuse. Basic and Applied Ecology 5:225229. Hanski, I., and M. E. Gilpin. 1997. Metapopulation biology: Ecology, genetics, and evolution. Academic Press, San Diego, California. Hanski, I., M. Kuussaari, and M. Nieminen. 1994. Metapopulation structure and migration in the butterfly Melitaea Cinxia. Ecology 75:747-762. Hanski, I., A. Moilanen, T. Pakkala, and M. Kuussaari. 1996. The quantitative incidence function model and persistence of an endangered butterfly metapopulation. Conservation Biololgy 10:578-590. Hansson, L. 1991. Dispersal and connectivity in metapopulations. Biological Journal of the Linnean Society 42:89-103. Harrison, S. 1991. Local extinction in a metapopulation context. Biological Journal of the Linnean Society 42:73-888. Heinz, S. K., L. Conradt, C. Wissel, and K. Frank. 2005. Dispersal in fragmented landscapes: Deriving a practical formula for patch accessibility. Landscape Ecology 20:83-99. Heinz, S. K., C. Wissel, and K. Frank. 2006. The viability of metapopulations: individual dispersal behaviour matters. Landscape Ecology 21:77-89.
96
Hill, J. K., C. D. Thomas, and O. T. Lewis. 1996. Effects of habitat patch size and isolation on dispersal by Hesperia comma butterflies: implications for metapopulation structure. Journal of Animal Ecology 65:725-735. Hokit, D. G., B. M. Stith, and L. C. Branch. 1999. Effects of landscape structure in Florida scrub: a population perspective. Ecological Applications 9:124-134. Keeling, M. J. 2002. Using individual-based simulations to test the Levins metapopulation paradigm. Journal of Animal Ecology 71:270-279. Kindvall, O. 2000. Comparative precision of three spatially realistic simulation models of metapopulation dynamics. Ecological Bulletins 48:101-110. Lande, R. 1993. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. The American Naturalist 142:911. Lande, R., S. Engen, and B. E. Saether. 1988. Extinction times in finite metapopulation models with stochastic local dynamics. Oikos 83:383-389. Levins, R. 1969. Some demographic and genetic consequences of environmental heterogeneity for biological control. Bulletin of the Entomological Society of America 71:237-240. Marsh, D. 2008. Metapopulation viability analysis for amphibians. Animal Conservation 11:463465. McCoy, E. G., and H. R. Mushinsky. 2007. Estimates of minimum patch size depend on the method of estimation and the condition of habitat. Ecology 88:1401-1407. Moilanen, A. 1999. Patch occupancy models of metapopulation dynamics: Efficient parameter estimation using implicit statistical inference. Ecology 80:1031-1043. Moilanen, A. 2004. SPOMSIM: software for stochastic patch occupancy models of metapopulation dynamics. Ecological Modelling 179:533-550. Moilanen, A., and I. Hanski. 1998. Metapopulation dynamics: effects of habitat quality and landscape structure. Ecology 79:2503-2515.
97
Moilanen, A., A. T. Smith, and I. Hanski. 1998. Long-term dynamics in a metapopulation of the American pika. The American Naturalist 152:530-542. Morris, W. F., and D. F. Doak. 2002. Quantitative conservation biology: Theory and practice of population viability analysis. Sinauer Associates Inc., Sutherland, MA. Nowicki, P., A. Pepkowska, J. Kudlek, P. Skorka, M. Witek, J. Settele, and M. Woyciechowski. 2007. From metapopulation theory to conservation recommendations: Lessons from spatial occurrence and abundance patterns of Maculinea butterflies. Biological Conservation 140:119-129. North, A., and O. Ovaskainen. 2007. Interactions between dispersal, competition and landscape heterogeneity. Oikos 116:1106-1119. Ogle, D.H. 2009. Functions to support fish stock assessment textbook (v.0.0-8). Accessed online (January 12, 2009): http://www.ncfaculty.net/dogle/ Oksanen, J. 2004. Incidence function model in R. Accessed online (February 12, 2009): http://cc.oulu.fi/~jarioksa/opetus/openmeta/metafit.pdf OMNR (Ontario Ministry of Natural Resources). 2007. Stream assessment protocol for southern Ontario. Ontario Ministry of Natural Resources, Picton, Ontario. Ovaskainen, O. 2004. Habitat-specific movement parameters estimated using mark-recapture data and a diffusion model. Ecology 85:242-257. Ovaskainen, O., and I. Hanski. 2002. Transient dynamics in metapopulation response to perturbation. Theoretical Population Biology 61:285-295. Ovaskainen, O., and I. Hanski. 2004. From individual behavior to metapopulation dynamics: Unifying the patchy population and classic metapopulation models. American Naturalist 164:364-377. Poos, M.S., and D.A. Jackson. 2009. Conservation by consensus: Reducing uncertainties in modeling the distribution of an endangered species using habitat-based ensemble models. Ecological Applications (Submitted 2009-06-08: #09-1012).
98
R Development Core Team. 2008. A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. (ISBN: 3-900051-07-0). Reid, S. M., N. E. Jones, and G. Yunker. 2008. Evaluation of single-pass electrofishing and rapid habitat assessment for monitoring redside dace. North American Journal of Fisheries Management 28:50-56. Resetarits, W. J., C. A. Binckley, and D. R. Chalcraftet. 2005. Habitat selection, species interactions, and processes of community assembly in complex landscapes: a metacommunity perspective. Pages 374-398 in M. Holyoak, M. A. Leibold, and R. D. Holt, editors. Metacommunities - spatial dynamics and ecological communities. The University of Chicago Pressz, Chicago, Illinois. Revilla, E., and T. Wiegand. 2008. Movement ecology special feature: Individual movement behavior, matrix heterogeneity, and the dynamics of spatially structured populations. Proceedings of the National Academy of Science 105:19120-19125. Robert, A. 2009. The effects of spatially correlated perturbations and habitat configuration on metapopulation persistence. Oikos (Online Early, doi: 10.1111/j.16000706.2009.17818.x). Roitberg, B. D., and M. Mangel. 1997. Individuals on the landscape: Behavior can mitigate landscape differences among habitats. Oikos 80:234-240. Schtickzelle, N., G. Mennechez, and M. Baguette. 2006 Dispersal depression with habitat fragmentation in the bog fritillary butterfly. Ecology 87:1057-1065. Shreeve, T. G., R. L. H. Dennis, and H. VanDyck. 2004. Resources, habitats and metapopulations - whither reality? Oikos 106:404-408. Sjogren-Gulve, P. 1991. Extinctions and isolation gradients in metapopulations: the case of the pool frog (Rana lessonae). Biological Journal of the Linnean Society 42:135-147. Thomas, C. D., and I. Hanski. 1997. Butterfly metapopulations. Pages 359-386 in I. Hanski and M. E. Gilpin, editors. Metapopulation Biology: Ecology, Genetics and Evolution. Academic Press, London.
99
Thomas, C. D., R. J. Wilson, and O. T. Lewis. 2002. Short-term studies underestimate 30generation changes in a butterfly metapopulation. Proceedings of the Royal Society of London B 269:563-569. Tischendorf, L. 2001. Can landscape indices predict ecological processes consistently? Landscape Ecology 16:235-254. Tischendorf, L., and L. Fahrig. 2001. On the use of connectivity measures in spatial ecology. A reply. Oikos 95:152-155. Vandewoestijne, S., T. Martin, S. Liégeois, and M. Baguette. 2004. Dispersal, landscape occupancy and population structure in the butterfly Melanargia galathea. Basic and Applied Ecology 5:581-591. Verboom, J., J. A. J. Metz, and E. Meelis. 1993. Metapopulation models for impact assessment of fragmentation. Pages 172-192 in C. C. Vos and P. Opdam, editors. Landscape ecology of a stressed environment. Chapman and Hall, London. Verboom, J., K. Lankester, and J. A. J. Metz. 1991. Linking local and regional dynamics in stochastic metapopulation models. Biological Journal of the Linnean Society 42:39-55. Vos, C. C., J. Verboom, P. F. M. Opdam, and C. J. F. ter Braak. 2001. Toward ecologically scaled landscape indices. American Naturalist 157:24-41. Walsh, M. G., and D. L. Winkelman. 2004. Anchor and visible implant elastomer tag retention by hatchery rainbow trout stocked into an Ozark stream. North AMerican Journal of Aquaculture 24:1435-1439. Zar, J.H. 1999. Biostatistical analysis. Upper Saddle River, New Jersey: Prentice Hall.
100
Appendices
Appendix 4.1 –The endangered redside dace (Clinostomus elongatus). Photo credit: Mark Poos.
Appendix 4.2 –Visual implant elastomer (VIE) tag inserted subcutaneously on the ventral surface of the endangered redside dace (Clinostomus elongatus). Photo credit: Mark Poos.
101 Appendix 4.3 – Model parameters used in the stochastic patch-occupancy models. Shown are the number of emigrants (mi), the mean probability of detection (PD), the number of emigrants adjusted for potentially missed tags (Mi), the number of immigrants needed to start a new sub-population, and the rate of extinction (Ei). Other parameters include: the incidence function for Leslie Tributary (dI = 210, x = 0.4926, e = 3.685, y = 6.12), and Berczy Creek (dI = 150, x = 0.5652, e = 4.187, y = 7.01). Not shown: dij given it is a pairwise estimate rather than unique for each pool.
Patch
mi
PD
Mi
Ij, 0.5
Ei
Leslie Tributary L6 L7 L8 L9 L10 L11 L12 L13 L14 L15 Overall
35 7 38 10 37 38 14 8 0 3 190
0.710 0.619 0.767 0.740 0.647 0.780 0.727 0.792 0.000 0.800 0.718
60 11 67 17 61 68 24 14 0 3 326
4.89 2.48 5.48 8.21 5.10 5.74 8.03 2.88 7.69 7.69
0.2593 0.2901 0.1956 0.2724 0.2583 0.1870 0.3167 0.2912 0.3383 0.3304
Berczy Creek B6 B7 B8 B9 B10 B11 B12 B13 B14 B15 B16 B17 B18 Overall
48 116 23 8 17 34 9 88 2 5 8 0 49 407
0.676 0.564 0.582 0.667 0.836 0.741 0.800 0.715 0.882 0.800 0.744 0.667 0.661 0.656
80 181 36 13 31 59 16 151 4 9 14 0 81 670
1.91 7.43 2.78 9.20 2.85 4.91 7.43 5.91 10.94 10.94 1.91 7.43 2.78
0.2645 0.1918 0.2231 0.1394 0.1639 0.1683 0.2026 0.2025 0.2110 0.2726 0.1569 0.2388 0.1417
102
Section II: Reducing uncertainty from methodological choices using consensus methods
103
Chapter 5: Reducing uncertainties in modeling the distribution of endangered species using habitatbased ensemble models
Abstract Modeling the distribution of endangered species is often problematic due to their rarity and the consequential statistical issues. Recent studies have demonstrated the utility of using ensemble models to reduce the uncertainty of singular predictions when modeling species distributions. This study provides a quantitative evaluation of the efficacy of ensemble models for improving the prediction of endangered species and their habitats using the endangered fish, the redside dace (Clinostomus elongatus), as a model organism. Specifically this study asks: 1) how well do ensemble models improve modeling metrics (e.g. specificity, sensitivity and overall classification); 2) how many singular methods are needed to build the optimal ensemble; and, 3) what scale(s) and type(s) of habitat are most important for modeling this species. For this evaluation, five ensemble models were compared based on seven singular approaches. Habitat variables were derived from 200 sites measured from 1997-2007 and divided hierarchically into fine-, intermediate-, and broad-scale habitat. In all cases, the ensemble models were equal to or better than any singular method across all modeling metrics, although there was large variation with certain combinations of initial models. This study demonstrates how comparative analysis of modeling types, ensemble approaches and scales can be useful for reducing uncertainty in the modeling of endangered species and their habitats. Keywords: ensemble model, consensus model, conservation, biodiversity, biogeography, species distributions.
104
Introduction Predictive models of species distributions have become important tools in modeling changes in biogeography and biodiversity (Guisan & Zimmermann, 2000; Guisan & Thuiller, 2005; Sharma and Jackson 2008). Studies have used predictive models in a wide range of applications including to evaluate species distributions in relation to climate change (Thuiller 2004; Araújo et al. 2005), evaluate the establishment and spread of expanding invasive species (Hartley et al. 2006), and model habitat suitability of endangered species (Rodriguez et al. 2007; Marmion et al. 2009; Franklin et al. 2009). With the increasing availability of remote sensing data and advances of geographic information systems, researchers often are no longer data limited, and are expanding the use of predictive models to include more applications (Guisan and Zimmerman 2000; Marmion et al. 2009). Modeling existing and future habitats of endangered species has become a popular application of predictive models. One reason for this popularity is that faced with limited data on the distribution, abundance and dynamics of endangered species (Mace et al., 2005; Rodríguez, 2007), predictive models allow for the extrapolation of relatively few field samples to the entire potential range of a species. However, the application of predictive models to endangered species remains controversial (Loiselle 2003; Wilson et al. 2005; Thompson et al. 2007). One difficulty with the standard sets of predictive models used by the majority of ecologists, is that they are often inappropriate when analyzing data limited to a few sites and scales (Ellison and Agrawal 2005; Araújo & Guisán, 2006). This, in turn, produces data sets that have many statistical limitations, including zero-inflated bias, increased collinearity between variables, and inflated coefficient of variation (Graham 2003; Dixon et al. 2005; Edwards et al. 2006; Guisan et al. 2006; Dormann et al. 2008). However, uncertainties in predictive models of endangered species may arise during all stages of modeling including obtaining species level data, obtaining accurate species counts, and, linking species to landscapes/habitats (Elith et al. 2002; Loiselle 2003; Heikkinen et al. 2006). Therefore, approaches that reduce uncertainty in predictive models are continually being sought, as are improvements to current approaches (Elith et al. 2006; Austin 2007). One way to overcome the difficulties of using singular predictive models has been the use of ensemble-based approaches, which combine several predictive models (Aruajo and New 2007;
105
Marimon et al. 2008). The goal of ensemble models is to reduce the uncertainty in any singular approach by combining their predictions. The idea of ensemble approaches dates back to Laplace (1820) who demonstrated that the probability of error will rapidly decline with the inclusion of additional predictors. However, in comparison to other disciplines, the application of ensemble models to ecological issues remains in its infancy (Bates and Granger 1969; Arujo et al. 2005b). In particular, the use of ensemble models have been restricted to determining distributions of species under climate-change scenarios (Thuiller 2004; Araujo et al. 2005b; Thuiller et al 2005; Buisson and Grenouillet 2008), with less emphasis on establishing habitat relationships. However, as the first step of building an ensemble model is developing singular predictive models of initial conditions as filters, studies regarding uncertainties with the existing singular approaches from which the ensembles are based remain relatively understudied. The use of ensemble models for improving model predictions of species distributions and their habitats has received little attention. Comparative studies of efficiencies of ensemble models are scarce (Thuiller 2004; Marimon et al. 2008), and few studies provide guidance on how to build the most suitable ensemble (e.g. how do models improve prediction, how many initial filters are needed?). Given the paucity of both ensemble models and data with most endangered species (Mace et al. 2005; Rodriguez et al. 2007), an evaluation is needed of the application of ensemble models for identifying habitat characteristics of endangered species. The objective of this paper is to provide a quantitative comparative analysis of both singular predictive models and ensemble approaches for identifying existing habitat for the endangered species, the redside dace (Clinostomus elongatus). Specifically, this study seeks to determine: 1) whether habitat-based ensemble models improve modeling metrics (e.g. sensitivity, specificity, and overall classification); 2) to what degree ensemble models actually improve predictive success; and, 3) how many models should be used to build an optimal ensemble?
106
Methods Study Area and Species The redside dace is a pool-dwelling minnow found only in Ontario within its Canadian range (Figure 5.1). Redside dace is in decline due to changes in land-use through urbanization, especially in the Greater Toronto Area, which includes 80% of its Canadian range, (COSEWIC 2007). To evaluate the efficacy of predictive models for correctly classifying redside dace locations and their habitats, a dataset of 100 known redside dace locations sampled from 19972007 was compiled from historical records from various government agencies, universities and conservation authorities. Absence data were obtained from a larger dataset of 560 locations obtained from the same agencies and randomly reduced to 100 locations as a means of developing a balanced approach for assessing modeling metrics (Olden et al. 2004). A hierarchical approach was used to predict the presence and absence of redside dace. In total 28 variables were used across multiple scales; fine-scale (i.e. site level habitat), intermediatescale (i.e. landscape level) and broad-scale (i.e. geologic and geomorphic variables) to help quantify the importance of redside dace habitat (Table 5.1). These variables included factors thought to influence redside dace, such as urban land cover (Scott and Crossman 1973; McKee and Parker 1982; Daniels and Wisniewski 1994; Novinger and Coon 2000; COSEWIC 2007), ranging to factors thought to influence stream fish in general (Grossman and Freeman 1987). Fine-scale habitat features were derived using the Ontario Stream Assessment Protocol (OMNR 2007) and included: habitat type (e.g. percent pool, riffle, run); substrate type (e.g. percent cobble, gravel or fines); stream depth (m); type of in-stream cover (e.g. percent flat rock, round rock, wood, bank, or macrophytes); characteristics of the stream bank (e.g. percent eroding, vulnerable, protected or depositional); overland temperature (degrees Celsius), and amount of adjacent riparian cover (Table 5.1). Intermediate-scale habitat variables included percent coverage by urban, forest, cropland, pasture and wetland landcover and were obtained from most recent satellite imagery (2001) that was converted in shapefiles using a geographic information system. Various spatial buffer sizes were used to characterize adjacent land-use; however a one kilometer buffer was deemed most appropriate given a sensitivity analysis using various buffer sizes. Broad-scale habitat variables were assessed using various buffer sizes in a similar fashion
107
to landscape variables. The broad-scale variables were derived from The Canadian System of Soil Classification and included type of soil geology in the vicinity (Newmarket Till, Elma Till, Halton Till), as well as slope (Soil Landscapes of Canada Working Group 2007). Soil categories that occurred at < 5% of sites were removed in order to not over-parameterize the models.
Lake Huron
Greater Toronto Area
Ontario, CANADA Lake Ontario
New York, 78°20’5 U.S.A.
Lake Erie
42°41’2
Figure 5.1 – Distribution of sampling locations between 1997-2007. Closed circles indicate redside dace occurrences, whereas, open circles indicate where redside dace were absent.
108 Table 5.1 – Summary of the hierarchical habitat-based model used for predicting the presence of endangered minnow, redside dace (Clinostomus elongatus). Seven predictive models were used including: LR (logistic regression), CT (classification trees), MARS (multivariate adaptative regression splines), ANN (artificial neural networks), DA (discriminant analysis), RF (random forest), and BR (boosted regression trees). Variables were derived using forward-selection procedures on five independent datasets and are shown as a percentage of datasets where each variable was selected (parentheses indicate negative associations). Variables selected from only one dataset were omitted. A prioi predictions (ap) based on habitat predictors thought to influence the decline of the species are shown for reference, where + indicates a positive correlation, - negative correlation, 0 none.
Scale Fine (Site)
Intermediate (Landscape)
Broad (Geologic)
Variable
ap
LR
Width Pool Glides Fast Riffles Slow Riffles Fine substrate Gravel substrate Cobble substrate Shallow depth Intermediate depth Deep Flat rock cover Round rock cover Wood cover Bank cover Macrophyte cover Eroding banks Vulnerable banks Protected banks Depositional banks Urban Cropland Pasture Forest Wetland Newmarket Till Halton Till Elma Till Slope Temperature
+ 0 0 0 0 0 0 0 0 0 0 0 0 + 0 + 0 0 0 0 0 0 0 0 0 0 0 0
[100] 100
CT 60
MARS ANN [80] 100
DA
RF
BR
[100] 80
[80] 100
100
100
[80] [80]
[60]
[80]
[40] 40
80
80
40 40
[80] [100]
[60]
[100]
[100]
[100]
80
40
80
60
[100]
[80] 80 40
109
Building Individual Models Previous comparisons of individual models have demonstrated that model comparisons are needed for understanding the relative strengths and weakness of each modeling approach, and this is especially true for modeling endangered species where greater uncertainty exists (Olden and Jackson 2001; 2002a). In this stud, five ensemble methods were compared for improving the predictive success of the redside dace and its habitat. For this assessment, seven modeling approaches were used as initial filters to identify environmental variables linked with the occurrence of redside dace (Figure 5.2). These modeling approaches were: logistic regression (LR); classification trees (CT); multivariate adaptive regression splines (MARS); artificial neural networks (ANN); discriminant analysis (DA); random forest (RF); and, boosted regression trees (BR). Each of these modeling approaches has been shown to be useful for describing species occurrences (Thuiller 2003; Elith et al. 2006; Thuiller et al. 2006; Marimon et al. 2008) and were used to develop the output models (Figure 5.2). The modeling methods used represent a continuum of use from logistic regression, the most prevalent and widespread statistical method for modeling binary data (Hosmer and Lemeshow 1989; Pampel 2000), to newer methodologies that may perform better, such as artificial neural networks, boosted regression trees and random forest (Olden and Jackson 2001; Olden et al. 2004; Thuiller et al. 2004). In general, the modeling methods can be broadly broken down into three groups: two regression-based methods (LR, MARS); three machine-learning methods (BR, RF, ANN); and, two classification methods (CT and DA; Elith et al. 2006; Marrimon et al. 2009).
110
Step 1: Partition Data
A123 B456 C789
Step 2: Analyze Training Data Logistic Regression (LG)
Multivariate Classification Adaptive Trees Regression (CT) Splines (MARS)
Artificial Discriminant Neural Analysis Networks (DA) (ANN)
Random Forest (RF)
Boosted Regression (BR)
Step 3: Produce Output Models Using Predictions of Independent Test Data
A=1 B=0 C=1
A=1 B=1 C=0
A=1 B=0 C=0
Step 4: Build Ensemble Forecast (e.g. Mean, PCA, WA)
A=0 B=0 C=1
A=1 B=1 C=0
A=0 B=0 C=0
A=0 B=1 C=1
A=1 B=0 C=0
Figure 5.2 – Summary of ensemble forecasting approach. Five-fold cross validation was used by partitioning data into five individual test and validation datasets. Using each training dataset, data was analyzed individually across seven initial approaches: logistic regression (LR), classification trees (CT), multivariate adaptive regression splines (MARS), artificial neural networks (ANN), discriminant analysis (DA), random forest (RF) and boosted regression (BR). Five ensemble models were built from initial seven predictions, including: a consensus model (CM), principal components analysis (PCA), weighted average using overall classification (WA), mean (Mn) and median (Md). Modeling metrics (specificity, sensitivity and overall classification) were obtained by comparing the individual or ensemble predictions on the independent test dataset.
111
Models were configured in the following ways. LRs were run using a logit function using maximum likelihood (Allison 1999; Olden and Jackson 2002a). CTs used a chi-square distance of 0.05 to determine significant group cut-offs. The number of terminal nodes in the leaf structure were optimized by iteratively running all terminal nodes and choosing the leaf structure with the smallest error (in this case a terminal structure of 5 was determined to have the lowest error). This step was important as too many terminal nodes may artificially inflate correct classification, whereas too few terminal nodes may not provide meaningful species groupings (Vayssiéres et al. 2000). These models were developed in the SAS programming language v.9.1. MARS provides an alternative regression based method using piecewise linear fits rather than a smoothing parameter (e.g. as in general additive models; Friedman 1991; Elith et al. 2006). ANNs used a one-hidden layer feed-forward network trained by the back-propagation algorithm (Bishop 1995). This type of network is considered a universal approximator of any continuous function, has low associated rates of error and is used most often in ecological studies (Rumelhart et al. 1986; Hornik and White 1989; Olden and Jackson 2002b, Olden et al. 2004). ANNs were optimized (optimal referring to minimizing the trade-off between network bias and variance) for the number of hidden neurons in the neural network by determining empirically the number of hidden neurons that produces the lowest misclassification rate (Bishop 1995); which in this case was a layer that contained seven nodes. Random forests were run for each node of the tree, randomly using m variables on which to base the decision at that node. The best split was calculated based on these m variables in the training set. Two group (presence/absence) discriminant function analysis was developed using a linear model. Models 3-6 were built using the statistical program Statistica v.7. BR prediction was based on an accuracy-weighted vote across estimated classifiers (Ridgeway 1999) and run using code provided by Elith et al. (2008) in the R programming language v.2.8.0 (R Development Core Team, 2008) Habitat variables were chosen using forward step-wise selection methods; except for the machine-learning methods (ANN, RF, BR), which iteratively fit their own responses to the habitat data (Olden and Jackson 2002b; Dormann et al. 2008). This approach was necessary to compare the ability of individual models for selecting habitat variables with high explanatory power. Although bias can be introduced by variable selection methods, recent studies have demonstrated that selection procedures may actually have very small impact on resultant models (Maggini et al. 2006; Meynard and Quinn 2007; Dormann et al. 2008). A forward-selection
112
procedure (p < 0.05) was used to determine variables with high explanatory power. In the cases of machine-learning methods, the single ANN model with the highest overall classification (e.g. highest area under the curve) was retained. With the RF and BR approaches the importance of variables was determined by estimating the relative influence of each variable reducing the loss function, based on a square error algorithm (Thuiller, 2003). Model thresholds were calculated using receiver operator characteristic curves (ROC) with thresholds balanced between misclassification of species presence and species absence (Olden and Jackson 2001). It should be noted that the choice of model thresholds can impact the results (Jiménez-Valverde & Lobo 2007; Lobo et al. 2008). For example, researchers must weigh the need for correct classification of either presence or absence at the expense of inflated error (i.e. misclassification). One could choose model-selection thresholds that prioritize correct classification, at the expense of increased misclassification of species absences (or vice versa). Such a decision may be well suited for conservation applications, where one would wish to emphasize correct predictions of the true occurrence of an endangered species, but were willing to accept higher commission rates (Loiselle et al. 2003). For the purpose of this study, a balanced approach was used, where models were select that had equal likelihood of misclassification of species presence and absence. This approach was well suited for this study as model were evaluated based on their specificity, sensitivity and overall classification, which require equal consideration (Hartley et al. 2006). In addition, previous studies (e.g. Marrimon et al. 2008), which have used area-under-the-curve operations to evaluate predictive performance of ensemble models may inappropriately produce models that have high overall classification due solely to model specificity, which would not otherwise be known (Loiselle et al 2003; Jiménez-Valverde. & Lobo 2007; Lobo et al. 2008). For all of the model comparisons, data were split into two types: 80% of the data were used to train each model, and 20% of the data were used to validate and test the predictive capability of each model (Figure 5.2). A 5-fold cross validation procedure was used where five separate models were calculated for each statistical method to build a complete validation data set (Figure 5.2). The validation data is an important component of model evaluation, as it provides an independent unbiased evaluation of the predictive capability of each trained model (Olden et al. 2002a).
113
Evaluating Individual Models Predictive models were compared in a number of ways. First, model outputs were compared to a set of a prioi predictions based on habitat predictors thought to influence the decline of the species (Scott and Crossman 1973; McKee and Parker 1982; Daniels and Wisniewski 1994; Novinger and Coon 2000). Second, predictive models were evaluated by three modeling metrics: model sensitivity (the ability of each model to correctly predict species presence); model specificity (the ability of each model to correctly predict species absence); and, overall classification (the ability of each model to correctly classify both species presence and absence). The use of these metrics provides an alternative means to evaluate each model, their comparative successes, their associated errors, and their relationship to ecological relevant variables, for both species presence and absence (Olden and Jackson 2002a). The distinction between modeling metrics is an important one as models for predicting imperiled species are often difficult to assess because they tend to artificially inflate correct classification of species absence (model specificity), where the majority of sites are associated and the habitats do not reflect the ecology of the model species, whereas, model sensitivity measures the ability to classify species presence, which may be more relevant for conservation applications. Finally, predictive outputs were compared to the observed data using unweighted pair-group method (UPGMA) cluster analysis using phi similarity. Cluster analysis has been used widely in ecological literature (Legendre and Legendre 1997) as the resultant output produces a dendrogram which connects closely matching neighbors. For example, statistical methods connected by branches proximal to one another match closer than methods connected by branches more distal (Podani 2000). The phi coefficient is the binary correlation coefficient and is not biased by frequency of occurrence as has been shown for other coefficients (Jackson et al 1989).
Building Ensemble Models Ensemble models allow for a decrease in predictive uncertainty of singular models, by using a combination of their predictions (Figure 5.2). Various ensemble models have been proposed, including those that use selective algorithms or basic mathematical functions such as the mean and median predictions (Araújo et al. 2005; Marrimon et al. 2008). Here, five ensemble approaches currently in favor were compared: a consensus model (CM), Principal Components Analysis (PCA), weighted average using overall classification (WA), Mean (Mn) and Median (Md) ensemble approaches. Each of these approaches have been used extensively in building
114
ensemble models (e.g. Gregory et al., 2001; Johnson and Omland 2004; Thuiller, 2004; Araújo et al 2005a; 2005b, 2006; Thuiller et al 2005; Araújo and New 2007; Goswami and O’Connor, 2007; Marimon et al. 2008) All ensemble models were built from reducing a matrix of model predictions (j) by sites (i) to a vector of ensemble predictions. For example, the consensus ensemble model used a majority rules criterion, where all model outputs were compared and the majority output (presence or absence) was retained. For the mean ensemble model, all continuous model outputs were averaged and the final output was rounded to either presence or absence. Similarly, the median ensemble model took the median value of all seven individual outputs (either presence or absence). The principal components analysis used the dominant axis of an eigenvalue decomposition of a covariance matrix programmed in the R v.2.8.0 (R Development Core Team, 2008). The resultant eigenvalue was scaled to presence absence using the ROC approach highlighted earlier. Finally, the weighted average ensemble model used the overall classification (across all predictions) and multiplied it to each prediction to produce a singular ensemble prediction.
Evaluating Ensemble Models Ensemble models were evaluated using identical metrics as the individual models: specificity, sensitivity and overall classification, and cluster analysis. The number of initial models needed to build the most appropriate ensemble model was also evaluated. For this, a re-sampling approach was developed using a sample-based rarefaction routine (Colwell et al. 2004), which was coded using R v.2.8.0 (R Development Core Team, 2008). This re-sampling approach randomly selected predictions at a given site from a matrix of model predictions by site. Each prediction was matched to the observed data and model specificity, sensitivity and overall classification was calculated across study sites. The set of predictions was then randomly permuted and a new prediction was added from a competing model and evaluated as to whether this additional prediction increased overall classification (Figure 5.2). This randomization procedure was repeated 1000 times to determine 95% confidence intervals for all combinations of ensemble models.
115
Results I) Individual Models The majority of the initial predictive models identified similar fine-scale habitat features, including: positive associations with deep, pool habitats with gravel substrate, and negative associations with stream width and fast riffles (Table 5.1). Only a negative association with urban land cover and a positive association with Elma till, came out as a strong indicator at a broader scale. These findings are in agreement with previous studies of current knowledge of this species and its habitat (Scott and Crossman 1973; McKee and Parker 1982; Daniels and Wisniewski 1994; Novinger and Coon 2000; COSEWIC 2007; Table 5.1). All individual models generally performed well (> 80% overall classification). Performance was generally better for model specificity (range 84-90%; Table 5.2) than for model sensitivity (7587%). When individual models were compared, logistic regression was most closely related to the observed data, followed by the learning-based approaches, such as boosted regression, random forest, and artificial neural networks (Fig. 5.3A). These methods also had the highest rates of model specificity and generally the highest rates of model sensitivity, although artificial neural networks dropped off in this area (Table 5.2). Alternatively, classification-based methods, such as classification trees and discriminant analysis were least associated with the observed data, followed by multivariate adaptive regression splines.
116 Table 5.2 – A comparison of model sensitivity (correct classification of species presence), specificity (correct classification of species absence), and overall classification (correct classification of both species presence and absence) for redside dace (Clinostomus elongatus). Single models are: LR (logistic regression); CT (classification trees); BR (boosted regression trees); MARS (multivariate adaptive regression splines); ANN (artificial neural networks); DA (discriminant analysis); and, RF (random forest). Ensemble forecasts are: consensus model (CM); principal component analysis (PCA); weighted average using overall classification (WA); mean (Mn); and, median (Md).
Singular Models
Ensemble Models
LG
CT
BR
MARS
ANN
DA
RF
CM
Sensitivity (%)
87
75
83
81
75
84
83
Specificity (%)
90
86
90
88
89
84
Overall Classification (%)
88
81
87
85
82
84
PCA
WA
Mn
Md
87
88
90
88
87
89
94
91
90
89
92
86
91
90
90
89
90
117
II) Ensemble Models Ensemble models provided as good as or improved model performance over singular methods. Overall classification increased from 1-10% with ensemble models (Table 5.2). On average, ensemble models improved model sensitivity over model specificity. For example, whereas, only logistic regression had model sensitivity above 85%, all ensemble models had rates of model sensitivity at 87% or above and had model sensitivity equal or above to all the individual models. With the exception of the mean ensemble result, all ensemble models provided equal or superior specificity of 89% rather than the 90% achieved by the best individual model. Overall, the various ensemble methods produced similar levels of model sensitivity, model specificity and overall classification; however, the consensus ensemble model provided the highest rate of classification, due to having both higher values of model sensitivity and specificity. The various ensemble models show close association (i.e. phi correlation measure) to one another (Fig. 5.3) and to the observed data, which were similar to the association between the observed data and the best individual model (i.e. logistic regression).
Figure 5.3- Cluster analysis showing the relationship with the observed distribution (Obs.) of redside dace (Clinostomus elongatus) and: A) the seven individual modeling approaches alone, and B) with ensemble forecasts included. Model short forms are carried over from Table 5.2.
118
There was large variation when comparing all possible combinations of consensus ensemble models. Whereas, model sensitivity and overall classification increased (on average) with the addition of more initial models, the average model specificity decreased for combinations of three and five initial models, before improving at the final seven-input ensemble model (Fig. 5.4). In addition, the variation (e.g. 25th and 75th percentiles, as well as minima and maxima) in model sensitivity, specificity and overall classification increased when larger combinations of ensembles were considered (e.g. there were 35 combinations of ensembles for three-input models, and 21 for five-ensemble models). For the three-model consensus, the combination of LG, MARS, RF had the highest model specificity, CT, ANN, BR had the highest model sensitivity and LG, RF, BR had the highest overall classification. Alternatively, the combination of learning-based methods ANN, RF, BR had the lowest model specificity and overall classification, whereas, the combination of LG, CT, RF showed the lowest model sensitivity (Appendix 5.1). Similarly, for the five-model consensus models, combinations with all three learning-based methods performed the worst, where the combination of MARS, ANN, RF, BR, DA had the lowest sensitivity and overall classification and LG, MARS, ANN, RF, BR had the lowest specificity (Appendix 5.1).
119
A)
B)
C) Figure 5.4 – Box and whisker plots showing the variability in consensus ensemble forecasts for predicting the presence of an endangered redside dace. Consensus ensemble models were compared across all combinations of one (n =7), three (n =35), five (n =21) and seven (n=1) input models (x-axis). Boxes are 25th and 75th percentiles, horizontal lines indicate the median, vertical lines indicate the upper and lower values, diamonds indicate the mean and are connected by dashed lines. Modeling metrics were: A) model sensitivity (i.e. correct classification of species presence); B) model specificity (i.e. correct classification of species absence); and, C) overall classification. Dashed lines indicate 95% confidence intervals.
120
Discussion Ensemble models provide a useful method for reducing uncertainty when modeling distributions of endangered species. Endangered species are not only rare, but they are also often difficult to capture and enumerate, thereby, complicating the evaluation of their habitat. Thus methods which can be used to reduce the uncertainty of model type should allow for improvement in the management of endangered species. In this study, ensemble models improved model specificity for the endangered redside dace (Clinostomus elongatus) by 6.9%, model sensitivity by 3.2% and overall classification by 5.3% (on average); over individual models. Previous studies have demonstrated that ensemble approaches can reduce uncertainty and improve model fit for predicting future distributions (Thuiller, 2004; Thuiller et al., 2005; Araujo and New 2007; Marrimon et al. 2008). Here, this study provides the first evidence that ensemble models can be equally applied to modeling current distribution and provides increased model performance over singular approaches. Ensemble models may be especially useful for identifying the potential importance of habitat for endangered species, and at the scale in which they function. Modeling the distributions of endangered species and their habitats is an activity filled with uncertainty. Using the endangered redside dace as an example, this study demonstrates that although predictive models varied in their ability to identify correctly existing habitat, in all cases an ensemble approach improved model prediction. In addition, whereas previous studies have demonstrated the utility of ensemble models for identifying habitats at broader scales (e.g. landscape), this study demonstrates the utility in developing a multi-scaled approach. With the increase in availability of remote-sensing data and advances in geographic information systems, researchers are left with a plethora of data from which to model distributions (Guisan and Zimmerman 2000). Often climatic, topographic and landuse data are available readily (Guralnick et al. 2007), whereas information on fine-scale variables is more scarce (Austin 2007; Dorman et al. 2008). This study demonstrates a situation where fine-scale habitat models had a greater ability to describe factors relevant to the redside dace than broader scale factors (Table 5.1). The need for comparative approaches using many statistical methods has also been highlighted as another way of reducing uncertainty with modeling endangered species (Guisan and Zimmerman 2000; Olden and Jackson 2002a). The choice of modeling approach has been
121
shown to have severe consequences for the application to conservation decisions (Loiselle 2003; Wilson et al. 2005; Rodriguez et al. 2007; Marrimon et al. 2008). For example, Pearson et al. (2006) showed that distribution changes of South African plant species varied from 92% loss to a 322% gain depending on the model they used. In addition, Dorman et al (2008) demonstrated that of several uncertainties in modeling species distributions, including variable selection and collinearity between variables, the choice of model type had the largest impact. In this study, there was variation in how each model type fitted predictions, with logistic regression most closely resembled the observed distribution, followed closely by the three machine learning methods: boosted regression, random forest and artificial neural networks. These results are in general agreement with other studies showing the benefits of machine-learning methods such as random forest, and general boosted regression (Cutler et al., 2007). One reason for the superior performance of logistic regression is that more complex models (which iteratively fit a solution) can be prone to overfitting (Olden and Jackson 2002a). As previous quantitative comparative analyses have demonstrated, the success of predictive modeling approaches can be largely data dependent and there is no clear indication of the preeminence of any singular approach (Olden and Jackson 2002a; Araujo and New 2007). Comparative analyses can be used to identify problems related to modeling approaches given the available data and to determine situations where comparative analysis may work better (Guisan and Zimmerman 2000; Olden and Jackson 2002a). Comparative approaches are needed not just to ensure that appropriate models are being developed, but to ensure the best ensemble models are being built. Comparisons of ensemble models are rare (e.g. Araujo et al. 2005; Marrimon et al. 2008) and only one (Marimon et al. 2008) has evaluated the relationship between ensemble models and predictive performance. In this study, there was variation in the ability of ensemble models to improve prediction over singular approaches. The consensus ensemble model (i.e. vote counting) provided the highest model specificity and the best estimated overall classification, whereas the weighted average ensemble model provided the highest model sensitivity (Table 5.2). The remaining ensemble models performed similarly, including those based on principal components analysis, which have been used most frequently in previous studies (Thuiller 2004; Araujo et al. 2005b). The reason for this similarity is likely due to the lack of independence in the information provided by the initial models. Future research is needed into understanding the study settings in which different ensemble methods are likely to perform best (Marrimon et al. 2007).
122
The evaluation of the efficacy of combining initial models into ensemble models has received little attention and remains an important knowledge gap (Loiselle et al. 2003; Wilson et al. 2005; Marrimon et al. 2008). Decisions such as how many initial models are needed to build the most useful ensemble model have not been addressed. Ensemble models work under the assumption that the initial models provide some independent information (Araujo and New 2007): at some stage, adding more initial models beyond those that provide independent information may actually decrease the utility of the ensemble approach. Here a simple sub-sampling procedure (Colwell et al. 2004) was used to demonstrate that an optimal ensemble model could be developed using as few as three initial models (Figure 5.3). As understanding the contributions of initial models to ensemble forecasts remains a limitation to the ensemble approach (Araujo and New 2007), the use of sample based rarefaction techniques may allow substantial insight into ensuring that appropriate ensemble models are being produced. The application of ensemble models to ecological issues remains in its infancy (Araujo et al. 2005b). Previous studies on ensemble models have demonstrated their utility in reducing uncertainty with singular methods (Thuiller et al. 2004; Araujo et al. 2005b; Marrimon et al. 2007), however several challenges remain to ensure ensemble models are applied appropriately, especially for modeling distributions of endangered species. First, researchers should have thorough understanding of predictive modeling, their uncertainties and conditions under which they should be applied (Elith et al. 2002; Loiselle et al. 2003). Second, researchers should be mindful of the limitations of working with endangered species, which provide additional statistical problems (Thompson 2003; Ellison and Agrawal 2005). Third, comparative analyses of not only singular methods, but ensemble methods are needed to ensure that the most appropriate models are being constructed. The development of appropriate validation data (Araujo and New 2007), decision of modeling metrics (e.g. AUC, sensitivity, specificity), and consideration of how many initial models are needed to build the ensemble, are all needed to ensure that ensemble models are actually improvements to singular methods. Fourth, the goals of the ensemble model need to be explicitly testable. Such testing should include the use of models to corroborate hypothesis with either a priori prediction or existing knowledge and ensuring that habitat is modeled using appropriate scale/s.
123
Conclusion This study demonstrates that ensemble models can provide marked improvements to singular approaches in situations where there is uncertainty exists, such as modeling suitable habitat of endangered species. Using the endangered redside dace as a model organism, this study shows that an ensemble approach can improve model sensitivity, specificity and overall classification. These findings have important consequences for improving species distribution models of endangered species. In addition, this study demonstrates the utility of sub-sampling procedures for determining how many initial models are needed to build an optimal ensemble. As ensemble models are dependent on their initial inputs, improvements to future ensemble models can be expected with greater consideration of how initial models perform best. This study provides an example of how comparative analyses across many scales, initial modeling types and ensemble approaches can be used to improve the prediction of endangered species and their habitats.
Acknowledgements Funding was provided by NSERC Canada and OGS Scholarships to M.S.P., an NSERC Discovery Grant to D.A.J., Interdepartmental Recovery Fund #1410 provided by Fisheries and Oceans (DFO), the Ontario Ministry of Natural Resources (OMNR), and the University of Toronto. The following agencies and individuals were helpful for provided data for this research: Royal Ontario Museum, Fisheries and Oceans, Ontario Ministry of Natural Resources, Toronto Region Conservation Authority (TRCA), Conservation Halton (CH), Lower Lake Simcoe Conservation Authority (LLSCA), Credit Valley Conservation Authority, D Forder (Ontario Streams), J Anderson (LLSCA), L Stanfield (OMNR), E Holm (ROM), D Lawrie (TRCA), S Jarvie (TRCA), S Watson-Leung (CH), and S Reid (DFO). In addition, M Neff and G Rawnsley provided field assistance. This manuscript benefited from discussions with C Harpur and A Drake. Finally anonymous reviewers were helpful for providing comments on early drafts of this paper.
124
References Allison, P. D. 1999. Logistic Regression Using the SAS System. Cary, N.C. SAS Institute. Araujo, M. B., and M. New. 2007. Ensemble forecasting of species distributions. TRENDS in Ecology and Evolution 22:42-47. Araújo, M. B., R. G. Pearson, W. Thuiller, and M. Erhard. 2005a. Validation of species–climate impact models under climate change Global Change Biology 11:1504-1513. Araújo, M. B., R. J. Whittaker, R. J. Ladle, and M. Erhard. 2005b. Reducing uncertainty in projections of extinction risk from climate change. Global Ecology and Biogeography 14:529-538. Araújo, M. B., and A. Guisán. 2006. Five (or so) challenges for species distribution modelling. Journal of Biogeography 33:1677-1688. Austin, M. 2007. Species distribution models and ecological theory: a critical assessment and some possible new approaches. Ecological Modelling 200:1-19. Bishop, C. M. 1995 Neural Networks for Pattern Recognition. Oxford University Press, New York, U.S.A. Burgman, M., D. B. Lindenmayer, and J. Eltih. 2005. Managing landscapes for conservation under uncertainty. Ecology 86:2007-2017. Buisson, L., and G. Grenouillet. 2009. Contrasted impacts of climate change on stream fish assemblages along an environmental gradient. Diversity and Distributions 15:613-626. Burnham, K. P., and D. R. Anderson. 2002. Model selection and multimodel inference: a practical information-theoretic approach, 2nd edition. Springer-Verlag, New York, New York. Cade, B. S., B. R. Noon, and C. H. Flather. 2005. Quantile regression reveals hidden bias and uncertainty in habitat models. Ecology 86:786-800.
125
Colwell, R. K., C. X. Mao, and J. Chang. 2004. Interpolating, extrapolating, and comparing incidence-based species accumulation curves. Ecology 85:2717-2727. COSEWIC. 2007. COSEWIC Assessment and update status report on the redside dace Clinostomus elongatus in Canada. Committee on the Status of Endangered Wildlife in Canada, Ottawa, ON. Culter, D. R., T. C. Edwards, K. H. Beard, A. Cutler, K. T. Hess, J. Gibson, and J. J. Lawler. 2007. Random forests for classification in ecology. Ecology 88:2783-2792. Cunningham, R. B., and D. B. Lindenmayer. 2005. Modeling count data of rare species: Some statistical issues. Ecology 86:1135-1142. Daniels, R. A., and S. J. Wisniewski. 1994. Feeding ecology of redside dace, Clinostomus elongatus. Ecology of Freshwater Fish 3:176-183. Dixon, P. M., A. M. Ellison, and N. J. Gotelli. 2005. Improving the precision of estimates of the frequency of rare events. Ecology 86:1114–1123. Dormann, C. F., O. Purschke, J. R. G. Márquez, S. Lautenbach, and B. Schröder. 2008. Components of uncertainty in species distribution analysis: A case study of the great grey shrike. Ecology 89:3371-3386. Edwards, T. C., D. R. Cutler, N. E. Zimmermann, L. Geiser, and J. Alegriae. 2005. Model-based stratifications for enhancing the detection of rare ecological events. Ecology 86:10811090. Elith, J., M. A. Burgman, and H. M. Regan. 2002. Mapping epistemic uncertainties and vague concepts in predictions of species distributions. Ecological Modelling 157:313-329. Elith, J., C. H. Graham, R. P. Anderson, M. Dudík, S. Ferrier, A. Guisan, R. J. Hijmans, F. Huettmann, J. R. Leathwick, A. Lehmann, J. Li, L. G. Lohmann, B. A. Loiselle, G. Manion, C. Moritz, M. Nakamura, Y. Nakazawa, J. M. Overton, A. T. Peterson, S. J. Phillips, K. Richardson, R. Scachetti-Pereira, R. E. Schapire, J. Soberón, S. Williams, M. S. Wisz, and N. E. Zimmermann. 2006. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29:129-151.
126
Elith, J., J. R. Leathwick, and T. Hastie. 2008. A working guide to boosted regression trees. Journal of Animal Ecology 77:802-813. Ellison, A. M., and A. A. Agrawal. 2005. The statistics of rarity. Ecology 86:1079-1080. Goswami, M., and K. M. O’Connor. 2007. Real-time flow forecasting in the absence of quantitative precipitation forecasts: a multi-model approach. Journal of Hydrology 334:125-140. Graham, M. H. 2003. Confronting multicollinearity in ecological multiple regression. Ecology 84:2809-2815. Gregory, A. W., G. W. Smith, and J. Yetman. 2001. Testing for forecast consensus. Journal of Business and Economic Statistics 19:34-43. Grossman, G. D., and M. C. Freeman. 1987. Microhabitat use in a stream fish assemblage. Journal of Zoology 212:151-176. Guisan, A., O. Broennimann, R. Engler, M. Vust, N. G. Yoccoz, A. Lehmann, and N. E. Zimmermann. 2006. Using niche-based distribution models to improve the sampling of rare species. Conservation Biology 20:501-511 Guisan, A., and W. Thuiller. 2005. Predicting species distribution: offering more than simple habitat models. Ecology Letters 8:993-1009. Guisan, A., and N. E. Zimmermann. 2000 Predictive habitat distribution models in ecology. Ecological Modelling 135:147-186. Guralnick, R. P., A. W. Hill, and M. Lane. 2007. Towards a collaborative, global infrastructure for biodiversity assessment. Ecology Letters 10:663-672. Hartley, S., R. Harris, and P. J. Lester. 2006. Quantifying uncertainty in the potential distribution of an invasive species: climate and the Argentine ant. Ecology Letters 9:1068-1079.
127
Heikkinen, R. K., M. Luoto, M. B. Araújo, R. Virkkala, W. Thuiller, and M. T. Sykes. 2006. Methods and uncertainties in bioclimatic envelope modelling under climate change. Progress in Physical Geography 30:751-777. Hornik, K., and H. White. 1989. Multilayer feedforward networks are universal approximators. Neural Networks 2:359-366. Hosmer, D. W., and S. Lemeshow. 1989. Applied Logistic Regression. John Wiley and Sons, New York. Jiménez-Valverde, A., and J. M. Lobo. 2007. Threshold criteria for conversion of probability of species presence to either–or presence–absence. Acta Oecologica 31:361-369. Johnson, J. B., and K. S. Omland. 2004. Model selection in ecology and evolution. TRENDS in Ecology and Evolution 19:101-108. Laplace, P. S. 1820. Théorie analytique des probabilités. Courcier, Paris. Legendre, P., and L. Legendre. 1998. Numerical Ecology, 2nd edition. Elsevier Science BV, Amsterdam, Neitherlands (853 pp.). Lobo, J. M., A. Jiménez-Valverde, and R. Real. 2008. AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography 17:145-151. Loiselle, B. A., C. A. Howell, C. H. Graham, J. M. Goerck, T. Brooks, K. G. Smith, and P. H. Williams. 2003. Avoiding pitfalls of using species distribution models in conservation planning. Conservation Biology 17:1591-1600. Mace, G. M., H. P. Possingham, and N. Leader-Williams. 2005. Prioritizing choices in conservation. Pages 17-34 In D. MacDonald and K. Services (ed.). Key Topics in Conservation Biology. Blackwell Publishers, Oxford. MacNally, R. 2000. Regression and model-building in conservation biology, biogeography and ecology: The distinction between - and reconciliation of – ‘predictive’ and ‘explanatory’ models. Biodiversity and Conservation 9:655-671.
128
Maggini, R., A. Lehmann, N. E. Zimmermann, and A. Guisan. 2006. Improving generalized regression analysis for the spatial prediction of forest communities. Journal of Biogeography 33:1729-1749. Marmion, M., M. Parviainen, M. Luoto, R. K. Heikkinen, and W. Thuiller. 2009. Evaluation of consensus methods in predictive species distribution modelling. Diversity and Distributions 15: 59-69. McKee, P. M., and B. J. Parker. 1982. The distribution biology and status of the fishes Campostoma anomalum Clinostomus elongatus Notropis photogenis Cyprinidae and Fundulus notatus Cyprinodontidae in Canada. Canadian Journal of Zoology 60:13471358. McCullagh, P., and J. A. Nelder. 1989. Generalized Linear Models. Chapman and Hall, London, England. Meynard, C. N., and J. F. Quinn. 2007. Predicting species distributions: a critical comparison of the most common statistical models using artificial species. Journal of Biogeography 34:1455-1469. Moisen, G. G., E. A. Freeman, J. A. Blackard, T. S. Frescino, N. E. Zimmermann, and T. C. Edwards. 2006. Predicting tree species presence and basal area in Utah: A comparison of stochastic gradient boosting, generalized additive models, and tree-based methods. Ecological Modelling 199: 176-187. Novinger, D. C., and T. G. Coon. 2000. Behavior and physiology of the redside dace, Clinostomus elongatus, a threatened species in Michigan. Environmental Biology of Fishes 57:315-326. Olden, J. D., and D. A. Jackson. 2001. Fish habitat relationships in lakes: gaining predictive and explanatory insight by using artificial neural networks. Transactions of the American Fisheries Society 130:878-897. Olden, J. D., and D. A. Jackson. 2002a. A comparison of statistical approaches for modelling fish species distributions. Freshwater Biology 47:1976-1995.
129
Olden, J. D., and D. A. Jackson. 2002b. Illuminating the 'black box': Understanding variable contributions in artificial neural networks. Ecological Modelling 154:135-150. Olden, J. D., M. K. Joy, and R. G. Death. 2004. An accurate comparison of methods for quantifying variable importance in artificial neural networks using simulated data. Ecological Modelling 178:389-397. OMNR. 2007. Stream assessment protocol for southern Ontario. Ontario Ministry of Natural Resources, Picton, Ontario. Pampel, F. C. 2000. Logistic regression: A primer. Sage Publications, London. Pearson, R. G., W. Thuiller, M. B. Araújo, E. Martinez-Meyer, L. Brotons, C. McClean, L. Miles, P. Segurado, T. P. Dawson, and D. C. Lees. 2006. Model-based uncertainty in species’ range prediction. Journal of Biogeography 33:1704-1711. Podani, J. 2000. Introduction to the exploration of multivariate biological data. Backhuys Publishers, Leiden, The Netherlands. R Development Core Team. 2008. A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. (ISBN: 3-900051-07-0). Ridgeway, G. 1999. The state of boosting. Computing Sciences and Statistics 31:172-181. Rodríguez, J. P., L. Brotons, J. Bustamante, and J. Seoane. 2007. The application of predictive modelling of species distribution to biodiversity conservation. Diversity and Distributions 13:243–251. Rumelhart, D. E., G. E. Hinton, and R. J. Williams. 1986. Learning representations by backpropagation errors. Nature 323:533-536. Scott, W. B., and E. J. Crossman. 1973. Freshwater fishes of Canada. Bulletin 184 Fisheries Research Board of Canada, Ottawa, ON.
130
Sharma, S. and D.A. Jackson. 2008. Predicting smallmouth bass incidence across North America: Comparison of statistical approaches. Canadian Journal of Fisheries and Aquatic Sciences 65: 471-481. Soil Landscapes of Canada Working Group, 2007. Soil Landscapes of Canada v3.1.1 Agriculture and Agri-Food Canada. (digital map and database at 1:1 million scale). Thomson, J. R., E. Fleishman, R. MacNally, and D. S. Dobkin. 2007. Comparison of predictor sets for species richness and the number of rare species of butterflies and birds. Journal of Biogeography 34:90-101. Thuiller, W. 2003. BIOMOD: optimizing predictions of species distributions and projecting potential future shifts under global change. Global Change Biology 9:1353-1362. Thuiller, W. 2004. Patterns and uncertainties of species’ range shifts under climate change. Global Change Biology 10:2020-2027. Thuiller, W., S. Lavorel, M. B. Araújo, M. T. Sykes, and I. C. Prentice. 2005. Climate change threats to plant diversity in Europe. Proceedings of the National Academy of Sciences 102:8245-8250. Thuiller, W. F., G. Midgley, M. Rougeti, and R. Cowling. 2006. Predicting patterns of plant species richness in megadiverse South Africa. Ecography 29:733-744. Vayssiéres, M. P., R. E. Plant, and B. H. Allen-Diaz. 2000. Classification Trees: An Alternative Non-Parametric Approach for Predicting Species Distributions. Journal of Vegetation Science 11:679-694. Whittingham, M. J., P. A. Stephens, R. B. Bradbury, and R. P. Freckleton. 2006. Why do we still use stepwise modelling in ecology and behaviour? Journal of Animal Ecology 75:11821189. Wilson, K. A., M. I. Westphal, H. P. Possingham, and J. Elith. 2005. Sensitivity of conservation planning to different approaches to using predicted species distribution data. Biological Conservation 122:99-112.
131
Appendices Appendix 5.1 – Model metrics for all combinations of consensus ensemble models. Models are: LR (logistic regression), CT (classification trees), MARS (multivariate adaptive regression splines), RF (random forest), ANN (artificial neural networks), BR (boosted regression trees, and DA (discriminant analysis). Overall Model/s Included Sensitivity Specificity Classification LG 0.8700 0.9000 0.8800 CT 0.7500 0.8600 0.8100 MARS 0.8100 0.8800 0.8500 RF 0.8300 0.8900 0.8600 ANN 0.7500 0.8900 0.8200 BR 0.8300 0.9000 0.8700 DA 0.8400 0.8400 0.8400 LG,CT,MARS 0.8822 0.8516 0.8669 LG,CT,ANN 0.8428 0.8650 0.8539 LG,CT,RF 0.8878 0.7922 0.8400 LG,CT,BR 0.8784 0.8254 0.8519 LG,CT,DA 0.9350 0.8198 0.8774 LG,MARS,ANN 0.8640 0.8796 0.8718 LG,MARS,RF 0.9318 0.8188 0.8749 LG,MARS,BR 0.9222 0.8250 0.8736 LG,MARS,DA 0.9088 0.8228 0.8658 LG,ANN,RF 0.8222 0.8646 0.8434 LG,ANN,BR 0.7900 0.8730 0.8315 LG,ANN,DA 0.8594 0.8396 0.8495 LG,RF,BR 0.9164 0.8656 0.8909 LG,RF,DA 0.9286 0.8102 0.8694 LG,BR,DA 0.9118 0.8610 0.8864 CT,MARS,ANN 0.9058 0.8532 0.8795 CT,MARS,RF 0.9150 0.8292 0.8721 CT,MARS,BR 0.8836 0.8408 0.8622 CT,MARS,DA 0.9296 0.8400 0.8848 CT,ANN,RF 0.9056 0.8626 0.8841 CT,ANN,BR 0.8688 0.8894 0.8791 CT,ANN,DA 0.8980 0.8804 0.8892 CT,RF,BR 0.8670 0.8280 0.8475 CT,RF,DA 0.9174 0.8106 0.8640 CT,BR,DA 0.8970 0.8372 0.8671 MARS,ANN,RF 0.7368 0.8718 0.8043
132
MARS,ANN,BR MARS,ANN,DA MARS,RF,BR MARS,RF,DA MARS,BR,DA ANN,RF,BR ANN,RF,DA ANN,BR,DA RF,BR,DA LG,CT,MARS,ANN,RF LG,CT,MARS,ANN,BR LG,CT,MARS,ANN,DA LG,CT,MARS,RF,BR LG,CT,MARS,RF,DA LG,CT,MARS,BR,DA LG,CT,ANN,RF,BR LG,CT,ANN,RF,DA LG,CT,ANN,BR,DA LG,CT,RF,BR,DA LG,MARS,ANN,RF,BR LG,MARS,ANN,RF,DA LG,MARS,ANN,BR,DA LG,MARS,RF,BR,DA LG,ANN,RF,BR,DA CT,MARS,ANN,RF,BR CT,MARS,ANN,RF,DA CT,MARS,ANN,BR,DA CT,MARS,RF,BR,DA CT,ANN,RF,BR,DA MARS,ANN,RF,BR,DA
0.7182 0.7430 0.7808 0.7982 0.8052 0.6836 0.7234 0.7018 0.7818 0.8744 0.8654 0.8624 0.9064 0.9194 0.8908 0.8574 0.8656 0.8432 0.8888 0.8232 0.8792 0.8238 0.9072 0.8094 0.8762 0.9140 0.8804 0.9078 0.8504 0.7368
0.8446 0.8284 0.8068 0.8018 0.8274 0.8766 0.8670 0.8750 0.8424 0.8826 0.8892 0.8558 0.8502 0.8344 0.8730 0.9074 0.8660 0.9006 0.8256 0.8244 0.8366 0.8686 0.8344 0.8720 0.8834 0.8562 0.8944 0.8462 0.8802 0.8424
0.7814 0.7857 0.7938 0.8000 0.8163 0.7801 0.7952 0.7884 0.8121 0.8785 0.8773 0.8591 0.8783 0.8769 0.8819 0.8824 0.8658 0.8719 0.8572 0.8238 0.8529 0.8462 0.8708 0.8407 0.8798 0.8851 0.8874 0.8770 0.8653 0.7896
133
Appendix 5.2 – R Code for testing configurations of 1,3, 5 and 7 prediction consensus models data<-read.table ("F:\\R\\Predict.txt",h=T) iters<-10000
# load prediction matrix # define number of iterations
model<-matrix(nrow=iters,ncol=8) for (i in 1:iters){
# start loop
(data[1:200,1:8][sample(200,1,replace=TRUE),])->temp 1 row (with replacement)
# from sample of 200, randomly pick
as.matrix(temp)->temp2 temp2->model[i,] number of random sites }
# make temporary matrix # fill this matrix with model outputs from n
model[1:iters,1]->OBS model[1:iters,2]->LG model[1:iters,3]->CT model[1:iters,4]->MARS model[1:iters,5]->RF model[1:iters,6]->ANN model[1:iters,7]->BR model[1:iters,8]->DA
# define where obs is (e.g. col 2) # name model in column 3 # name nmodel in column 4 # name model in column 5 # name model in column 6 # name model in column 7 # name model in column 8 # name model in column 9
# define combination # THREE MODEL COMBINATIONS combo<-cbind(LG,CT,GLM) #combo<-cbind(LG,CT,ANN) #combo<-cbind(LG,CT,RF) #combo<-cbind(LG,CT,BR) #combo<-cbind(LG,CT,DA) #combo<-cbind(LG,MARS,ANN) #combo<-cbind(LG,MARS,RF) #combo<-cbind(LG,MARS,BR) #combo<-cbind(LG,MARS,DA) #combo<-cbind(LG,ANN,RF) #combo<-cbind(LG,ANN,BR) #combo<-cbind(LG,ANN,DA) #combo<-cbind(LG,RF,BR) #combo<-cbind(LG,RF,DA) #combo<-cbind(LG,BR,DA)
# end loop
134
#combo<-cbind(CT,MARS,ANN) #combo<-cbind(CT,MARS,RF) #combo<-cbind(CT,MARS,BR) #combo<-cbind(CT,MARS,DA) #combo<-cbind(CT,ANN,RF) #combo<-cbind(CT,ANN,BR) #combo<-cbind(CT,ANN,DA) #combo<-cbind(CT,RF,BR) #combo<-cbind(CT,RF,DA) #combo<-cbind(CT,BR,DA) #combo<-cbind(MARS,ANN,RF) #combo<-cbind(MARS,ANN,BR) #combo<-cbind(MARS,ANN,DA) #combo<-cbind(MARS,RF,BR) #combo<-cbind(MARS,RF,DA) #combo<-cbind(MARS,BR,DA) #combo<-cbind(ANN,RF,BR) #combo<-cbind(ANN,RF,DA) #combo<-cbind(ANN,BR,DA) #combo<-cbind(RF,BR,DA) # FIVE MODEL COMBINATIONS #combo<-cbind(LG,CT,MARS,ANN,RF) #combo<-cbind(LG,CT,MARS,ANN,BR) #combo<-cbind(LG,CT,MARS,ANN,DA) #combo<-cbind(LG,CT,MARS,RF,BR) #combo<-cbind(LG,CT,MARS,RF,DA) #combo<-cbind(LG,CT,MARS,BR,DA) #combo<-cbind(LG,CT,ANN,RF,BR) #combo<-cbind(LG,CT,ANN,RF,DA) #combo<-cbind(LG,CT,ANN,BR,DA) #combo<-cbind(LG,CT,RF,BR,DA) #combo<-cbind(LG,MARS,ANN,RF,BR) #combo<-cbind(LG,MARS,ANN,RF,DA) #combo<-cbind(LG,MARS,ANN,BR,DA) #combo<-cbind(LG,MARS,RF,BR,DA) #combo<-cbind(LG,ANN,RF,BR,DA) #combo<-cbind(CT,MARS,ANN,RF,BR) #combo<-cbind(CT,MARS,ANN,RF,DA) #combo<-cbind(CT,MARS,ANN,BR,DA) #combo<-cbind(CT,MARS,RF,BR,DA) #combo<-cbind(CT,ANN,RF,BR,DA) #combo<-cbind(MARS,ANN,RF,BR,DA) # SEVEN MODEL COMBINATION
135
#combo<-cbind(LG,CT,MARS,RF,ANN,BR,DA) combo1<-apply(combo,1,mean) combined (e.g. combo#) as.matrix(combo1)->t matrix round(t)->RoundAll matrix(nrow=iters,ncol=1)->OC classificaiton (OC) matrix(nrow=iters,ncol=1)->presence matrix(nrow=iters,ncol=1)->absence
# take mean of models being
# load mean predictions from combos as # round predictions (consensus) # make matrix for overall # make matrix for overall presence # make matrix for overall absence
cbind(model[,1],RoundAll)->temp3 # combine obs with ensemble rowSums(temp3)->temp4 # sum rows (obs + ensemble), if = 2 then true presence, if = 1 then disagree, if = 0 then true absence as.matrix(temp4)->tempsum # place into matrix called tempsum for (k in 1:iters){
# start loop
OC[k,]<-if (RoundAll[k,]==model[k,1]) 1 else 0 # compare combo (ie. ensemble) to obs for each observation in matrix (iters by predictions), if they agree then 1, otherwise 0 presence[k,]<- if (tempsum[k,]==2) 1 else 0 absence[k,]<-if (tempsum[k,]==0) 1 else 0
# make vector of true presences # make vector of true absences
} (sum(OC/iters))->TrueOC # rate of overall classification sd(combo1)->SDModels # SD of ensemble #error <- qt(0.95,df=length(combo1$vals)-1)*sd(combo1$vals)/sqrt(length(combo1$vals)) (sum(presence/iters*2))->Truepr (sum(absence/iters*2))->Trueab Truepr Trueab TrueOC
# sensitivity rate # specificity rate
136
Chapter 6: Using consensus methods to identify (and reduce) sensitivity from methodological choices when measuring functional diversity
Abstract Functional diversity indices have become important tools for measuring variation in species characteristics that are relevant for ecosystem services. Recently, a popular method for measuring functional diversity, FD, was shown to be sensitive to methodological choices in its calculation. The objective of this study was to determine whether consensus methods can be used to identify situations where methodological choices may be an issue when measuring dendrogram-based functional diversity, FD. To calculate FD, a distance measure and a clustering method must be chosen. Using data from natural communities, this study demonstrates that consensus methods were able to determine instances where the choice of distance measure (Euclidean and cosine) and clustering method (UPGMA, complete and single linkage) produced qualitatively different relationships across communities and markedly different dendrogram topologies. In particular this study highlights how consensus methods may aid in the choice of a particular index of functional diversity with the hope that such discussions may improve biodiversity-ecosystem studies.
Keywords: functional diversity, clustering, distance measures, community composition, biodiversity, ecosystem productivity, ecological organization; index.
137
Introduction The association between ecosystem properties and levels of species diversity is well studied (Tilman 1997, Symstad 1998, Loreau et al. 2001). This association was thought to be driven by the tendency for species-rich communities to have wider variation in functional traits (Diaz and Cabido 2001, Heemsbergen et al. 2004, Hooper et al. 2005). The importance of functional traits has led to the development of indices of functional diversity aimed at quantifying functional trait variation; the evaluation of functional diversity indices continues to be an area of active research (Walker 1999, Symstad 2000, Petchey and Gaston 2002, Mason et al. 2003, 2005, Naeem and Wright 2003, Botta-Dukat 2005, Mouillot et al. 2005). Several studies (Cardinale et al. 2000, Petchey and Gaston 2002) have concluded that ecosystem function tend to correlate more strongly with functional diversity indices than with species diversity indices. Results such as these have spurred interest in developing new and improved functional diversity indices that incorporate ecosystem functions (Wright et al. 2006). To calculate most functional diversity indices, a method is required for quantifying interspecific differences in functional traits. In cases where there is only one trait of interest, simple approaches may be appropriate, such as the weighted-trait variation (FDVar; Mason et al. 2003) or the functional evenness (known as functional regularity; Mouilet et al. 2005). However, the flexibility to use more than one trait often is required to understand even simple natural systems and in such cases, the inclusion of trait matrices, distance measures, and sometimes dendrograms, is required. Unfortunately, the use of these multivariate statistical procedures introduces complications that require researchers to make several key decisions for data analysis. Ultimately, these decisions should have minimal effect on patterns of species characteristics as they relate to ecosystem function. However, recent studies (e.g. Poos et al. 2009) have shown that these methodological decisions may be more important than thought here to fore. Not surprisingly there are countless ways to incorporate multivariate functional variation into measures of diversity. A popular index, known as FD (Petchey and Gaston 2002), measures functional diversity as the total branch length of a dendrogram based on functional traits. To produce a dendrogram several decisions need to be made. First, the number and type of traits important to ecosystem function need to be identified. Second, a distance measure needs to be chosen that characterizes the
138
relative differences among species based on their traits. Finally, a clustering algorithm is needed to produce a dendrogram that hierarchically segregates species into functional groups based on their relative distances (Petchey and Gaston 2002). This method has been criticized recently for having the additional subjective step of clustering traits onto a dendrogram (Podani and Schmera 2006; Poos et al. 2009). Although standard methods may provide one way to reduce subjectivity, it is unlikely that a single distance measure or clustering algorithm can be used in all circumstances (Poos et al. 2009). Therefore, methods should be sought which can be used to identify sensitivity. Despite the common use of FD as a measure, claims related to whether FD can be used to derive ecologically robust conclusions have only been quantitatively evaluated recently (Podani and Schmera 2006; Poos et al. 2009). For example, Poos et al. (2009) showed that the probability of two random assemblages showing contrasting levels of functional diversity ranged from 0 to as high as 97.6%. Recently, consensus methods have been suggested as a means for providing a standardized approach for dealing with variation in methodological issues (Mouchet et al. 2008; Poos and Jackson, submitted). However, there are difficulties when employing consensus methods to methodological choices in functional diversity, such as how many initial methods are needed to build the optimal consensus (Poos et al. 2009), and to what degree adding poorly fitted models may undermine a consensus approach (Poos and Jackson, submitted). Therefore the objective of this study is to determine whether consensus methods can be used to understand when sensitivity in measuring FD may be an issue. For this purpose, sensitivity is defined in two ways. First, a sensitive index is defined as one where the same qualitative trends in functional diversity across communities do not persist despite methodological choices (e.g. distance measure and clustering algorithm). This definition is meant to compare the broad ecological consequences of applying FD to ecological communities, with the implicit assumption that a robust index should provide qualitatively repeatable trends. Second, a sensitive index is defined as one where the identified dendrogram topologies persist. This definition requires that the produced patterns of species groupings are maintained. If FD is not robust in this sense, it would suggest that the implication of using FD may be unclear. Explicit recognition of the effects of using a dendrogram, and the decisions needed to get there (e.g. choosing a distance measure and clustering algorithm) need to be better understood, so appropriate guidelines for making these decisions can be formulated.
139
Methods In this study, the same five data sets used in previous studies of FD (Petchey and Gaston 2002, Podani and Schmera 2006; Poos et al. 2009). These datasets represent a variation in the number and type of species (from 13 to 37), and the number and type of functional traits (from 6 to 27). For example, the three vertebrate datasets use characteristics ranging from foraging behavior to the consumption of prey species as their functional traits (Holmes et al. 1979, Jaksic and Medel 1990, Munoz and Ojeda 1997), whereas, the remaining two datasets rely on vegetative characteristics, such as rooting depth and herbivore palatability, of the plants being studied (Golluscio and Sala 1993, Chapin et al. 1996). The measure functional diversity, FD, is based on the total branch length of a dendrogram of functional traits. To obtain this dendrogram, species traits must be assigned a distance (or resemblance) measure and clustering algorithm. Distance measures quantify the association between two entities based on their characteristics (e.g. species based on their functional traits). There are a large number of distance measures from which to choose depending on the data (Jackson et al. 1989; Legendre and Legendre 1998; Podani 1999). Two distance measures were used: Euclidean distance as suggested by Holmes et al. (1979), and cosine distance. Cosine distance was used because it down-weights the potential over-fit created by covarying traits (Legendre and Legendre 1998), a problem often encountered when analyzing functional traits of species (Petchey and Gaston 2006), whereas, Euclidean distance emphasizes larger values, in particular where positive covariance exists between traits. All trait matrices were standardized so that all traits have a mean = 0 and variance =1 (i.e. z-scores; Holmes et al. 1979, Gaston and Petchey 2002). Variability in ecological data is often associated with just a few entities of which clustering into key groupings can provide insight, such as the clustering of species based on functional traits (Legendre and Legendre 1998; Podani 1999). Three clustering algorithms were used in this analysis, unpaired pair group method with arithmetic mean (UPGMA), single linkage (i.e. nearest neighbor) and complete linkage (i.e. maximum or farthest neighbor). These algorithms represent natural endpoints across a methodological continuum of hierarchical clustering algorithms, where single linkage lies on one end, complete linkage on the other, and UPGMA lies somewhere in the middle (Gordon 1999, Podani and Schmera 2006; Poos et al. 2009).
140
Using Consensus Methods to Identify Uncertainty when Measuring FD To determine whether distance measure or clustering algorithm influenced dendrogram topologies of FD, a routine was developed (in MatLAB version 7.1) to randomize the removal of n species from the dataset and recalculate FD for each species combination, clustering algorithm and distance measure. Each level of n species was replicated 1000 times. Functional diversity (FD) was measured at each species richness interval as the total distance of branches in the dendrogram. As FD measures the total branch lengths of a functional dendrogram, which relies on clustering method and distance measure, all dendrograms were rescaled to value between 0-1 using the full species model. The range in FD at the full species level at each different clustering method and distance measure was summarized. The initial dendrogram was compared to each variation in distance measure using consensus trees (Margus and McMorris 1981, Rohlf 1982) in NT-SYS (Rohlf 1997). Dendrograms were compared using the consensus index CI(C) (Rohlf 1982, 1997). Unlike cophenetic correlation (e.g. Mouchet et al. 2008), which compares a dendrogram to the un-modeled raw data, the consensus index compares the similarity of dendrograms based on their cluster membership (Sokal and Rohlf 1962, Shao and Soskal 1986, Legendre and Legendre 1998). The 50% majority rules consensus index was used where a value of one indicates all subgroups share at least 50% membership (i.e. the consensus tree is completely bifurcated indicating similar topology between the original trees) and a value of zero indicates no subgroups are shared (Jackson et al. 1989, Lapointe and Legendre 1990). Although a more strict measure of consensus can be used, the use of a 50% majority rules leads to a more liberal assessment of the similarity between trees than a strict measure would provide.
141
Results The Relationship between FD, Distance Measure & Clustering Algorithm The relationship between FD and community type changed with distance measure and clustering algorithm. In particular, a change in distance measure caused communities to be ranked differently with respect to FD (Figure 6.1). Differences in levels of FD between communities changed with either a change in distance measure or a change in clustering method. For example, FD for predatory vertebrate communities was more similar to Patagonian forb communities using UPGMA and Euclidean distance, but more similar to intertidal fish communities when the clustering algorithm was changed to complete linkage (Figure 6.1). Similarly, the Patagonian forb communities showed similar levels of FD with three different communities, depending on the distance measure and clustering method (Figure 6.1). Rescaling FD did not change these conclusions. Altering the clustering algorithm had a greater impact on the measured variation of FD than altering of distance measure. On average, the overall choice of distance measure and clustering algorithm accounted for a range of 27.4% in the measured amount of functional diversity. At maximum species richness measured functional diversity ranged between a minimum of 21% and a maximum of 61% (mean 34.2%) across clustering methods. Similarly, at maximum species richness the measured amount of functional diversity ranged between a minimum of 12% and a maximum of 41% (mean 20.5%) based on distance measure (Figure 6.2). Qualitatively, single linkage showed the smallest amount of functional diversity, whereas, complete linkage showed the largest.
142
Figure 6.1 - The relationship between distance measure (Euclidean or cosine) and clustering algorithm (SL = single linkage / nearest neighbor, UPGMA = unweighted pair group method with arithmetic mean, CL = complete linkage / farthest neighbor) with FD using five community data sets: A) Arctic vegetation (Chapin et al. 1996); B) Insectivorous birds (Holmes et al. 1979); C) Patagonian forbs (Golluscio and Sala 1993); D) Intertidal fish (Munoz and Ojeda 1997); and, E) Predatory vertebrates (Jaksic and Medel 1990). FD values were re-scaled relative to the Arctic vegetation data, which has the highest FD values. This standardization leads to the appearance of a constant outcome for the Arctic dataset, but this consistency is solely an artifact of using it as the reference point rather than the outcomes not differing depending on the resemblance measure or clustering algorithm.
There was large variation in the measured amount of functional diversity as it relates to the removal of species (Figure 6.2). Not surprisingly, FD was related to species richness, with communities with less species showing smaller ranges of FD. In general, there was a greater similarity between the fitted functional curves for FD and clustering algorithms than the functional curves for FD and similar distance measures. There were two exceptions to this generality. First, the predatory vertebrate dataset showed functional relationships more closely related to similar distance measures (Figure 6.2D), whereas, the curve fitted with UPGMA clustering algorithm and Euclidean distance of the intertidal fishes dataset was more closely related to the curve based on cosine distance and single linkage, than a similar clustering algorithm or distance measure. Overall, conclusions regarding the qualitative relationship of functional diversity among communities were not robust to methodological choices (i.e. FD was not consistent across communities).
143
A)
B)
C)
D)
E) Figure 6.2 - The relationship between distance measure (solid lines = Euclidean distance, dashed lines = cosine distance) and building a dendrogram using a clustering algorithm (1 = complete linkage / farthest neighbor, 2 = unweighted pair group method with arithmetic mean / UPGMA, 3 = single linkage / nearest neighbor) where species are individually removed when calculating FD. Five community data sets are shown: A) Arctic vegetation (Chapin et al. 1996), B) Insectivorous birds (Holmes et al. 1979), C) Patagonian forbs (Golluscio and Sala 1993), D) Intertidal fishes (Munoz and Ojeda 1997), and E) Predatory vertebrates (Jaksic and Medel 1990). Shown inset are 50% majority rule consensus trees demonstrating lack of between group fidelity of species where calculating functional diversity using different distance measures, but the same clustering approach.
144
Identifying Sensitivity in FD Using Consensus Methods The overall low value of the CI(C) index indicates that the decision of choosing a distance measure influences the overall dendrogram to such an extent that there was little resemblance between the dendrogram based on Euclidean distance and the dendrogram based on cosine distance (Table 6.1; Figure 6.2 inset). A comparison of consensus values indicates that FD identified different dendrogram groups (only 46-51% similarity) depending on the distance measure or clustering algorithm (unpublished results) used (Table 6.1). Table 6.1 – Consensus measures of dendrogram group fidelity across distance measures (Euclidean and cosine) for each clustering algorithm: single linkage, unweighted pair group method with arithmetic means (UPGMA), and complete linkage. Group fidelity was determined by majority rules consensus tress using CI(C) consensus index.
Data set Insectivorous birds Intertidal fish Patagonian forbs Predatory vertebrates Arctic vegetation Average
No. species 22 13 24 11 37 ---
Complete Linkage 0.55 0.55 0.45 0.22 0.69 0.49
UPGMA 0.35 0.45 0.50 0.55 0.71 0.51
Single Linkage 0.35 0.55 0.41 0.33 0.66 0.46
The clustering algorithm did not improve the similarity between functional topologies. For example, single linkage, unweighted pair group methods using arithmetic mean, and complete linkage all showed similar rates of consensus tree resemblance, regardless of the size of the tree or the dataset used (Table 6.1). Overall, dendrogram topologies were not robust to the choice of distance measure or clustering algorithm.
Discussion Functional diversity has become an important, but controversial focus of research at the boundary between community and ecosystem ecology (Tilman 2000, Mason et al. 2003, Leps et al. 2006; Poos et al. 2009). This study focused on a popular measure of functional diversity: the total branch length of a functional dendrogram, known as FD (Petchey & Gaston 2002). The results demonstrate that consensus methods were able to identify instances where FD was not robust to the choice of distance measure or clustering algorithm. Specifically, consensus methods indicated that both definitions of robustness were not supported. First, the qualitative
145
relationships did not persist across communities (Figure 6.1). Second, the dendrogram topologies of communities measured using FD differed with the decision of clustering method and distance measure (Figure 6.2). These results were all in agreement with previous quantitative analysis of FD (Poos et al. 2009). The ability of consensus methods to determine qualitative differences of measuring functional diversity across communities have ecological repercussions to many biodiversity-ecosystem functioning studies. It is commonly assumed that diversity is a relative concept; that is, the diversity of a community can only be judged vis-à-vis another community (Magurran 2004). In this case, the conclusions reached by comparing communities based on functional diversity would depend on the distance measure used to create the dendrogram and, therefore, relationships among communities would be altered (Figure 6.1). For example, the FD for the Patagonian forbs community was identical to predatory vertebrate community using single linkage and Euclidean distance. However, if the distance measure was altered to cosine distance, the Patagonian forbs community more closely resembled the insectivorous bird communities, whereas, a change in clustering algorithm would show the Patagonian forbs community more closely resembled the intertidal fish community. In this study, a difference in clustering algorithm or distance measure altered community relationships altogether, and these relationships persisted regardless of the scaling used. These differences have a high likelihood of altering ecological interpretations of functional diversity across communities that, in turn, will potentially confound ecosystem-based studies. Different methods of calculating FD can lead to different dendrograms, and consequently different measures of functional diversity. For example, data analyzed using different distance measures and clustering algorithms varied by a range of 27% in the measured amount of functional diversity at maximum species richness (Figure 6.2). Therefore, studies which use dendrogram-based methods of measuring functional diversity may drastically under or over estimate the amount of functional diversity, assuming some “true” value can be determined. Therefore one may assume two communities are more similar using one method when they may not be when based on another measure. Finally, decisions are required regarding how the data are treated - for example, depending on how tied values in a similarity matrix are treated, a number of different dendrograms can be produced (Jackson et al. 1989). As the clustering technique will produce dendrograms (with the underlying goal of determining group structure)
146
whether or not true groups exist (Jackson et al. 1989), the cumulative effect of these decisions to the relevance of the identified groups may be unknowingly large. For example, this study suggests that the decision of adding distance measure and clustering algorithm to produce a dendrogram will alter the range in the measure of FD (e.g. 34.2% for clustering algorithm and 20.5% for distance measure) and the ecological conclusions reached (Fig. 6.1). Identifying areas where sensitivity may be an issue can aid in developing more robust or representative indices of functional diversity. Recently, studies have attempted to use consensus methods, which take the average of several different methodological approaches, to reconcile differences in FD from methodological issues (Mouchet et al. 2008). However, previous applications of consensus methods may not be sufficient in determining where sensitivity may be an issue. In this paper, only 46-51% of the topologies were in agreement (Table 6.1) and there was large variation in the measurement of FD across methodological choices and species (Figure 6.2). Although consensus methods were useful in identifying issues of sensitivity, using consensus methods as a standardized approach may not be appropriate because averaging three poorly associated models together may actually produce a worse outcome (Poos and Jackson submitted). In addition, previous applications of consensus methods provide several more decisions – such as how many initial models are needed to build the optimal consensus, and which consensus method is most appropriate (Poos et al. submitted) – which may further impact FD (Poos et al. 2009). There is considerable debate regarding the most appropriate measure of functional diversity and the qualities that metric should possess (Loreau et al. 2001, Mason et al. 2003, Ricotta 2005, Leps et al. 2006). Clearly, regardless of the index used, any index of functional diversity must be robust (e.g. qualitatively similar across methods and dendrograms) to decisions inherent in its calculation or one must decide upon a common statistical methodology in order to permit comparisons amongst studies. Quantitative comparisons of how functional diversity indices differ are rare (e.g. Petchey et al. 2004; Walker et al. 2008), and evaluations of other functional diversity indices are needed. In calculating FD, the decisions inherent in its calculation represent two additional difficulties aside from previous criticisms of which species, what kind of diversity, and which ecosystem function (Bengtsson 1998, Symstad et al. 1998, Cardinale 2000, Jax 2005) are to be included and therefore should also define the choice of similarity measure and clustering algorithm. Further criticisms, such as how many functional traits (Walker et al. 1999, Podani and Schmera 2006), what qualifies as a functional
147
group (Petchey and Gaston 2006), and what type of consensus approach, also apply. Explicit recognition and justification of each of these decisions is warranted for improving functional diversity research.
Acknowledgements Funding NSERC & OGS Scholarships to M.S.P and S.C.W., NSERC Discovery Grant to D.A.J., Ontario Ministry of Natural Resources, and University of Toronto.
References Bengtsson, J. 1998. Which species? What kind of diversity? Which ecosystem function? Some problems in studies of relations between biodiversity and ecosystem function. Applied Soil Ecology 10: 191-199. Botta-Dukat, Z. 2005. Rao's quadratic entropy as a measure of functional diversity based on multiple traits. Journal of Vegetation Science 16: 533-540. Cardinale, B. J., K. Nelson, and M.A. Palmer. 2000. Linking species diversity to the functioning of ecosystems: On the importance of environmental context. Oikos 91: 175-183. Chapin, F. S. I., M.S. Bret-Harte, S.E. Hobbie, and Z. Hailan. 1996. Plant functional types as predictors of transient responses of arctic vegetation to global change. Journal of Vegatation Science 7: 347-358. Diaz, S., and M. Cabido. 2001. Vive la difference: plant functional diversity matters to ecosystem processes. Trends in Ecology & Evolution 16: 646-655. Golluscio, R. A., and O.E. Sala. 1993. Plant functional types and ecological strategies in Patagonian forbs. Journal of Vegatation Science 4: 839-846. Gordon, A. D. 1999. Classification, 2nd ed. Chapman and Hall, London, United Kingdom. Heemsbergen, D. A., M.P. Berg, M. Loreau, J.R. van Hal, J.H. Faber, and H.A. Verhoef. 2004. Biodiversity effects on soil processes explained by interspecific functional dissimilarity. Science 306: 1019.
148
Holmes, R. T., R.E.J. Bonney, and S.W. Pacala. 1979. Guild structure of the Hubbard Brook bird community: a multivariate approach. Ecology 60: 512-520. Hooper, D. U., F.S. Chapin, J.J. Ewel, A. Hector, P. Inchausti, S. Lavorel, J.H. Lawton, D.M. Lodge, M. Loreau, S. Naeem, B. Schmid, H. Setala, A.J. Symstad, J. Vandermeer, and D.A. Wardle. 2005. Effects of biodiversity on ecosystem functioning: A consensus of current knowledge. Ecological Monographs 75: 3-35. Jackson, D. A., K.M. Somers, and H.H. Harvey. 1989. Similarity coefficients: measures of cooccurrence and association or simply measures of occurrence? American Naturalist 133: 436453. Jaksic, F. M., and R.G. Medel. 1990. Objective recognition of guilds: testing for statistically significant species clusters. Oecologia 82: 87-92. Jax, K. 2005. Function and "functioning" in ecology: What does it mean? Oikos 111: 641-648. Lapointe, F.J., and P. Legendre. 1990. A statistical framework to test the consensus of two nested classifications. Systematic Zoology 39: 1-13. Legendre, P., and Legendre, L. 1998. Numerical ecology, 2nd ed. Elsevier B.V. Leps, J., F., de Bello, S. Lavorel, and S. Berman. 2006. Quantifying and interpreting functional diversity of natural communities: practical considerations matter. Preslia 78: 481-501. Loreau, M., S. Naeem, P. Inchausti, J. Bengtsson, J.P. Grime, A. Hector, D.U. Hooper, M.A. Huston, D. Raffaelli, B. Schmid, D. Tilman, and D.A. Wardle. 2001. Biodiversity and ecosystem functioning: Current knowledge and future challenges. Science 294: 804-808. Margurran, A. E. 2004. Measuring Biological Diversity. Blackwell Publishing, Oxford, United Kingdom. Margus, T., and F.R. McMorris. 1981. Consensus n-trees. Bulletin of Mathematical Biology 43: 239-244.
149
Mason, N. W. H., K. MacGillivray, J.B. Steel, and J.B. Wilson. 2003. An index of functional diversity. Journal of Vegetation Science 14: 571-578. Mason, N. W. H., D. Mouillot, W.G. Lee, and J.B. Wilson. 2005. Functional richness, functional evenness and functional divergence: the primary components of functional diversity. Oikos 111: 112-118. Mouillot, D., W.H.N. Mason, O. Dumay, and J.B. Wilson. 2005. Functional regularity: a neglected aspect of functional diversity. Oecologia 142: 353-359. Munoz, A. A., and F.P. Ojeda. 1997. Feeding guild structure of a rocky intertidal fish assemblage in central Chile. Environmental Biology of Fishes 49: 471-479. Naeem, S., and J.P. Wright. 2003. Disentangling biodiversity effects on ecosystem functioning: deriving solutions to a seemingly insurmountable problem. Ecology Letters 6: 567-579. Petchey, O. L., and K.J. Gaston. 2002. Functional diversity (FD), species richness and community composition. Ecology Letters 5: 402-411. Petchey, O. L., and K.J. Gaston. 2006. Functional diversity: back to basics and looking forward. Ecology Letters 9: 741-758. Petchey, O.L., A. Hector and K.J. Gaston. 2004. How do different measures of functional diversity perform. Ecology 85: 847-857. Podani, J. 1999. Extending Gower's general coefficient of similarity to ordinal characters. Taxon 48: 331-340. Podani, J., and D. Schmera. 2006. On dendrogram-based measures of functional diversity. Oikos 115: 179-185. Poos, M. S., and D. A. Jackson. Conservation by consensus: reducing uncertainties in modeling the distribution of an endangered species using habitat-based ensemble models. Ecological Applications (submitted).
150
Poos, M. S., S. W. Walker, and D. A. Jackson. 2009. Functional-diversity indices can be driven by methodological choices and species richness. Ecology 90:341-347. Ricotta, C. 2005. A note on functional diversity measures. Basic and Applied Ecology 6: 479486. Rohlf, F. J. 1982. Consensus indices for comparing classifications. - Mathematical Biosciences 59: 131-144. Rohlf, F. J. 1997. NTSYSpc Numerical Taxonomy and Multivariate Analysis System. Vol. 2.0. Exeter Software, Setauket, New York. Shao, K. and R.R. Sokal. 1986. Significance tests of consensus indices. Systematic Zoology 35: 582-590. Simpson, E. H. 1949. Measurement of diversity. Nature 163: 688. Sokal, R. R., and F.J. Rohlf. 1962. The comparisons of dendrograms by objective methods. Taxon 11:33-40. Symstad, A. J., D. Tilman, J. Willson, and J.M.H. Knops. 1998. Species loss and ecosystem functioning: Effects of species identity and community composition. Oikos 81: 389-397. Symstad, A. J. 2000. A test of the effects of functional group richness and composition on grassland invasibility. Ecology 81: 99-109. Tilman, D., J. Knops, D. Wedin, P. Reich, M. Ritchie, and E. Siemann. 1997. The Influence of functional diversity and composition on ecosystem processes. Science 277: 1300-1303. Tilman, D. 2000. Causes, consequences and ethics of biodiversity. Nature 405: 208-211. Walker, B.H., A. Kinzig, and J. Langridge. 1999. Plant attribute diversity, resilience, and ecosystem function: the nature and significance of dominant and minor species. Ecosystems 2: 95-113. Walker, S.C, M.S. Poos, and D.A. Jackson. 2008. Functional rarefaction: estimating functional diversity from field data. Oikos 117: 286-296.
151
Wright, J. P., S. Naeem, A. Hector, C. Lehman, P.B. Reich, B. Schmid, and D. Tilman. 2006. Conventional functional classification schemes underestimate the relationship with ecosystem functioning. Ecology Letters 9: 111-120.
152
Appendices Appendix 6.1 – MatLAB Code for calculating total branch lengths of dendrograms from various species combinations function branches = branch_lengths(Z) linkages = size(Z,1); species = linkages+1; branches = zeros(linkages, 2); for i = 1:linkages if (Z(i,1) < (1+species)) & (Z(i,2) < (1+species)) branches(i,1) = Z(i,3); branches(i,2) = Z(i,3); elseif (Z(i,1) > species) & (Z(i,2) < (1+species)) branches(i,1) = Z(i,3) - Z((Z(i,1)-species),3); branches(i,2) = Z(i,3); elseif (Z(i,1) < (1+species)) & (Z(i,2) > species) branches(i,1) = Z(i,3); branches(i,2) = Z(i,3) - Z((Z(i,2)-species),3); else branches(i,1) = Z(i,3) - Z((Z(i,1)-species),3); branches(i,2) = Z(i,3) - Z((Z(i,2)-species),3); end end
153
Chapter 7: General Conclusions This thesis provides novel advancement in determining the influence of methodological choices in conservation-based models and on how consensus methods may reduce uncertainty in such models. In particular, this thesis demonstrates that regardless of the scale (local as in Chapter 5 or regional in Chapter 6), of the species, or study system in question, methodological choices have the ability to dramatically impact resultant analyses and ecological inquiry. For example, in Chapter 2, methodological choices had the ability to provide levels of sensitivity over 97% in the measure of functional diversity. Similarly, in Chapter 4, methodological choices changed estimates of population viability by several orders of magnitude. In Chapter 3, the addition or removal of rare species did not impact multivariate analyses as much as choice of distance measure or multivariate method, which had the ability to drastically alter bioassessments. As these findings are novel, this chapter will highlight some conclusions regarding the importance of methodological choices and attempt to provide recommendations for minimizing methodological impacts.
A) Conservation-based Models in General The importance of understanding the impact of methodological choices is not only timely, but essential. With complex statistical software readily available, ecologists now have several dozen approaches to choose from when developing conservation-based models. Perhaps not surprisingly, this thesis demonstrates, as have previous authors, the importance of model type in impacting results (e.g. Jackson et al. 1989; Jackson 1993; Guissan and Zimmerman 2000; Thuiller 2005; Elith et al. 2006; Dormann et al. 2008; Marrimon et al. 2009). For example, in Chapter 3, multivariate technique impacted analyses more than choice in distance measure or the removal of rare species. This is in agreement with others who have demonstrated that rare species may provide useful information (e.g. Cao et al 1998; 2001) and that selection of multivariate method may influence interpretation strongly (e.g. Jackson et al. 1989; Jackson 1993). However, this thesis also demonstrates that despite laborious efforts in developing modeling comparisons (Olden and Jackson 2002; Thuiller et al. 2005; Elith et al. 2006; Sharma and Jackson 2008; Marrimon et al. 2009; this thesis), there has been no evidence of the preeminence of any singular methodological approach. Issues related to methodological impacts
154
will vary depending on the dataset and objectives in question, such as which distance measure, clustering algorithm, ordination (or other). Therefore, methods that can be used to reduce the influence of methodological choices, or highlight when methodological choices may be an issue - such as the use of consensus models - will help recognize and reduce the impact of methodological choices.
B) Functional Diversity Incorporating functional traits can be another way of improving modeling approaches for species with conservation concern. Species’ functional characteristics strongly influence ecosystem properties (Loreau et al. 2001; Hooper et al. 2005), and the understanding of the relationships between functional diversity and community structure has been important for identifying mechanisms of biodiversity effects. To date, the focus on functional traits has been on dominant species (Grime 1998); although rare (presumably including imperiled) species can have a large influence on ecosystem processes (Power et al. 1996). Methodological choices in the measure of functional diversity remain a controversial topic (Mason et al. 2003, Ricotta 2005, Podani and Schmera 2006). For example, this thesis (Chapters 2 and 6) demonstrates that the measure of functional diversity is greatly complicated by methodological choices. Therefore, discussion is warranted on the future use of FD as a metric of functional diversity and some of the qualities this metric poses. Ultimately, the choice of what to do regarding the impact of methodological choices on measuring functional diversity may rely on the advantages and disadvantages of the features of each approach. Many of these features have been discussed previously (Mason et al. 2003, Hooper et al. 2005, Mouillot et al. 2005, Ricotta 2005, Leps et al. 2006, Petchey and Gaston 2006; Schmera et al. 2009a; 2009b) and here, properties of metrics of functional diversity are discussed plus are new insights into areas worthy of future research. One issue with metrics of functional diversity is whether or not they increase monotonically with species richness (known as ‘set monotonicity’; Ricotta 2005; Schmera et al. 2009a). For example, unlike Rao’s quadratic entropy (Rao 1982, Botta-Dukat 2005), FD has the intuitive property of a monotonic relationship with species richness (Petchey and Gaston 2006). FD cannot decrease when a species is added to a community, and when a species is removed FD
155
cannot increase (Petchey and Gaston 2002). However, in some cases FD violates this feature (Podani and Schmera 2006) and Walker et al. (1999, 2008) have noted that functional attribute diversity (FAD) previously has been misidentified as lacking this feature. Recent advancements in measuring FAD (Schmera et al. 2009a) have helped to clarify the issue of set monotonicity in FAD; however, future research into whether set monotonicity is a requirement for robust estimation of functional diversity remains to be determined. For example, Rao’s quadratic entropy is gaining currency as a flexible index of functional diversity as it is an intuitive extension of Simpson’s index of diversity (Simpson 1949, Botta-Dukat 2005, Leps et al. 2006), yet does not have set monotonicity (Pavoine and Bonsall 2009). Given that a change in distance measure caused a related change in the qualitative relationships of functional diversity across communities (in Chapter 6), it bears asking to what degree a distance measure is required. Although the choice of distance measure involves subjectivity that may influence the analysis, the inclusion of a distance measure allows for a continuous segregation of multiple species based on multiple functional traits (Legendre and Legendre 1998). Therefore, the question of whether or not to include a distance measure depends to what extent ecological diversity is based on the trait dissimilarity among species in a community (Tilman 1997, Petchey and Gaston 2002) and to what degree a distance measure can distinguish those traits and/or species. Few studies relate trait dissimilarity among species in a community (Petchey and Gaston 2002, Garnier et al. 2004, Heemsbergen et al. 2004, Leps et al. 2006), and unfortunately, even fewer studies determine the degree from which distance measure can distinguish those trends. One mechanism to lower the subjectivity of the calculation in the choice of distance measure is to provide standards based on the purpose of the study. Where traits have a mixture of data types the use of Gower similarity may more applicable (Podani and Schmera 2006), although quantitative comparisons with other distance measures are still needed. Furthermore, unlike Euclidean distance, Gower similarity can be used with missing values and has the advantage of not being influenced by the unit of measure (Gower 1971, Legendre and Legendre 1998, Podani 1999). Future research should focus on quantitative comparisons of pairwise distance measures and ensuring they retain strong linkages to raw data on species traits. The treatment of functionally redundant species in functional diversity indices remains controversial. Central to this debate is whether an index of functional diversity should change if a species is added or lost that is identical to one already present; a feature known as functional
156
redundancy (Walker et al. 2008; Schmera et al. 2009b). Petchey and Gaston (2006) contend that a good index of functional diversity should not increase if a redundant species is added. However, this feature may be a product of the number of traits used to segregate species, and the ecosystem function in question, and not necessarily because the species provide no additional function (Rosenfield 2002). More recently, improvements to the calculation of functional diversity indices have allowed for comparison of functional diversity across communities without sampling bias (Walker et al. 2008). This feature may provide additional benefit to functional diversity indices as it deemphasizes the current trend in identifying functionally redundant species and emphasizes the need to compare functional traits across broader assemblages (Rosenfeld 2002). There is a growing consensus that functional diversity is likely to be the component of biodiversity most relevant to ecosystem function (Wright et al 2006). Recently, biodiversity theory and management perspectives have converged, where each has embraced the need for incorporating species-specific biology into models (Srivastava and Velland 2005). For example, large-bodied species that occupy high trophic positions in food webs and occur at low abundance are thought to be particularly vulnerable to extinction (Lawton and May 1995). The application of functional traits, such as these, to improve models is becoming more widespread (Cardillo and Bromham 2001; Olden et al. 2006), and their ability to improve models for species with conservation concern needs to e assessed. The consideration of both diversity and evenness of species traits may be an important consideration as evenness measures not only the range of functional variation, but how much of the functional variation is filled within that range (Mason et al. 2003 Mouillot et al. 2005). Newer approaches that disassemble functional evenness from functional richness (e.g. Walker et al. 2008) and allow the researcher to distinguish one from the other, may allow for better application of functional diversity indices to research on species with conservation concerns.
157
Recommendations •
Methodological choices should be closely scrutinized when modeling species with conservation concern. As this thesis (Section I) demonstrates, regardless of scale, study system or ecological question, methodological choices had the ability greatly to impact results.
•
The use of singular statistical methods should be avoided when dealing with species where data are deficient (as often is the case of species with conservation concern). This dependency was true regardless of the scale - community level as in Chapters 2, 3, and 6, population level as in Chapter 5, or the metapopulation level in Chapter 4. The need for comparative approaches using multiple statistical methods has been highlighted as one means to reduce problems related to conservation-based models (Guisan and Zimmerman 2000; Olden and Jackson 2002).
•
In the case of modeling species with conservation concerns, singular statistical methods should not be interpreted in isolation. As Chapter 6 demonstrated, several methods provided contrasting explanatory relationships. These relationships, if taken in isolation, may bias future conservation efforts in areas where species are not showing strong habitat relationships.
•
This thesis demonstrates that consensus methods may provide reduced uncertainty in modeling species with conservation concerns. The advantages of consensus modeling over singular approaches are numerous : o It is an intuitive extension of modeling using singular approaches; o Advancements in computing power have accelerated the use of new statistical models, which can be added to this approach (e.g. random forest, boosted regression); o It reduces biases based on choosing singular statistical approaches, and the resultant limitations of fitting data to that approach (especially if assumptions are not met);
158
o It produces a prioritized model output that can be used to address conservation priorities, identify areas of high conservation value … etc., even if there are differences across the models and data is noisy; o Model thresholds can be adjusted readily to prioritize the correct classification of species presence or absence (e.g. a balance approach was used here in); and, o Misclassifications can be identified and plotted spatially to identify areas where error rates were high and data can be collected to refine the models.
References Araújo, M. B., and A. Guisán. 2006. Five (or so) challenges for species distribution modelling. Journal of Biogeography 33:1677-1688. Araújo, M. B., and M. New. 2007. Ensemble forecasting of species distributions. Trends in Ecology and Evolution 22:42-47. Araújo, M. B., R. J. Whittaker, R. J. Ladle, and M. Erhard. 2005. Reducing uncertainty in projections of extinction risk from climate change. Global Ecology and Biogeography 14:529-538. Botta-Dukat, Z. 2005. Rao's quadratic entropy as a measure of functional diversity based on multiple traits. Journal of Vegetation Science 16: 533-540. Cao Y., D.D. Williams, and N.E. Williams. 1998. How important are rare species in community ecology and bioassessment. Limnology and Oceanography 43: 1403–1409. Cao Y., D.P. Larsen, and R.S. Thorne. 2001. Rare species in multivariate analysis for bioassessment: some consideration. Journal of the North American Benthological Society 20: 144–153. Cardillo, M., and L. Bromham. 2001. Body size and risk of extinction in Australian mammals. Conservation Biology 15: 1435-1440. Cunningham, R.B., and D.B. Lindenmayer. 2005. Modeling count data of rare species: some statistical issues. Ecology 86: 1135-1142.
159
Dormann, C. F., O. Purschke, J. R. G. Márquez, S. Lautenbach, and B. Schröder. 2008. Components of uncertainty in species distribution analysis: a case study of the great grey shrike. Ecology 89:3371-3386. Elith, J., C. H. Graham, R. P. Anderson, M. Dudík, S. Ferrier, A. Guisan, R. J. Hijmans, F. Huettmann, J. R. Leathwick, A. Lehmann, J. Li, L. G. Lohmann, B. A. Loiselle, G. Manion, C. Moritz, M. Nakamura, Y. Nakazawa, J. M. Overton, A. T. Peterson, S. J. Phillips, K. Richardson, R. Scachetti-Pereira, R. E. Schapire, J. Soberón, S. Williams, M. S. Wisz, and N. E. Zimmermann. 2006. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29:129-151. Ellison, A.M., and A.A. Agrawal. 2005. The statistics of rarity. Ecology 86: 1079-1080. Grime, J.P. 1998. Benefits of plant diversity to ecosystems: immediate, filter and founder effects. Journal of Ecology 86: 902-910. Garnier, E., J. Cortez, G. Billes, M. L. Navas, C. Roumet, M. Debussche, G. Laurent, A. Blanchard, D. Aubry, A. Bellmann, C. Neill, and J. P. Toussaint. 2004. Plant functional markers capture ecosystem properties during secondary succession. Ecology 85:26302637. Gower, J.C. 1971. A general coefficient of similarity and some of its properties. Biometrics 27: 857_874. Guisan, A., and N. E. Zimmermann. 2000 Predictive habitat distribution models in ecology. Ecological Modelling 135:147-186. Green, R.H. and R.C. Young. 1993. Sampling to detect rare species. Ecological Application 3: 351-366. Heemsbergen, D. A., M.P. Berg, M. Loreau, J.R. van Hal, J.H. Faber, and H.A. Verhoef. 2004. Biodiversity effects on soil processes explained by interspecific functional dissimilarity. Science 306: 1019. Hooper, D. U., F. S. Chapin, J. J. Ewel, A. Hector, P. Inchausti, S. Lavorel, J. H. Lawton, D. M. Lodge, M. Loreau, S. Naeem, B. Schmid, H. Setala, A. J. Symstad, J. Vandermeer, and
160
D. A. Wardle. 2005. Effects of biodiversity on ecosystem functioning: A consensus of current knowledge. Ecological Monographs 75: 3-35. Jackson, D.A. 1993. Multivariate analysis of benthic invertebrate communities: the implication of choosing particular data standardizations, measures of association, and ordination methods. Canadian Journal of Fisheries and Aquatic Sciences 50: 2641-2651. Jackson D. A., K.M. Somers, and H.H. Harvey. 1989. Similarity coefficients: measures of cooccurrence and association or simply measures of occurrence? American Naturalist 133: 436453. Laplace, P. S. 1820. Théorie analytique des probabilités. Courcier, Paris. Lawton, J.H. and R.M. May (eds.). 1995. Extinction rates. Oxford University Press, London, United Kingdom. Legendre, P. and E.D. Gallagher. 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129: 271-280 Legendre, P., and L. Legendre. 1998. Numerical Ecology, 2nd ed. - Elsevier B.V. Leps, J., F. de Bello, S. Lavorel, and S. Berman. 2006. Quantifying and interpreting functional diversity of natural communities: practical considerations matter. Preslia 78: 481-501. Loreau, M., S. Naeem, P. Inchausti, J. Bengtsson, J. P. Grime, A. Hector, D. U. Hooper, M. A. Huston, D. Raffaelli, B. Schmid, D. Tilman, and D. A. Wardle. 2001. Biodiversity and ecosystem functioning: current knowledge and future challenges Science 294: 804-808. Marmion, M., M. Parviainen, M. Luoto, R. K. Heikkinen, and W. Thuiller. 2009. Evaluation of consensus methods in predictive species distribution modelling. Diversity and Distributions 15: 59-69. Mason, N. W. H., K. MacGillivray, J.B. Steel, and J.B. Wilson. 2003. An index of functional diversity. Journal of Vegetation Science 14: 571-578.
161
Mouillot, D., W.H.N. Mason, O. Dumay, and J.B. Wilson. 2005. Functional regularity: A neglected aspect of functional diversity. Oecologia 142: 353-359. Olden, J. D., and D. A. Jackson. 2002. A comparison of statistical approaches for modelling fish species distributions. Freshwater Biology 47:1976-1995. Olden, J.D., N.L. Poff, and B.P. Bledsoe. 2006. Incorporating ecological knowledge into ecoinformatics: An example of modeling hierarchically structured aquatic communities with neural networks. Ecological Informatics 1: 33-42 Pavoine, S., and M. B. Bonsall. 2009. Biological diversity: distinct distributions can lead to the maximization of Rao’s quadratic entropy Theoretical Population Biology 75:153-163. Petchey, O. L. and K.J. Gaston. 2002. Functional diversity (FD), species richness and community composition. Ecology Letters 5: 402-411. Petchey, O. L. and K.J. Gaston. 2006. Functional diversity: back to basics and looking forward. Ecology Letters 9: 741-758. Podani, J. 1999. Extending Gower's general coefficient of similarity to ordinal characters. Taxon 48:331-340. Podani, J. and D. Schmera. 2006. On dendrogram-based measures of functional diversity. Oikos 115: 179-185. Power, M.E., D. Tilman, J.A. Estes, B.A. Menge, W.J. Bond, L.S. Mills, G. Daily, J.C. Castilla, J. Lubchenco, and R.T. Paine. 1996. Challenges in the quest for keystones. BioScience 46: 609-620. Rao, C.R. 1982. Diversity and dissimilarity coefficients—a unified approach. Theoretical Population Biology 21: 24-43. Ricotta, C. 2005. A note on functional diversity measures. Basic and Applied Ecology 6: 479486.
162
Rodríguez, J. P., L. Brotons, J. Bustamante, and J. Seoane. 2007. The application of predictive modelling of species distribution to biodiversity conservation. Diversity and Distributions 13:243–251. Rosenfeld, J. S. 2002. Functional redundancy in ecology and conservation. Oikos 98:156-162. Schmera, D., T. Eros, and J. Podani. 2009a. A measure for assessing functional diversity in ecological communities. Aquatic Ecology 43:157-167. Schmera, D., J. Podani, and T. Eros. 2009b. Measuring the contribution of community members to functional diversity. Oikos 118:961-971 Sharma, S. and D.A. Jackson. 2008. Predicting smallmouth bass incidence across North America: Comparison of statistical approaches. Canadian Journal of Fisheries and Aquatic Sciences 65: 471-481. Simpson, E. H. 1949. Measurement of diversity. Nature 163:688. Srivastava, D. S., and M. Vellend. 2005. Biodiversity-ecosystem function research: Is it relevant to conservation? Annual Review of Ecology, Evolution, and Systematics 36: 267-294. Thuiller, W. 2004. Patterns and uncertainties of species’ range shifts under climate change. Global Change Biology 10: 2020-2027. Thuiller, W., S. Lavorel, M. B. Araújo, M. T. Sykes, and I. C. Prentice. 2005. Climate change threats to plant diversity in Europe. Proceedings of the National Academy of Sciences 102:8245-8250. Tilman, D., J. Knops, D. Wedin, P. Reich, M. Ritchie, and E. Siemann. 1997. The Influence of Functional Diversity and Composition on Ecosystem Processes. Science 277: 1300-1303. Walker, B.H., A. Kinzig, and J. Langridge. 1999. Plant attribute diversity, resilience, and ecosystem function: the nature and significance of dominant and minor species. Ecosystems 2: 95-113.
163
Walker, S. W., M. S. Poos, and D. A. Jackson. 2008. Functional rarefaction: estimating functional diversity from field data. Oikos 117:286-296. Wright, J. P., S. Naeem, A. Hector, C. Lehman, P. B. Reich, B. Schmid, and D. Tilman. 2006. Conventional functional classification schemes underestimate the relationship with ecosystem functioning. Ecology Letters 9:111-120.
