MODELING THE IMPACTS OF BARRIER REMOVAL ON GREAT LAKES SEA
LAMPREY
By
Alexander James Jensen

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Fisheries and Wildlife - Master of Science
2017

ABSTRACT
MODELING THE IMPACTS OF BARRIER REMOVAL ON GREAT LAKES SEA
LAMPREY
By
Alexander James Jensen
Barriers in the Great Lakes represent an effective form of control for the invasive sea
lamprey (Petromyzon marinus) by blocking large extents of river habitat and subsequently
eliminating the need for the lampricide treatments in these upstream areas. With increasing
pressure for barrier removals, the availability of suitable sea lamprey habitat above these barriers
and the expected population response to dam removals represent key uncertainties in decisionmaking. The development and evaluation of models to predict larval habitat quantities using
readily-available, reach-scale landscape predictors improved our understanding of common
influences on stream habitat, but failed to reliably predict habitat proportions upstream of
barriers in the Lake Michigan drainage basin. Subsequent simulation-based modeling of the
Lake Michigan sea lamprey population revealed a disproportionate, exponential response to
increasing habitat availability, driven in part by decreasing overall lampricide treatment
frequencies under a fixed control budget. The same modeling approach was used to generate sea
lamprey population predictions associated with projected removal of Grand River’s Sixth Street
Dam under a suite of alternative management actions and biological assumptions. Based on all
simulation results, barrier removals appear to necessitate a substantial increase in annual
lampricide control costs to prevent disproportionate increases in sea lamprey abundance across
the Lake Michigan basin.

ACKNOWLEDGMENTS

I would like to first acknowledge my advisor Dr. Michael Jones for his constant support
and willingness to enable my tendency to pursue often unrelated lines of research. Those
opportunities may have added some time to my stay at Michigan State University, but they were
absolutely worthwhile. I’d also like to thank my thesis committee members, Dr. Travis Brenden
and Dr. Michael Wagner, for their guidance and constructive comments throughout this research
project. GIS support and landscape analytic advice from Dr. Dana Infante and Kyle Herreman
set the stage for all habitat modeling efforts, and ongoing support from the Quantitative Fisheries
Center was crucial during the development and testing of the modified sea lamprey model,
notably from Dr. Norine Dobiesz. I would also like to acknowledge the involvement of
numerous researchers and managers from the Great Lakes Fishery Commission, many of them
from the Larval Assessment and Barrier Task Forces, in informing and improving my research
efforts. This work was facilitated by contributions from multiple Fisheries and Wildlife
undergraduate students. I would also be rather remiss in failing to acknowledge the support of
my colleagues and fellow graduate students within the Department of Fisheries and Wildlife
throughout my degree. Finally, the Great Lakes Fishery Commission and the Dr. Howard A.
Tanner Fisheries Excellence Fellowship provided funding for this research.

iii

TABLE OF CONTENTS

LIST OF TABLES ......................................................................................................................... vi
LIST OF FIGURES ...................................................................................................................... vii
THESIS INTRODUCTION .............................................................................................................1
CHAPTER 1: COMPARISON OF MODEL PERFORMANCE IN PREDICTING REACHSCALE LARVAL SEA LAMPREY HABITAT ABOVE LAKE MICHIGAN BARRIERS ........8
INTRODUCTION .....................................................................................................................9
METHODS ..............................................................................................................................14
Historical Data ..................................................................................................................14
Landscape Variable Reduction Techniques .......................................................................17
Pre-Modeling Data Transformations ................................................................................17
Fitting the Original Type I and Type II Response Data ....................................................18
Hierarchical Logit Regression Model .........................................................................18
Hierarchical Zero/One Inflated Beta Regression Model .............................................21
Boosted Generalized Additive Model...........................................................................23
Fitting the Derived Response Variables ............................................................................24
Assessing Variable Influence and Model Fit .....................................................................24
Evaluating Model Predictive Ability on Out-of-Sample Data ...........................................26
RESULTS ................................................................................................................................30
Data Transformations ........................................................................................................30
Variable Reduction ............................................................................................................30
Final Dataset Size ..............................................................................................................34
Convergence Diagnostics ..................................................................................................34
Variable Coefficients .........................................................................................................35
Model Fit ............................................................................................................................40
Sampling Success and Model Predictive Ability ................................................................45
DISCUSSION ..........................................................................................................................49
APPENDICES .........................................................................................................................57
APPENDIX A: WinBUGS Script for the Type I Habitat Logit Regression .....................58
APPENDIX B: Partial Effects Plots for Derived Response Variables ..............................60
CHAPTER 2: SIMULATING THE SEA LAMPREY RESPONSE TO DAM REMOVALS IN
THE LAKE MICHIGAN DRAINAGE BASIN ............................................................................62
INTRODUCTION ...................................................................................................................63
METHODS ..............................................................................................................................68
Modifying the SLaMSE Model ...........................................................................................68
Assessing Population Responses to Systematically Increasing Habitat Availability ........69
Modeling Barrier Removal on the Grand River ................................................................72
RESULTS ................................................................................................................................80
Systematic Habitat Increases .............................................................................................80

iv

Case Study: Barrier Removal on the Grand River ............................................................87
DISCUSSION ..........................................................................................................................91
MANAGEMENT RECOMMENDATIONS ...............................................................................102
LITERATURE CITED ................................................................................................................105

v

LIST OF TABLES

Table 1.1: Descriptive key of model parameters describing the Bayesian hierarchical logit
regression and beta regression models. .........................................................................................20
Table 1.2: Summary of aggregated class-based predictor variables considered for model
inclusion. .......................................................................................................................................32
Table 1.3: Final set of predictors considered for modeling. Bolded terms indicate those variables
not considered collinear according to the GVIF analysis, and therefore included in the final
modeling efforts. LB, NB, LC, NC, and reach refer to local buffer, network buffer, local
catchment, network catchment, and reach scales, respectively. Further details on the resolution
and year of creation for these predictor datasets can be obtained in Wang et al. 2011. ...............33
Table 1.4: Variable selection frequency in the boosted generalized additive model fit to Type I
habitat, Type II habitat, total larval habitat, and the Type II to Type I habitat ratio. The bolded
values indicate the three most frequently selected predictor variables for each response. ...........37
Table 1.5: Summary of model fit metrics across the different modeling approaches and habitat
types. The root mean squared deviation (RMSD) is the observed average deviation of the
predicted values from the observed. Intercept (a), slope (b), and R2 were calculated from a linear
regression of observed against predicted values. (D) indicates the model was fit using derived
response variables, and * indicates the intercept/slope parameter fell within the expected range
based on equivalence testing. ........................................................................................................41
Table 1.6: Model predictive fit for Type I and Type II habitat, among each of the modeling
approaches. (D) indicates the model was fit to derived response variables, and * indicates the
intercept/slope parameter fell within the expected range based on equivalence testing. ..............47
Table 2.1: Base model inputs for a single standardized habitat block in the systematic habitat
additions. .......................................................................................................................................69
Table 2.2: Base model inputs for the new treatment units in the Grand River, upstream of the
Sixth Street Dam. Treatment units are ordered top to bottom by increasing upstream distance
from the Sixth Street Dam. ...........................................................................................................76

vi

LIST OF FIGURES

Figure 1.1: Distribution of GLFC sampling locations around the Lake Michigan drainage basin
for the 2004-2008 time period. .....................................................................................................15
Figure 1.2: Distribution of 2016 sampling effort across the Lake Michigan basin, and within each
of the sampled river systems. ........................................................................................................28
Figure 1.3: Distributions of transect-scale values of logit-transformed Type I (a), Type II (b), and
total larval habitat (c) proportions, in addition to the distribution of Type II to Type I habitat
proportion ratios (d). .....................................................................................................................31
Figure 1.4: Plots of the coefficient mean values (points) and 95% HPD intervals (bars), generated
from both the logit and beta linear regressions, for Type I (a) and Type II (b) habitat. ...............36
Figure 1.5: Partial effects plots, generated from the generalized additive model for Type I
habitat, showing the marginal, nonlinear influence of each continuous covariate on the
response. ........................................................................................................................................38
Figure 1.6: Partial effects plots, generated from the generalized additive model for Type II
habitat, showing the marginal, nonlinear influence of each continuous covariate on the
response. ........................................................................................................................................39
Figure 1.7: Observed model fit for the Type I (a, c, e) and Type II (b, d, f) habitat values
generated from the logit regression (a, b), beta regression (c, d), and boosted additive model (e,
f). The line represents the linear fit between observed and predicted values. .............................42
Figure 1.8: Observed model fit for Type I (a, b) and Type II (c, b) transect-scale (a, c) and reachscale (b, d) habitat values generated by the boosted generalized additive model fit to the two
derived response variables. The line represents the linear fit between observed and predicted
values. ...........................................................................................................................................44
Figure 1.9: Boxplots of observed covariate variation among reaches sampled in each river
system. ..........................................................................................................................................46
Figure 1.10: Plots of predictive fit for the logit regression model (a), beta regression model (b),
boosted additive model (c), and boosted additive model fit to the derived habitat quantities (d).
Closed circles and solid lines represent data points and linear fits, respectively, between observed
and predicted values for Type I habitat; open circles and dotted lines represent the same for Type
II habitat. .......................................................................................................................................48
Figure 1.11: Partial effects plots, generated from the generalized additive model for total larval
habitat, showing the marginal, nonlinear influence of each continuous covariate on the
response. ........................................................................................................................................60

vii

Figure 1.12: Partial effects plots, generated from the generalized additive model for Type II to
Type I habitat ratios, showing the marginal, nonlinear influence of each continuous covariate on
the response. ...................................................................................................................................61
Figure 2.1: Map of the Grand River mainstem (a) and the modeled Grand River system between
the Sixth Street Dam and North Lansing Dam (b). Numbers in (b) correspond to numbered
treatment unit names in Table 2.2. ................................................................................................73
Figure 2.2: Visual representation of the modeling scenario combinations considered in
simulating the Sixth Street Dam removal. ....................................................................................78
Figure 2.3: Adult sea lamprey abundance trends with increasing habitat availability, assuming
habitat units are added within a single stream unit. Boxes, whisker bars, and open circles
represent the 25th and 75th, 10th and 90th, and 5th and 95th percentiles of simulated adult
abundances, respectively. Solid horizontal lines and black circles represent corresponding
median and mean values, respectively, and the dashed line represents the exponential model fit to
the mean abundance values. The asterisk indicates mean lamprey abundance from Scenario #1
of the Grand River case study, with an assumed 10% habitat use (see Figure 2.7). Gray squares
indicate proportions of simulations with abundances greater than the status quo 90th
percentile. ......................................................................................................................................81
Figure 2.4: Adult sea lamprey abundance trends with increasing habitat availability, assuming
habitat units are added as independent stream units. Boxes, whisker bars, and open circles
represent the 25th and 75th, 10th and 90th, and 5th and 95th percentiles of simulated adult
abundances, respectively. Solid horizontal lines and black circles represent corresponding
median and mean values, respectively, and the dashed line represents the exponential model fit to
the mean abundance values. The asterisk indicates mean lamprey abundance from Scenario #1
of the Grand River case study, with an assumed 10% habitat use (see Figure 2.7). Gray squares
indicate proportions of simulations with abundances greater than the status quo 90th
percentile. ......................................................................................................................................82
Figure 2.5: Changing model characteristics with increasing habitat. Lines and polygons
represent the median and 10th and 90th percentiles, respectively, across all simulations (a, b) or
treatment units (c). The dotted line and lighter polygon illustrate the effect of adding a single,
large unit, and the solid line and darker polygon illustrate the addition of habitat units as
multiple, discrete treatment units, respectively. ............................................................................85
Figure 2.6: Average annual budget expenditure on the original treatment units with increasing
habitat availability. The dotted and solid lines illustrate the response when habitat is added as a
single, ever-larger system and multiple, discrete treatment units, respectively. ...........................86
Figure 2.7: Expected sea lamprey abundances for each of the management scenarios. Scenarios
#1 and #2 exclude the Looking Glass River, while Scenarios #3 and #4 account for its influence.
New treatment units are treated by the SLCP in Scenarios #1 and #3, and ignored in Scenarios #2
and #4. Boxes, whisker bars, and open circles represent the 25th and 75th, 10th and 90th, and 5th

viii

and 95th percentiles of simulated adult abundances, respectively. The solid horizontal lines and
black circles represent median and mean values, respectively. Numbers above the upper whisker
bars indicate the proportion of simulations greater than the status quo 90th percentile. ...............88
Figure 2.8: Expected sea lamprey abundances when the Sixth Street Dam is removed, the
Webber Dam is modified to block sea lamprey, lampreys are assumed to use 10% (a) or 50% (b)
of maximum potential river length, and the new treatment units are allocated control efforts with
a steadily increasing Lake Michigan control budget. Boxes, whisker bars, and open circles
represent the 25th and 75th, 10th and 90th, and 5th and 95th percentiles of simulated adult
abundances, respectively. The solid horizontal lines and black circles represent median and
mean values, respectively. ............................................................................................................90

ix

THESIS INTRODUCTION

The parasitic sea lamprey (Petromyzon marinus) has played an important contributing
role in the ecology and economics of the Great Lakes since its invasion in the early 20th century.
In conjunction to high fishing pressure and fish habitat loss, peak sea lamprey abundances in the
1950s and 1960s contributed to the declines of numerous fish populations and commercial
fisheries across the Great Lakes; the commercial catch of lake trout (Salvelinus namaycush)
catch declined from 2268 t to 76 t between 1938 and 1954 in Lake Huron, and similar declines
were reported in both Lakes Michigan and Superior (Smith and Tibbles 1980). Other species,
including burbot (Lota lota), lake whitefish (Coregonus clupeaformis), walleye (Sander vitreus),
and rainbow trout (Oncorhynchus mykiss) also experienced declines during this time period
(Smith 1972). In addition to the co-occurring rise and fall of sea lamprey and teleost
abundances, respectively, evidence for the effect of sea lampreys on Great Lakes fish includes
the explicit linking of lamprey wounding rates, mortality rates, and abundance for multiple
species of Great Lakes fish, including lake whitefish and lake trout (Moore and Lychwick 1980;
Spangler et al. 1980; Wells 1980). The removal of top predators also led to longer-term,
persistent changes in lake-wide ecology. Released from predation pressure, multiple species of
native ciscoes (Coregonus spp.) experienced a temporary surge in abundance before shifting
parasitic pressure by sea lamprey caused subsequent depletions in both Lakes Michigan and
Huron (Smith 1968).
The sea lamprey, native to the northern Atlantic Ocean, invaded the upper Great Lakes in
the early 20th century after the creation of the Welland Canal. This artificial waterway
presumably allowed sea lampreys, previously constrained to Lake Ontario, to bypass Niagara

1

Falls and begin spreading throughout the Great Lakes (Hubbs and Pope 1937). The first
evidence of sea lamprey spawning in Lake Erie was confirmed in 1932, and subsequent
observations of sea lamprey in Lakes Huron and Michigan followed in 1937 and 1936,
respectively (Hubbs and Pope 1937; Applegate 1950). The earliest potential observed sea
lamprey in Lake Superior was reported in 1938, with subsequent observations in 1946
(Applegate 1950; Smith and Tibbles 1980). By 1948, sea lampreys were firmly established, with
confirmed spawning runs in 34, 68, and 4 tributaries in Lakes Huron, Michigan, and Superior,
respectively (Applegate 1950).
Sea lampreys are an anadromous species in their native range, and their pattern of
migrating to coastal rivers for spawning persists in the Great Lakes. After the parasitic juvenile
stage finishes feeding within one of the Great Lakes, individuals commence migration to find
suitable rivers for spawning (Applegate 1950). Unlike many anadromous species, sea lampreys
appear to exhibit a lack of natal homing (Bergstedt and Seelye 1995). Stream selection is instead
based on river-specific cues including discharge and the concentration of pheromone-based
larval cues (Moore and Schleen 1980; Morman et al. 1980; Mullett et al. 2003; Sorensen and
Vrieze 2003). Higher river discharge presumably helps migrating individuals sense and orient to
river plumes, while larval cues indicate suitable habitat based on the occurrence of previous
successful spawning (Wagner et al. 2009; Meckley et al. 2014). Upon selecting a river and
moving upstream, adults find suitable gravel habitat for spawning, construct mounded nests,
deposit and fertilize eggs in the nests, and subsequently die (Applegate 1950; Manion and
Hanson 1980). Successfully hatched sea lamprey larvae from these eggs will leave the nest
approximately three weeks after spawning, drift downstream, and locate suitable burrowing
habitat (Applegate 1950).

2

Sea lamprey ammocoetes, or larvae, are largely sedentary and possess a different
morphology from their subsequent life stages. Ammocoetes are relatively unselective filterfeeders that live in fine substrate sediments (Beamish 1980). They lack any semblance of the
oral disc found in the parasitic juvenile stage, and instead use an oral hood to direct water into
the pharynx for respiration and filter feeding (Youson and Potter 1979). Larvae prefer fine
substrate habitat, like silt and sand, for burrowing, and will live and grow within these habitats
for an average of approximately 5 years (Potter 1980). Over the course of their larval stage,
ammocoetes can disperse widely from their original spawning sites; full siblings can be found
more than 1 km apart within one year after hatch (Derosier et al. 2007). Once larvae reach a
threshold length and lipid content, they undergo metamorphosis, in which they develop larger
and more complex eyes, a toothed-oral disc, a rasping tongue, and more prominent dorsal fins,
among other features (Youson and Potter 1979; Potter 1980). Sea lampreys then migrate
downstream to the Great Lakes, where they typically spend 12 to 20 months feeding on teleosts
using their newly developed oral disc and piston-like rasping tongue (Applegate 1950). The
commencement of the spawning migration, after feeding stage is done, then marks the start of a
new cycle.
Since their invasion in the early 20th century, sea lamprey abundances have been
controlled to less than 10% of peak abundances due primarily to the application of sea lampreyspecific pesticides (“lampricides”). Thousands of chemicals were tested by the Great Lakes
Fishery Commission (GLFC) for their selective effects on larval sea lampreys since 1953, and
two chemicals were eventually selected: 3-trifluoromethyl-4-nitrophenol (TFM) and 5,2’dichloro-4’nitrosalicylaniiide (Bayluscide) (Howell et al. 1980). These two chemicals were
applied in conjunction to Great Lakes rivers, starting in 1958, to reduce numbers of larval sea

3

lampreys and subsequently decrease the number of out-migrating parasites. The application of
lampricides has had an observable effect on sea lamprey populations: spawning runs have
become smaller, wounding rates on teleosts have dropped, and abundance of important fish
stocks have rebounded (Smith and Tibbles 1980). Currently, rivers typically are selected
annually for lampricide application based on treatment efficiency; this system is otherwise
known as the Empiric Stream Treatment Ranking (ESTR) system (Christie et al. 2003).
Decisions to treat specific river systems around each Great Lakes basin are guided by
assessment surveys to characterize larval densities and habitat quality. Streams are regularly
surveyed for relative measures of sea lamprey abundance, using the rapid assessment approach,
to rank streams for treatment (Hansen and Jones 2008). Backpack electrofishers are used to
sample optimal larval habitat, and the number of expected larvae greater than 100 mm in length
(“large” larvae) at the end of the growing season is calculated from sampled larvae quantities
using known length frequencies and growth rates. Streams are then ranked according to cost per
expected large larva killed, and the highest ranked streams are selected for lampricide
applications. Larval habitat quantity is also irregularly sampled alongside these surveys. Habitat
is sampled by measuring the proportion of larval sea lamprey-specific habitat types along stream
transects collected around fixed access points and aggregating these habitat measures to the
spatial scale of treatment units (Slade et al. 2003). These measures of habitat can help extrapolate
sampled larval densities to total larval quantities across stream reaches.
Dams play a complementary and well-recognized role in sea lamprey control within the
Great Lakes. Incidental barriers (i.e., those in place before the invasion of sea lamprey) block
sea lamprey from accessing abundant quantities of both spawning and larval habitat in tributaries
(Lavis et al. 2003). The presence of effective barriers has a two-fold benefit for sea lamprey

4

control by reducing their production from streams and eliminating the need for lampricide
treatment in select systems (i.e., freeing up resources for use in other tributaries). An explicit
strategy of the GLFC in sea lamprey control is the construction of new, purpose-built barriers
and maintenance of existing dams in Great Lakes tributaries (Great Lakes Fishery Commission
2011).
Intentional fragmentation in river systems also plays a useful conservation role in systems
other than the Great Lakes. Purposefully restricting aquatic connectivity within watersheds can
yield multiple benefits, including limiting the spread of invasive species, restricting the dispersal
of new diseases, blocking hybridization between wild and hatchery populations, and stopping
organisms from accessing suboptimal habitats (Rahel 2013). Barriers have been modified or
constructed to specifically block the spread of nonnative species like lake trout, brown trout
(Salmo trutta), and rainbow trout in Montana, California, and Arizona, respectively (Avenetti et
al. 2006; D’Angelo et al. 2010; Pister 2010).
Despite these examples for the positive role of barriers, dams can cause more harm than
good to aquatic ecosystems and increasingly are being considered for removal. Correspondingly,
there is growing opportunity and interest in barrier removals across the United States for a
multitude of reasons.

First, there is widespread concern over the public safety of existing dam

infrastructure. Approximately 80% of dams identified by the National Inventory of Dams will
be older than the average life expectancy of dams (i.e., 50 years) by 2020 (NRC 1992; ASDSO
2016). Furthermore, hundreds of hydroelectric dams will require relicensing from the Federal
Energy Regulatory Commission, which necessitates the decision to retain or remove the dam
(Pejchar and Warner 2001). Concurrently, there is increasing motivation to increase aquatic
connectivity and restore fish populations, as the declines of several migratory fish populations

5

have been linked to dams (Limburg and Waldman 2009). These opportunities and motivations,
among numerous others, have driven dam removals efforts in Great Lakes states, including
Michigan. The Michigan Department of Natural Resources has granted millions of dollars to aid
in dam repair and removal, and otherwise assisted in the removal of numerous additional barriers
(Michigan Department of Natural Resources 2016)
Both the benefits and costs of barrier removals in regions like the Great Lakes have to be
considered in decision-making, and these aspects can be informed by fisheries research. Dams
have both positive and negative effects on Great Lakes fisheries due to their role in
simultaneously blocking sea lamprey and other desirable, migratory species fish species (e.g.,
salmonids). The relative magnitudes of these effects, and the expected response of fish
populations upon barrier removals, need to be understood before making dam removal decisions.
The effects of barrier removals on sea lamprey populations represents an especially pertinent
aspect to understand due to their known effects on numerous Great Lakes fish species and a
detailed understanding of their population dynamics. An operating model, incorporating both
sea lamprey population dynamics and control efficacy, has already been developed to formalize
our understanding of this system and represents a potentially useful tool to inform barrier
removal decisions (Jones et al. 2009).
This thesis is intended to characterize the likely response of sea lamprey populations,
including the future effectiveness of sea lamprey control, to barrier removals in Lake Michigan.
Lake Michigan was selected as the focal representative of the Great Lakes based on the detailed
population dynamics information available for the system. In chapter 1, I compare the
performance of multiple modeling approaches in predicting the quantity of sea lamprey larval
habitat upstream of barriers in Lake Michigan tributaries using predominantly catchment-scale

6

landscape attributes. The availability of larval habitat is an important determinant of recruitment
success but is resource-intensive to measure across large spatial extents. Habitat modeling
efforts shed light on dominant influences of substrate-based stream habitat in Michigan and may
inform future population dynamics modeling. In chapter 2, I apply a modified version of a sea
lamprey operating model to evaluate the quantitative response of the Lake Michigan sea lamprey
population to both systematic increases in habitat availability and a specific barrier removal case
study. In doing so, I assess the influence of modeling assumptions, including habitat quality,
habitat use, and alternative management actions, on the expected response, and identify driving
factors behind simulated changes in abundance. Population modeling results revealed
unexpected influences of life history traits and spatial scale on a species’ response to barrier
removals and demonstrate the utility of similar modeling efforts to inform future barrier removal
decisions.

7

CHAPTER 1: COMPARISON OF MODEL PERFORMANCE IN PREDICTING REACHSCALE LARVAL SEA LAMPREY HABITAT ABOVE LAKE MICHIGAN BARRIERS

8

INTRODUCTION

Quantitative measures of stream habitat features, such as those describing substrate
characteristics, are often necessary for estimating potential fish production within rivers. For
example, the amount of larval habitat, physically defined as substrate dominated by fine
particles, is an established driver of Great Lakes sea lamprey (Petromyzon marinus) recruitment
success (Jones et al. 2003; Dawson and Jones 2009; Jones et al. 2009). However, physical
habitat measurements are often unavailable and difficult to collect at the appropriate extent and
resolution for management application, so an alternative is to develop predictive models capable
of estimating habitat quantities in streams.
The use of predictive models to estimate measures of stream habitat can save resources
that otherwise would be required for physical sampling. Stream habitat modeling typically
draws upon the hierarchy theory of streams, in which large-scale catchment characteristics
constrain local physical, chemical, and biotic components of stream, to predict fine-scale habitat
features relevant to aquatic resource management (Curry 1972; Hynes 1975). These established
hierarchies can be very useful in predicting local stream habitat because it can be difficult to
measure fine-resolution habitat features, like sea lamprey larval habitat, across broad spatial
extents (Brenden et al. 2007). Conversely, large-scale catchment predictors, including predictors
like land use and topography, are often readily available in databases like the National Land
Cover Database (Multi-Resolution Land Characteristics Consortium 2015). Many studies have
successfully applied hierarchy theory for local stream habitat predictions, including substrate
composition, using various catchment-scale predictors (Jeffers 1998; Mugodo et al. 2006;
Frappier and Eckert 2007; Brenden et al. 2007; Wang et al. 2013).

9

The quantification of sea lamprey larval habitat in stream reaches above dams in the
Great Lakes basin is a particularly important issue that can be informed by predictive modeling.
Sea lampreys in the Great Lakes are an invasive species predominantly controlled by the
application of selective pesticides, specifically 3-trifluoromethyl- 4’-nitrophenol (TFM) and 2’,
5-dichloro-4’- nitrosalicylanilide (Bayluscide), to streams to kill the resident larval sea lamprey
(Smith and Tibbles 1980). Dams also play a key role in sea lamprey control by blocking
lampreys from sections of stream that would otherwise provide additional spawning and larval
habitat area (Hunn and Youngs 1980). These same dams, however, are being considered for
removal to improve public safety and enhance fishery resources. Removal of dams could greatly
increase the amount of available habitat (Lavis et al. 2003) and possibly decrease the
effectiveness of current control efforts if funds available for lampricide treatments are
constrained. Knowledge of sea lamprey population responses to dam removals is needed to
allow more explicit cost-benefit analyses to inform decisions on barrier removals. Quantifying
the amount of available larval habitat above dams, as an important driver of recruitment success,
is needed to predict the overall response of sea lamprey populations to dam removals.
The Great Lakes Fishery Commission (GLFC) has conducted habitat sampling for larval
sea lamprey in currently infested stream reaches, in which they quantify the proportion of
different habitat categories based on substrate composition along stream transects (Larval
Assessment Work Group 2013). Larval habitat is divided into three primary categories (Type I,
II, and III), among which only Type I and II habitat are considered suitable for larvae due to the
presence of fine sediments for burrowing. More specifically, Type I, II, and III habitat are
defined as mixtures of fine sands and organic matter, medium to coarse shifting sands, and hardpacked surfaces (e.g., bedrock, gravel, hardpan clay), respectively (Applegate 1950; Slade et al.

10

2003). Type I habitat was observed to contain 93% of sampled sea lamprey larvae in Great
Lakes tributaries and has been found to supports significantly higher larval densities than Type II
(Mullett 1997). Because sea lamprey are not present in stream reaches above effective lamprey
barriers, little sampling has been conducted by the GLFC to characterize larval habitat in these
areas. Although previous studies have attempted to predict larval sea lamprey habitat using a
combination of river and landscape features (Neeson et al. 2007; Neeson et al. 2012), none have
developed models with the intent to predict habitat in new spatial extents, using landscape-scale,
widely available GIS-derived covariates.
Large-scale landscape predictors and fine-scale habitat elements are measured at very
different scales; hierarchical models are particularly well-suited to handling these spatial
discrepancies. Hierarchical models, also known as multilevel models, are a generalization of
linear regressions, in which model coefficients themselves are modeled with distinct parameters
(Gelman 2006a). The incorporation of random effects in these models, in which spatial scales
can be treated as nested levels, allow for prediction outside of the existing spatial units or levels;
this is a valuable feature for predicting habitat in areas outside of the data extent. Random
effects represent “a random sample of a larger set of potential effects” instead of representing
“all possible levels of a factor” (Wagner et al. 2006). These hierarchical models can be flexibly
implemented in a Bayesian framework, which also allow for probabilistic interpretation of
results (Midway et al. 2014). Multilevel models are therefore suitable for modeling transectlevel sea lamprey habitat features based on reach-level landscape predictors.
Proportional data can also be difficult to model using traditional modeling techniques.
Sea lamprey habitat data are reported as proportions; a logit transformation is commonly used to
re-scale proportion data to allow estimation within a linear modeling framework (Warton and

11

Hui 2011). An extension of the logit transformation is beta regression, based on the beta
distribution; this approach is suited to model response variables that take values between 0 and 1
(Ferrari and Cribari-Neto 2004). With a beta regression, a link function (typically logit, but
sometimes probit or cloglog) connects the response data to the explanatory variable, which are
then modeled assuming a beta distribution (Liu and Kong 2015) An extension of the beta
regression, namely zero/one inflated beta regression, is useful for cases where zeroes and ones
are common in the response variable data.
Additional problems in predictive modeling include incorrectly assuming linear
relationships exist between dependent and independent variables and overfitting from an excess
of predictor variables. Nonlinearity and thresholds can be expected in fluvial geomorphology,
wherein stream features are simultaneously affected by interacting factors at a multitude of
spatial and temporal scales (Shreve 1979). Generalized additive modeling using componentwise gradient boosting is an emerging method to avoid overfitting data and allow nonlinear
relationships among predictor and response variables (Hofner et al. 2014). Additive models,
introduced by Hastie and Tibshirani (1986), allow for smooth, nonlinear, component-wise fitting
of predictors to response data using specified linear components. Pairing additive modeling with
component-wise gradient boosting further allows for implicit variable selection by iteratively
improving the model one predictor variable at a time (Friedman 2001; Hofner et al. 2014). This
modeling approach was successfully applied in Wang et al. (2013) to predict multiple facets of
fine-scale stream habitat, including substrate composition.
I used larval habitat data collected by the GLFC between 2004 and 2008 and a wide
range of potential landscape-scale covariates attributed to stream reaches to predict the amount
of Type I and Type II larval habitat (i.e., suitable burrowing habitats) in stream reaches upstream

12

of existing dams. I first evaluated the ability of three modeling approaches, fit independently to
the original Type I and Type II habitat proportion response data, to predict both habitat types by
comparing their model fit, important predictor variables, and predictive capacity: 1) hierarchical
linear regression in a Bayesian framework using logit-transformed response variables, 2)
hierarchical zero/one inflated beta regression in a Bayesian framework using the raw
proportional response variables, and 3) generalized additive modeling, in a Gaussian likelihood
framework, with component-wise gradient boosting using the logit-transformed response
variables. The predictive capacity of each modeling approach was assessed using independent,
out-of-sample data collected in Lake Michigan tributaries, upstream of existing dams, in the
summer of 2016. Following these comparisons, I then evaluated the ability of the generalized
additive model to predict Type I and Type II habitat when the model was fit to two new, derived
response variables: the total larval habitat proportion (i.e., the sum of Type I and Type II habitat
proportions) and the ratio between the two habitat types. Combining the two new predicted
response variables represents a unique way to predict individual habitat quantities, in which
information is shared across the two response variable datasets.
The successful quantification of larval habitat in stream reaches upstream of relevant
barriers can inform simulations of sea lamprey response to dam removals and has the potential to
subsequently improve barrier removal decisions. Furthermore, an improved understanding of the
driving landscape influences on sea lamprey habitat can guide future habitat sampling and
estimation.

13

METHODS

Historical Data
Larval habitat data were collected by GLFC field crews between 2004 and 2008
according to their current larval assessment sampling protocol (Larval Assessment Work Group
2013). Within each treatment unit, defined by the GLFC as a section of stream where a unified
lampricide treatment can be applied, field crews sampled up to six access points, which were
simply points along the river where crews could physically access the water. Sampling locations
were broadly distributed around sea lamprey-infested watersheds in the Lake Michigan drainage
basin (Fig. 1.1). At each access point, control agents sampled two upstream and two
downstream cross-sectional transects, and divided the total transect length of each into segments
of one of the following three aquatic habitat types, identified using visual and tactile assessment:
Type I, Type II, and Type III. The proportion of Type I and Type II habitat along each transect
was calculated from these distances and used as the primary response variable in this analysis.
The GLFC uses total larval habitat areas, calculated from these surveys as the product of
treatment unit length, average width, and average habitat proportions, in combination with larval
abundance surveys to estimate total larval abundance within treatment units and inform control
efforts.

14

Figure 1.1: Distribution of GLFC sampling locations around the Lake Michigan drainage basin
for the 2004-2008 time period.

15

All landscape covariate data were obtained from a spatial framework developed by Wang
et al. (2011), in which landscape and stream network variables were assigned to stream reaches,
as well as their associated catchments and buffers, across the conterminous United States using
the National Hydrography Dataset Plus (1:100,000 scale; http://www.horizonsystems.com/nhdplus). A wide range of variables were selected to capture aspects of climate,
elevation, geology, soil, land cover, groundwater contribution, as well as river size and
connectivity. These variables were attributed to five possible spatial scales: the reach, local
buffer (LB), network buffer (NB), local catchment (LC), and network catchment (NC). Within
this spatial framework, a reach is defined as an inter-confluence section of a stream, a buffer as
the portion of the landscape within a set distance of the reach (90 m in this case), and a
catchment as the entirety of upstream land draining into the lowermost portion of a reach
(Brenden et al. 2006). Local and network refer to how the buffers and catchments are bounded;
local buffers and catchments only extend to the boundaries with the nearest neighboring reaches,
and network buffers and catchments are the totality of all buffer or catchment area upstream of
the reach in question (Brenden et al. 2006).
Because measures of the predictor variables were attributed to NHDPlus-defined reaches
and not those defined by the GLFC, the habitat measurements associated with transects were reassigned to NHDPlus-defined reaches using GPS coordinates recorded for each
downstream/upstream transect pair. This work was performed in ArcMap (PC ARC/GIS
Version 10.0. Environmental System Research Institute, Redlands, California,
http://www.esri.com/software/arcgis).

16

Landscape Variable Reduction Techniques
Variable reduction of the initial candidate set of approximately 300 variables included the
following steps: 1) remove variables with no observed variation, 2) aggregate detailed
categorical variables (e.g., attributes like land cover, surficial lithology) into more commonlyused, functional predictors, 3) remove variables with no proposed mechanistic relationship to
stream sediment composition, 4) choose just one spatial scale per variable where applicable (LB,
NB, LC, NC), and 5) remove collinear variables. The selection of combined variable categories,
removal of variables with no proposed relationships to sediment composition, and choice of
spatial scale were based on available scientific literature, professional judgment by experts in the
field of ecohydrology, and best personal judgment. Assessment of collinearity was conducted
using a generalized variance inflation factor (GVIF) analysis applied to the logit-transformed
Type I and Type II response variables (Fox and Monette 1992), where modified GVIF values
(GVIF1/2p) greater than 3, as a proxy for VIF, indicated the presence of troublesome collinearity
(Zuur et al. 2010). GVIF was used in place of the more traditional VIF due to the presence of
categorical variables in the analysis. The goal of variable reduction was to obtain a manageable
set of uncorrelated covariates, supported by personal and professional judgment, to predict sea
lamprey habitat.

Pre-Modeling Data Transformations
Prior to fitting the hierarchical logit regression and generalized additive models to the
individual Type I and Type II responses, I transformed the proportion data for both habitat types
using the logit-transformation recommended by Warton and Hui (2011). Because the log-based

17

transformation cannot inherently handle the presence of zeroes and ones in the dataset, the
following modified logit-transformation was used:
(𝑥+𝑀)

𝑙𝑜𝑔 (1−𝑥+𝑀)
where x is the proportional response data and M is the minimum non-zero value in the response
dataset being transformed. In contrast, untransformed proportion data were used in fitting the
zero/one-inflated beta regression model.
The same modified logit-transformation was applied to the transect-specific total larval
habitat proportions, which were simply the sum of transect-specific Type I and Type II habitat
proportions. The ratio of Type II to Type I habitat proportions was calculated using a similar
modification to handle zeroes and ones; the minimum non-zero value in each response dataset
was added to the Type I and Type II habitat proportions before the quotient calculation.
All landscape covariates were centered and scaled to increase the ease of model
convergence and to facilitate coefficient comparison within and among models (KĂŠry 2010).

Fitting the Original Type I and Type II Response Data
Hierarchical Logit Regression Model
I constructed this model with three levels in the hierarchy: transect-scale habitat
proportion measurements, the random effect of access point, and the random effect of reach. I
included the effect of access point to account for the fact that transects were not randomly
sampled throughout each reach; instead, transects were systematically sampled within “clusters”
(i.e., around the access points). The random effect of access point helps accounts for the
expected higher similarity among transects within these “clusters” than among transects from
separate clusters. Access point was nested within reach to match the spatial scale of the

18

landscape covariates, and the linear regression of the thirteen covariates was applied to the
random effect of reach (variables are defined in Table 1.1):
Hierarchy Level 1 (Transect):
𝑦𝑖𝑗𝑘 ~ 𝑁𝑜𝑟𝑚𝑎𝑙 (𝛼1 𝑗𝑘 , 𝜎1 2 )
Hierarchy Level 2 (Access Point):
𝛼1𝑗𝑘 ~ 𝑁𝑜𝑟𝑚𝑎𝑙(𝛼2 𝑘 , 𝜎2 2 )
Hierarchy Level 3 (Reach):
𝛼2 𝑘 ~ 𝑁𝑜𝑟𝑚𝑎𝑙(𝜇1 𝑘 , 𝜎3 2 )
𝜇1 𝑘 = 𝛽1 ∗ 𝑃1 𝑘 + 𝛽2 ∗ 𝑃2 𝑘 + . . . + 𝛽𝑛 ∗ 𝑃𝑛 𝑘 + 𝛼𝑝 + 𝛼3
The model was constructed in WinBUGS (Lunn et al. 2000) and run from R (R Core
Team 2016) using the R2WinBUGS package (Sturtz et al. 2005, Running ‘WinBUGS’ and
‘OpenBUGS’ from ‘R’ / ‘S-PLUS’, version 2.1-21, https://cran.r-project.org/package=coda).
The full WinBUGS script for hierarchical Bayesian model used to predict Type I habitat is
provided in Appendix A.

19

Table 1.1: Descriptive key of model parameters describing the Bayesian hierarchical logit
regression and beta regression models.
Variable
i
j
k
y
𝛼1
𝛼2
𝛼3
𝜇1
𝛽
𝑃
n
𝛼𝑝
p
𝜎12
𝜎22
𝜎32
𝑖
𝑗
𝑦
𝑝𝑖𝑗
𝑞𝑖𝑗
𝛼𝑖𝑗1 , 𝛼𝑖𝑗2
𝜇𝑖𝑗
𝑣𝑖𝑗
𝛽
𝑃
n
𝛼𝑝
p
𝛼1 , 𝛼3 , 𝛼4
𝛼2𝑗
𝜇1
𝜎12

Description
Logit Regression
Transect-scale
Access point-scale
Reach-scale
Logit-transformed transect-scale habitat measurement
Random effect of access point, nested within reach
Random effect of reach
Linear regression intercept
Mean reach effect
Slope coefficient for continuous predictor n
Value for continuous predictor
Number of continuous predictors
Coefficient for dam presence
Categorical predictor levels
Unexplained model variance
Within access point variance
Within reach variance
Beta Regression
Transect-scale
Reach-scale
Transect-scale habitat measurement
Probability yij=0
Conditional probability yij=1|yij≠0
Beta distribution shape parameters
Beta distribution mean
Beta distribution variance
Slope coefficient for continuous predictor n
Value for continuous predictor
Number of continuous predictors
Coefficient for categorical predictor
Categorical predictor levels
Regression intercept values
Random effect of reach
Mean reach effect
Within reach variance

20

In running this model, I used uninformative normal priors (N ~ [0, 1000]) for all slope,
fixed effect, and intercept terms. Similar priors have been used in other hierarchical, linear
modeling approaches (Thomas et al. 2006; Midway et al. 2015). Vague uniform priors were
applied for the three standard variance terms (Unif ~ [0, 20]), as recommended by Gelman
(2006b) for hierarchical models.
Recommended chain length was assessed by running a 1000 iteration chain with no burnin and applying Raftery and Lewis’ diagnostics. Raftery and Lewis’ convergence diagnostics
indicated the necessary chain and burn-in length for each evaluated model parameter to achieve a
certain probability of obtaining a model estimate within some margin of error around a set
quantile; I specified a 95% probability, 5% margin of error, and 50% quantile in my analyses
(Cowles and Carlin 1996). Following the burn-in and chain length recommendations from the
Raftery and Lewis’ diagnostics, I ran the full model and assessed model convergence using
Geweke’s criterion, in which the absolute value of a z score greater than 1.96, from the
comparison between the means of the first 10% and last 50% of the MCMC chain for every
evaluated parameter, indicated a failure to converge at the p=0.05 level of significance (Geweke
1992). I also performed visual examinations of the parameter trace and autocorrelation plots. I
evaluated the posterior probability distribution for every intercept, fixed effect, and slope
parameter, in addition to all variance terms and reach-scale model predictions. MCMC chain
diagnostics were assessed using the coda package (Plummer et al. 2015, Output Analysis and
Diagnostics for MCMC, version 0.18-1, https://cran.r-project.org/package=coda).
Hierarchical Zero/One Inflated Beta Regression Model
I intended the hierarchical zero/one inflated beta regression to account for the influence
of extreme response variable values (i.e., zeroes and ones), model the response using a non-

21

normal distribution, and account for the random effect of reach. In fitting the zero/one inflated
beta regression for the untransformed proportion data, I included the same set of covariates as
those applied in the hierarchical logit regression model, as well as the random effect of reach. I
also allowed for separate intercept terms to account for the zeroes and ones via zero/one
inflation. Unlike the logit regression, the random effect of access point was not accounted for
and predictors were fit to response data at the transect scale; these changes were due to inherent
limitations in the R package used to fit the model. The model equations for this beta regression
with zero/one inflation and the random effect of reach, modified from Liu and Kong (2015), are
provided below (variables are defined in Table 1.1):
Hierarchy Level 1 (Transect):
𝑝𝑖𝑗

𝑖𝑓 𝑦𝑖𝑗 = 0

𝑦𝑖𝑗 = { (1 − 𝑝𝑖𝑗 )𝑞𝑖𝑗
(1 − 𝑝𝑖𝑗 )(1 − 𝑞𝑖𝑗 )𝐵𝑒𝑡𝑎(𝛼𝑖𝑗1 , 𝛼𝑖𝑗2 )

𝑖𝑓 𝑦𝑖𝑗 = 1
𝑖𝑓 𝑦𝑖𝑗 ∈ (0,1)

𝜇𝑖𝑗 (0,1) = 𝛼𝑖𝑗1 (𝛼𝑖𝑗1 + 𝛼𝑖𝑗2 )−1
𝑣𝑖𝑗 = 𝛼𝑖𝑗1 + 𝛼𝑖𝑗2
𝑙𝑜𝑔𝑖𝑡(𝜇𝑖𝑗 (0,1) ) = 𝛽1 ∗ 𝑃1𝑗 + 𝛽2 ∗ 𝑃2𝑗 + . . . + 𝛽𝑛 ∗ 𝑃𝑗𝑛 + 𝛼1 + 𝛼𝑗𝑝 + 𝛼2𝑗
𝑙𝑜𝑔𝑖𝑡(𝑝𝑖𝑗 ) = 𝛼3
𝑙𝑜𝑔𝑖𝑡(𝑞𝑖𝑗 ) = 𝛼4
Hierarchy Level 2 (Reach):
𝛼2𝑗 ~𝑁𝑜𝑟𝑚𝑎𝑙(𝜇1 , 𝜎12 )
Default diffuse normal priors (N ~ [0,1000]) were used for all model parameters except
for the variance terms; uniform priors with the model’s default distribution (Unif ~ [0, 20]) were
used for the variance parameters. Recommended chain length and initial convergence were

22

again assessed using Raftery and Lewis’ convergence diagnostics run on two separate 1000
iteration chains with no burn-in. Following the burn-in and chain length recommendations from
the Raftery and Lewis’ diagnostics, I assessed final model convergence using Geweke’s criterion
and examination of the trace and autocorrelation plots.
The model was run from R (R Core Team 2016) using the zoib package (Liu and Kong
2015, Bayesian Inference for Beta Regression and Zero-or One Inflated Beta Regression, version
1.3.3, https://cran.r-project.org/package=zoib), which runs the MCMC chain using JAGS
software (Plummer 2014).
Boosted Generalized Additive Model
For this model, I assumed a Normal distribution for the logit-transformed habitat values,
and correspondingly applied a normal negative log-likelihood loss function within the model
framework. The same landscape predictors as in the previous models were applied in the
boosted generalized additive model; all variable reduction was conducted during the iterative
model fitting. No random effects were included due to the model’s inability to separate random
and fixed effects in the implicit variable reduction. The continuous covariate effects were
incorporated using additive P-spline base-learners (typical in generalized additive modeling)
with the default 4 degrees of freedom and 20 interior knots; categorical effects were modeled
using ordinary least squares base-learners (Schmid and Hothorn 2008; Wang et al. 2013).
Following the recommendations of Hofner et al. (2014), I determined the optimal number of
boosting iterations (mstop) in model fitting by fitting the full model (mstop=100), applying 10-fold
cross-validation, and identifying the mstop value that minimized predictive risk. The model was
run from R (R Core Team 2016) using the mboost package (Hofner et al. 2016, Model-Based
Boosting, version 2.6-0, https://cran.r-project.org/package=mboost).

23

Fitting the Derived Response Variables
I fit the same landscape predictors for Type I habitat models to the total larval habitat
proportions and ratio response variables within a generalized additive model framework. For the
total larval habitat response, I again assumed a Normal distribution for the logit-transformed
values and applied a normal negative log-likelihood loss function. Based on the observed
distribution of the ratio between the two habitat types, I decided to fit a Laplace (i.e., double
exponential) loss function to best characterize the ratio-based response variable. The Laplace
distribution exhibits exponential decline moving away from zero in a positive direction and, as
such, can be successfully applied to long-tailed error distributions (Hofner et al. 2014). The
same base-learners and procedure to identify an optimal number of boosting operations, as
described for fitting the individual habitat response data, were applied in these two models.
After fitting the model and generated predictions for derived quantities, individual predictions
for Type I and Type II habitat were calculated by dividing the predicted total larval habitat
proportion by the appropriate Type II to Type I ratio estimate.

Assessing Variable Influence and Model Fit
Variable importance and direction of covariate influence for the logit regression and
zero/one inflated beta regression models were assessed by examining the coefficient posterior
distribution, and specifically the 95% highest posterior densities (HPD). Highest posterior
density intervals reflect the smallest range of coefficient values that contains a specified
percentage of the coefficient estimates generated from the MCMC chain, and can be applied to
determine variable importance in a similar manner to frequentist confidence intervals. I chose to
interpret coefficients with a 95% HPD that did not overlap zero as a strong predictor of larval

24

habitat. The coefficient sign (positive/negative) determined the direction of variable influence.
For the boosted generalized additive models, I assessed each variable’s relative
importance in driving model fit using selection frequencies, representing the number of iterations
in the model fitting process in which each variable was selected as the best predictor. The
direction of variable influence was assessed through the use of partial effects plots, which
illustrate the nonlinear effect of each variable on the response with all other variables held at a
constant value.
Model fit was assessed by comparing observed and predicted habitat values, on both the
transect- and reach-scale larval habitat proportions. The transect scale is the original scale of our
models’ response variables, and assessing model fit on this scale reflects the realized amount of
variation among transects within each reach. Model fit was also assessed on the reach scale, as
this the scale of interest for future habitat predictions; observed reach-scale habitat proportions
were obtained by calculating the simple average of all available transect-scale observations
within each reach. For the hierarchical logit regression and zero/one inflated beta regression
models, habitat predictions were generated using just the mean of the model-estimated
coefficients and the original model data. Random effects were not incorporated because future
habitat predictions, in stream access points or reaches outside of those used to fit the models, can
only be calculated using these fixed coefficients. Habitat predictions for the boosted generalized
additive model were obtained using the package’s built-in prediction functionality.
Quantification of model fit for both habitat types was achieved by fitting a linear
regression model between observed and predicted values, and calculating the intercept, slope,
and model fit (R2; equivalent to the squared Pearson’s correlation coefficient between the two
datasets). I also performed regression-based equivalence tests to assess whether intercept and

25

slope terms from a linear regression fit to the observed and predicted data fell within some
specified range around the expected value (Robinson et al. 2005). The intercept test specifically
assessed whether the means of the observations and predictions are equivalent; the slope test
assesses whether the slope term is equal to the expected value of 1. I chose a relatively broad
equivalence range of 50% due to the large degree of stochasticity and dynamism typical in instream habitat features; a similar value was applied in Wang et al. (2013). Finally, based on the
recommendation of PiĂąeiro et al. (2008), I calculated the root mean squared deviation (RMSD)
for the association between the observed and predicted larval habitat proportions; lower RMSD
values indicate better model fit. The RMSD correctly assessed the deviation of predicted values
from the expected 1:1 line between observed and predicted values, as opposed to similar
measures like mean squared error (MSE).

Evaluating Model Predictive Ability on Out-of-Sample Data
As the modeling results of this study were intended to predict habitat quantities in stream
reaches upstream of dams, and more generally in spatial extents not included in the model
building dataset, I conducted habitat sampling in stream reaches upstream of existing dams in
order to create a dataset for external model validation. This dataset allowed me to assess
whether the models can successfully predict in these upstream systems, the habitats of interest
for this study. Sampling occurred between July 18, 2016 and August 11, 2016.
Sampling sites were determined by first selecting river systems in the Lake Michigan
with the following characteristics: (a) contained a sea lamprey barrier blocking an extensive
amount of potential lamprey habitat, (b) were previously sampled sometime between 2004-2008
downstream of the dam, (c) were listed as a high priority by sea lamprey control agents in the

26

GLFC’s Barrier and Larval Control Task Forces, and (d) represented the diverse spatial extent
captured by the original model dataset. A diverse spatial extent was crucial to capture the full
range of covariate variability within the model training dataset. Six streams were eventually
selected for sampling: Twin River, Peshtigo River, Whitefish River, Boardman River, Pentwater
River, and Kalamazoo River. The Twin and Peshtigo Rivers are in Wisconsin, the Whitefish
River is in Michigan’s Upper Peninsula, and the remaining three are in Michigan’s Lower
Peninsula (Fig. 1.2).

27

1

3
2
4
1
2
5

6

3

4

6

5

Figure 1.2: Distribution of 2016 sampling effort across the Lake Michigan basin, and within each
of the sampled river systems.

28

Based on resampling simulations of the 2004-2008 transect data (A. Jensen, Michigan
State University, East Lansing, Michigan, unpublished analysis), four access points were
determined to be desirable to capture the true expected amount of Type I and Type II habitat
within a given reach. This desired number of access points was selected by qualitatively
assessing the trade-off between minimizing the observed coefficient of variation for observed
reach-scale Type I and Type II habitat quantities, assuming that four transects were taken at each
access point, and sampling a maximum number of stream reaches.
For each of the six river systems, I aimed to sample a minimum of eight randomly
selected NHD reaches, upstream of the first effective lamprey barrier, with each reach having at
least four accessible road crossings. Low crossing densities and inaccessible stream crossings in
some systems necessarily resulted in stream reaches with less than four access points; a
minimum of two access points and eight transects were required to include reach data in model
validation. Furthermore, some reaches, when examined in the field, were stagnant (lentic), dry,
or otherwise inaccessible; when feasible, these reaches were replaced by opportunistic sampling
of nearby stream reaches.
The in-stream sampling protocol followed that of the GLFC Sea Lamprey Control Board
(Larval Assessment Work Group 2013), with one primary exception. If either upstream or
downstream transects could not be measured at an access point, compensatory additional
transects were collected in the accessible direction; a maximum of four transects were collected
at any one access point. I then used these data to assess the ability of each modeling approach to
predict reach-scale Type I and Type II habitat quantities using available landscape covariate data.
The same metrics used to evaluate model fit (e.g., RMSD, equivalence testing, R2) were used to
compare predictive capacity.

29

RESULTS

Data Transformations
The logit transformation established approximately normal distributions in the larval
habitat quantity variables between extreme high and low values, but did not account for the
dominance of these extreme values on both ends of the distribution (e.g., zeroes and ones in the
original proportional data) (Fig. 1.3). The habitat ratio response had a strongly right-skewed
distribution (Fig. 1.3d).

Variable Reduction
At each of the four covariate spatial scales, I aggregated 15 land cover and surficial
lithology variables into just 5 covariates following recommendations of local experts in Lake
Michigan limnology (Table 1.2; D. Infante, Michigan State University, East Lansing, MI,
personal communication, 2015; D. Lusch, Michigan State University, East Lansing, MI, personal
communication, 2015). Following all other variable reductions, sixteen variables were finally
considered for model inclusion prior to analysis of collinearity; this analysis resulted in the
removal of three covariates for each habitat type (Table 1.3). The final variable set for each
habitat type therefore consisted of the following thirteen variables: agriculture land cover,
impervious land cover, natural vegetation land cover, wetland land cover, catchment area,
catchment elevation, catchment slope, lithology (fine or coarse for Type I and II habitat,
respectively), precipitation, soil permeability, groundwater index, road crossing density, and dam
presence.

30

Frequency

a)

b)

Habitat Quantity (Log-Odds)

Habitat Quantity (Log-Odds)

c)

d)

Habitat Quantity (Log-Odds)

Habitat Ratio

Figure 1.3: Distributions of transect-scale values of logit-transformed Type I (a), Type II (b), and
total larval habitat (c) proportions, in addition to the distribution of Type II to Type I habitat
proportion ratios (d).

31

Table 1.2: Summary of aggregated class-based predictor variables considered for model
inclusion.
Aggregated Variable
% Agriculture

Constituent Variables
Pasture/Hay
Cultivated Crops

Justification
Agricultural land is expected to
serve as source of fine sediments

% Natural Vegetation

Deciduous Forest
Evergreen Forest
Mixed Forest
Shrubland - Shrub/Scrub
Herbaceous / Grassland

Riparian buffer zones are known
to filter run-off, including
sediments

% Wetland

Woody Wetlands
Emergent Herbaceous Wetlands

Wetlands, similar to riparian
zones, can filter run-off and serve
as a sediment trap

% Fine Lithology

Glacial Till - Clayey
Glacial Till - Loamy
Glacial Lake Sediment - Fine

Small substrate lithology is
expected to serve as a source of
fine sediments and inhibit
groundwater loading

% Coarse Lithology

Glacial Outwash - Coarse
Glacial Lake Sediment - Coarse
Eolian Sediment - Coarse

Coarse substrates are expected to
serve as a source of coarser
sediments and promote
groundwater loading

32

Table 1.3: Final set of predictors considered for modeling. Bolded terms indicate those variables
not considered collinear according to the GVIF analysis, and therefore included in the final
modeling efforts. LB, NB, LC, NC, and reach refer to local buffer, network buffer, local
catchment, network catchment, and reach scales, respectively. Further details on the resolution
and year of creation for these predictor datasets can be obtained in Wang et al. 2011.
Variable
Agricultural Land Cover
Impervious Surface Land Cover
Urban Land Cover
Wetland Land Cover
Catchment Area
Mean Annual Precipitation
Upstream Dam Presence within 100 km
Catchment Slope
Mean Elevation
Coarse Surficial Lithology (Type II)
Fine Surficial Lithology (Type I)
Soil Permeability Rate
Groundwater Contribution to Baseflow
Road Crossing Density by Stream Length
Natural Vegetation Land Cover
Stream Order
Stream Type

Unit
%
%
%
%
km2
mm yr-1
Binary (Y / N)
°
m
%
%
100*(inches hr-1)
%
# km-1
%
Categorical
Categorical

33

Scale
NC
NC
NC
NC
NC
NC
NC
LC
LC
LC
LC
LC
LC
LC
LB
Reach
Reach

Final Dataset Size
Only response data with corresponding data for all thirteen independent variables were
included in the final modeling efforts. This resulted in 4153 transect-scale measurements for
both Type I and Type II habitat. These 4153 transects represented a total of 941 unique access
points and 428 NHD stream reaches; these groupings were used in the hierarchical modeling and
the assessment of model fit at the stream reach-scale.

Convergence Diagnostics
Based on both the initial and final convergence diagnostics, the logit and beta regression
models for both habitat types appeared to converge satisfactorily. The Raftery and Lewis’
diagnostics suggested varying chain lengths across the habitat types and models, but consistently
low recommended burn-in values. To ensure adequate model convergence, I applied the same
10,000 iteration chain with a 1,000 iteration burn-for Type I and Type II habitat in the logit
regression model. Because model computation took significantly more computing power and
time, I was less liberal in total chain length for the zero/one inflated beta regression model: I
chose to run a 5,000 iteration, 500 burn-in model for Type I habitat and an 8,000 iteration, 500
burn-in model for Type II habitat.
For the resulting final models, Geweke’s convergence diagnostics for all intercepts, slope
coefficients, and variance terms indicated that all parameters converged satisfactorily; the worst
value across both modeling approaches and habitat types was -1.764, safely below the threshold
value of 1.96. There were no apparent issues with autocorrelation, as the observed degree of
autocorrelation for all model variables decreased either immediately or regularly with increasing

34

lag; examinations of parameter trace plots also revealed reasonably consistent mixing for all
parameters over the chain length.

Variable Coefficients
I used the distribution mean and 95% highest posterior density (HPD) intervals to
evaluate the relative influence of each covariate coefficient among the logit and beta regression
modeling approaches and two habitat types (Fig. 1.4). I selected the mean as the measure of
central tendency based on the highly symmetric parameter distributions and associated 95% HPD
intervals. Among the coefficients for both Type I and Type II habitat, there was a noticeable
difference between those generated from the logit and beta linear regressions; coefficients
generated from the beta regressions tended to have less extreme mean values and narrower 95%
HPD intervals (Fig. 1.4). This difference resulted in different predictors considered important
from each modeling approach. For Type I habitat, based on the 95% HPD intervals, the logit
linear regression indicated that six covariates (agriculture land cover, natural vegetation land
cover, catchment area, fine lithology, groundwater index, and dam presence) were strong
predictors; only four (natural vegetation land cover, catchment area, fine lithology, and dam
presence) had 95% HPD intervals that did not overlap with zero for the beta regression. For
Type II habitat, logit regression indicated that eight covariates (agriculture land cover, wetland
land cover, catchment area, catchment slope, coarse lithology, precipitation, groundwater index,
and road crossing density) were important predictors; beta regression indicated that seven
covariates (catchment area, catchment elevation, catchment slope, coarse lithology, precipitation,
soil permeability, and groundwater index) were strong predictors. The directions of covariate
influence, however, were highly similar between the two modeling approaches.

35

Covariate

a)

a)

b)

Coefficient Value
Figure 1.4: Plots of the coefficient mean values (points) and 95% HPD intervals (bars), generated
from both the logit and beta linear regressions, for Type I (a) and Type II (b) habitat.

36

Selection frequencies, which served as the best available measure of variable importance
for the boosted additive model, indicated that catchment elevation and mean precipitation were
the most utilized predictors for Type I and Type II habitat, respectively (Table 1.4). Catchment
area and precipitation were the most influential predictors for the derived response variables of
total larval habitat and Type II to Type I habitat ratios, respectively. The directional influence
of the landscape covariates for both habitat types were predominantly nonlinear, with the
apparent exceptions being natural vegetation land cover and fine lithology for Type I habitat and
wetland land cover, impervious land cover, and coarse lithology for Type II habitat (Figs. 1.5,
1.6). Similar nonlinear trends were observed for the two derived response variables (Appendix
B).

Table 1.4: Variable selection frequency in the boosted generalized additive model fit to Type I
habitat, Type II habitat, total larval habitat, and the Type II to Type I habitat ratio. The bolded
values indicate the three most frequently selected predictor variables for each response.
Variable
Agriculture Land Cover
Impervious Land Cover
Natural Vegetation Land Cover
Wetland Land Cover
Catchment Area
Catchment Elevation
Catchment Slope
Coarse Lithology
Fine Lithology
Precipitation
Soil Permeability
Groundwater Index
Road Crossing Density
Dam Presence

Selection Frequency
Type II
Total Habitat
0.124
0.108
0.050
0.028
0.024
0.026
0.074
0.092
0.144
0.176
0.036
0.034
0.030
0.128
0.012
N/A
N/A
0.048
0.208
0.144
0.122
0.108
0.050
0.126
0.050
0.046
0
0.012

Type I
0.068
0.058
0.038
0.135
0.100
0.171
0.118
N/A
0.058
0.086
0.048
0.028
0.078
0.012

37

Ratio
0.006
0.038
0.052
0.032
0.200
0.110
0
N/A
0.018
0.318
0.010
0.182
0
0.034

fpartial

Standardized Covariate Value
Figure 1.5: Partial effects plots, generated from the generalized additive model for Type I
habitat, showing the marginal, nonlinear influence of each continuous covariate on the response.

38

fpartial

Standardized Covariate Value
Figure 1.6: Partial effects plots, generated from the generalized additive model for Type II
habitat, showing the marginal, nonlinear influence of each continuous covariate on the response.

39

Model Fit
Across the three models fitting the original Type I and Type II habitat data, none of the
approaches predicted Type I habitat proportions well (Table 1.5; Fig. 1.7). Focusing on the
reach-scale as the scale of interest, the logit regression, beta regression, and additive model
explained just 14.1%, 14.4%, and 16.2% of variation in the observed habitat proportions,
respectively, and none of the models had both their slope and intercept terms fall within expected
equivalence regions. Although the intercepts and slopes of the linear fits between observed and
predicted values differed noticeably, none of the three models distinguished themselves from one
another in the other measures of model fit; all three had similar R2 and RMSD values. All three
models performed better in fitting Type II habitat observations (Table 1.5; Fig. 1.7). The
additive model explained the most variation at the reach-scale (37.0%), followed by the logit and
beta regressions at 26.0% and 24.1%, respectively. Although there was little difference among
the model’s RMSD values, equivalence test results suggested that both the logit regression and
additive model outperformed the beta regression; this is supported by the strong deviation of the
beta regression’s slope between observed and predicted values from the expected 1:1 line (Fig.
1.7d). All three models had higher observed model fits at the reach-scale than for individual
transect-scale observations.

40

Table 1.5: Summary of model fit metrics across the different modeling approaches and habitat
types. The root mean squared deviation (RMSD) is the observed average deviation of the
predicted values from the observed. Intercept (a), slope (b), and R2 were calculated from a linear
regression of observed against predicted values. (D) indicates the model was fit using derived
response variables, and * indicates the intercept/slope parameter fell within the expected range
based on equivalence testing.
Model
Logit Regression
Beta Regression
Boosted Additive
Boosted Additive (D)

Transect-Scale Model Fit
Reach-Scale Model Fit
RMSD a
b
RMSD
a
b
Type I Habitat
0.273 0.165
0.469
0.184
0.088
1.103*
0.240 -0.144* 1.725
0.156
-0.123* 1.644
0.257 0.110
0.954*
0.182
0.116
0.826*
0.234 -0.006* 0.984*
0.150
0.010*
0.966*

R2
0.141
0.144
0.162
0.193

Type II Habitat
Logit Regression
Beta Regression
Boosted Additive
Boosted Additive (D)

0.338
0.331
0.315
0.328

0.213*
-0.274*
0.174*
-0.014*

0.641*
1.627
0.704*
0.851*

41

0.251
0.249
0.240
0.257

0.199*
-0.214*
0.206*
0.059*

0.658*
1.450
0.659*
0.751*

0.260
0.241
0.370
0.307

Observed Reach-Scale Habitat (Proportion)

a)

b)

c)

d)

e)

f)

Predicted Reach-Scale Habitat (Proportion)
Figure 1.7: Observed model fit for the Type I (a, c, e) and Type II (b, d, f) habitat values
generated from the logit regression (a, b), beta regression (c, d), and boosted additive model (e,
f). The line represents the linear fit between observed and predicted values.

42

Among all modeling approaches, the additive model fit to the derived quantities of total
larval habitat and habitat ratio had the greatest observed model fit for Type I habitat predictions
and second best fit for Type II habitat (Table 1.5). This model explained 19.3% of variation in
observed, reach-scale Type I habitat, had the lowest observed RMSD among modeling
approaches, and had the nearest fit to the expected 1:1 line between observed and predicted
values (Fig. 1.8). For reach-scale Type II habitat, the model explained 30.7% of variation, and
had the second lowest RMSD value of 0.257. Again, there was a clear difference in the model fit
between the transect-scale and reach-scale observations, with noticeably more variation among
transects than reaches (Fig. 1.8).

43

Observed Reach-Scale Habitat (Proportion)

a)

b)

c)

d)

Predicted Reach-Scale Habitat (Proportion)
Figure 1.8: Observed model fit for Type I (a, b) and Type II (c, b) transect-scale (a, c) and reachscale (b, d) habitat values generated by the boosted generalized additive model fit to the two
derived response variables. The line represents the linear fit between observed and predicted
values.

44

Sampling Success and Model Predictive Ability
All six streams were sampled with varying degrees of success; between 4 and 10 reaches
were sampled within any given stream, with an average sampling intensity of 2.8 access points
and 10.3 transects per reach (Fig. 1.2). Sampling resulted in a total of 47 reaches, 130 access
points, and 483 transects, composing a dataset over 10% of the size of the model training
datasest (11.0%, 13.8%, and 11.6% for reaches, access points, and transects, respectively).
Drought conditions in late July and early August contributed to the difficulty of habitat sampling.
The distribution of sampling effort across the basin and among reaches within each stream
system, however, effectively captured the range of covariate variation within the original model
training dataset (Fig. 1.9).
Similar to the model fitting evaluation results, none of the three models fit to the original
response variables experienced particular success in predicting Type I habitat (Table 1.6; Fig.
1.10); all models were similar in explaining less than one percent of variation in the response.
Predictive performance was noticeably higher and variable among models for Type II habitat,
however. The boosted additive model performed best by explaining 24.6% of variation in the
response, followed by the logit and beta regression at 11.7% and 3.8%, respectively.
The model fit to derived habitat quantities had similarly limited success in predicting
Type I and Type II habitat quantities in the new spatial extents. The additive model explained
less than one percent of variation in Type I habitat and just 18.0% of variation in Type II habitat
(Table 1.6; Fig. 1.10).

45

Covariate Value

Dataset
Figure 1.9: Boxplots of observed covariate variation among reaches sampled in each river
system, where 1 through 6 correspond to the Twin, Peshtigo, Whitefish, Boardman, Pentwater,
and Kalamazoo River, respectively. LM refers to covariate variation in the original model
training dataset.

46

Table 1.6: Model predictive fit for Type I and Type II habitat, among each of the modeling
approaches. (D) indicates the model was fit to derived response variables, and * indicates the
intercept/slope parameter fell within the expected range based on equivalence testing.
Model
Logit Regression
Beta Regression
Boosted Additive
Boosted Additive (D)

Logit Regression
Beta Regression
Boosted Additive
Boosted Additive (D)

RMSD
0.159
0.160
0.145
0.201

Reach-Scale Model Fit
a
b
R2
Type I Habitat
0.139*
0.031
0.000
0.210
-0.280
0.003
0.154*
-0.073
0.001
0.160
-0.056
0.001
Type II Habitat
0.130*
0.500
-0.057* 0.793
0.143*
0.626
0.066
0.503*

0.254
0.289
0.235
0.301

47

0.118
0.040
0.242
0.181

Observed Reach-Scale Habitat (Proportion)

a)

b)

c)

d)

Predicted Reach-Scale Habitat (Proportion)
Figure 1.10: Plots of predictive fit for the logit regression model (a), beta regression model (b),
boosted additive model (c), and boosted additive model fit to the derived habitat quantities (d).
Closed circles and solid lines represent data points and linear fits, respectively, between observed
and predicted values for Type I habitat; open circles and dotted lines represent the same for Type
II habitat.

48

DISCUSSION

Among models fit to original habitat data, the boosted additive model clearly performed
best among the evaluated approaches, followed by the logit regression. The boosted model
explained the most variation in the training dataset for both habitat types, and was the most
accurate in predicting Type II habitat quantities. Despite the still significant amount of variation
left unexplained by the model, the ability of the boosted additive model to account for 24.2% of
variation in the out-of-sample data for Type II dataset supports its potential future utility for
predicting habitat quantities above barriers. The markedly greater performance of this approach
can be attributed to both its application of implicit variable reduction, through the use of iterative
model fitting, and built-in cross-validation techniques to decrease model size, as well as its
nonlinear nature. In fact, the improved ability of the nonlinear model to predict habitat supports
the theorized predominance of nonlinearities and thresholds in fluvial geomorphology and
complex river-landscape systems (Shreve 1979; Phillips 2003; Phillips 2006). Accounting for
the random effects of access point and reach in the logit and beta regressions did not appear to
overcome the restrictions imposed by the assumption of linearity. Furthermore, the use of a beta
distribution in the beta regression, even with zero-one inflation, did not improve model
performance; model predictions appeared to be condensed towards observed habitat averages,
resulting in low model and predictive fits.
The boosted generalized additive model fit to the derived habitat quantities, representing
an attempt to work around possible limitations of the individual Type I and Type II datasets, did
not improve the predictive fit for either habitat type. I primarily selected this modeling approach
to avoid the limited variation, lack of normality, and preponderance of zeroes or similarly small

49

values in modeling Type I habitat. The derived quantities of total larval habitat and ratio
between habitat types were meant to uniquely combine Type I and Type II information, fit new
models, and subsequently back-calculate improved predictions for the two habitat types. Instead
of improving predictive performance, this model resulted in model predictions equivalent or
worse than those for the additive model fit to the original habitat datasets. These results are not
entirely surprising, as the distribution of the habitat ratios and log-odds total habitat quantities
did not readily conform to standard statistical distributions and exhibited a continued dominance
of extreme values.
The apparent mismatch between the assessed model fit and predictive capacity, across
both habitat types and all four modeling approaches, speaks to the difficulty of predicting across
new spatial extents. The intent of this modeling effort, in essence, is simple extrapolation; I built
a model from habitat data collected downstream of dams in order to predict upstream habitat
quantities. The model fit values, especially those for Type II habitat, do compare favorably to
those reported in past studies attempting to predict substrate habitat, with observed R2 ranging
from 0.22 (Mugodo et al. 2006) and 0.32 (Neeson et al. 2012) all the way to 0.50 (Wang et al.
2013). However, the predictive ability for Type I habitat was non-existent and the observed
correlation between observed and predicted Type II habitat quantities declined when moving
from the training to validation datasets. The failure to predict above dams, by itself, might
indicate a complete inability to extrapolate Type I habitat quantity predictions to any new spatial
extents, or simply the inability to predict in river headwaters specifically; regardless, the models
are not useful for predicting Type I habitat in the context of barrier removal scenarios, and only
the boosted additive model appears able to predict Type II habitat in upstream regions.

50

The mismatch between predictive capacity and model fit also reinforces the importance
of out-of-sample data collection for model validation. If the validation results on the out-ofsample data are excluded from these analyses, a different story would have emerged: the logit
and beta regressions appear to perform similarly and all three models appear to explain some
degree of variation in Type I habitat. Instead, the zero/one inflated beta regression performs
exceptionally poorly in predicting both habitat types, and none of the models can predict Type I
habitat quantity above dams reliably, especially compared to the habitat variation explained by
similar studies. This example provides compelling motivation for testing models’ capacity for
spatial extrapolation when evaluating predictive methods against new data (Miller et al. 2004).
The differences in selected important predictors for predicting Type I habitat among the
three modeling approaches reflect the difficulty of predicting stream habitat quantities using only
landscape-scale covariates in a complex, fluvial system, as well as the importance of accounting
for nonlinear relationships in modeling efforts. The two regression models both agreed natural
vegetation land cover, catchment area, fine lithology, and dam effect were statistically important
in influencing fine-substrate habitat quantities

However, the boosted additive model identified

catchment elevation, wetland land cover, and catchment slope as the strongest predictors for
Type I habitat. The apparent lack of consistency between the variables selected by modeling
approaches may reflect a preponderance of purely correlative (i.e., not mechanistic) relationships
and explain the poor predictive capacity upstream of barriers. Given that low velocity regions in
streams, often created by fine-scale channel characteristics, are required for the deposition of the
fine sedimentary materials, exclusively catchment-scale characteristics may be insufficient to
explain trends in this habitat feature (Charlton 2008). Alternatively, simple nonlinearity within
the system may explain the difference between the modeling approaches; natural vegetation land

51

cover and fine lithology, both selected as important by the linear modeling approaches, were the
only two covariates with primarily linear relationships to the response (Fig. 1.5). Rather than
emphasizing real drivers of habitat, linear models may ignore important non-linear effects that
drive habitat variation.
The more consistent identification of important covariates across models for Type II
habitat seems to reflect the increased capacity to predict these habitat quantities and may help
researchers understand dominant influences on stream habitat features. The two linear models
identified catchment area, catchment slope, coarse lithology, precipitation, and groundwater
index as important predictors for Type II habitat, and the boosted additive model used
precipitation, catchment area, and groundwater index most frequently. Furthermore, it is
worthwhile to note that five of the six top predictors for total larval habitat and habitat ratio
response variables, representing related but unique measures of stream habitat, were also
catchment area, precipitation, and groundwater index. Past studies have also identified both
catchment area and catchment slope as strong predictors of stream substrate habitat (Jeffers
1998; Davies et al. 2000; Frappier and Eckert 2007). Indeed, the selection of these landscape
features reinforces our understanding of these characteristics’ influence on stream flow. The
medium- to coarse-grained sands comprising Type II habitat are typically transported as bedload,
and are more frequently transport limited rather than supply limited (Charlton 2008).
Subsequently, stream power, calculated as a function of catchment area and slope, both regulates
the capacity of a river to transport sediment and influences median substrate size (Bagnold 1980;
Whiting et al. 1999). In the case of these analyses, Type II habitat quantities appear to respond
positively to high catchment areas and negatively to high catchment slopes. Furthermore,
groundwater contribution to base flow is known to play a role in regulating extreme flow events

52

in rivers, and particularly in watersheds with limited anthropogenic disturbances like
urbanization or agricultural development (Resh et al. 1988; Allan 2004). It is possible that
increased flow stability in some Lake Michigan systems, indicated by high surficial drainage
rates related to coarse lithology and subsequently high groundwater contributions, promotes the
local accumulation of medium- to coarse-grain sand. Indeed, the direction of groundwater
influence, across all modeling approaches, is positive. Finally, although multiple landscape
attributes likely interact to influence river habitat features, it is worth noting that interactions
among variables were not explicitly considered in my modeling efforts due to a lack of prior
knowledge regarding the form these interactions might take.
This study reflects a unique attempt to predict GLFC-defined, sea lamprey larval habitat
quantities in new spatial extents, using only readily-available, landscape-scale, GIS-derived
covariates and a convenient predictive spatial scale for management. Past attempts to model to
model sea lamprey habitat have either been limited to very fine-scale spatial extents (Neeson et
al. 2007) or focused on evaluating model fit only within the spatial extent used to train the
models (Neeson et al. 2012).

By building the predictive models using only landscape-scale

covariates, the models have the capability to predict anywhere within the Lake Michigan basin,
and perhaps beyond if fundamental landscape-stream relationships happen to hold. Field
sampling-based “ground-truthing” is a newly applied methodology in the prediction of sea
lamprey habitat that assesses this broader modeling utility. Moreover, many previous stream
habitat modeling efforts have focused on strictly linear approaches, which have been shown by
this study’s inclusion of multiple modeling approaches to potentially produce spurious
correlations (Davies et al. 2000; Neeson et al. 2007). Finally, this study’s predictions are made
at the spatial scale of NHDPlus stream reaches; this base scale better facilitates the scaling of

53

predictions to larger extents for the purpose of management (Wang et al. 2013). Sea lamprey
assessment and control efforts are coordinated across Great Lakes basins, and habitat data have
to be accessible on similar scales.
Regrettably, the failure of the assessed modeling approaches to capably predict both Type
I and Type II habitat limits their utility in informing sea lamprey management. Given the
differing production potential among habitat types, the quantity of both habitat types should be
known in order to predict future lamprey production with established modeling approaches, like
the SLaMSE model for the Great Lakes (Jones et al. 2009; Dawson et al. 2016; see chapter 2).
Depending on the desired scale of analysis, predicting sea lamprey production in new rivers
without accompanying larval habitat data or subsequent model improvements may require
locally measured or region-specific average habitat quantities with some realistic measure of
uncertainty. This removes the need to rely on potentially unreliable habitat predictions, and
instead directly incorporates some measure of uncertainty around habitat quantities into model
predictions.
Limitations in the extent, accuracy, and utility of the available response and predictor
variable data can help explain the observed inadequate predictive capabilities, especially for
Type I habitat. First, the observed data forming the basis of modeling efforts may have been
collected with a significant amount of measurement error, typical in the sampling of local-scale
stream habitats (Wang et al. 2013). Second, the available data and selected spatial scale of
model fitting could have led to inaccurate estimates of observed reach-scale habitat quantities, as
over 50% of training dataset’s stream reaches were characterized by four or fewer transects.
Despite findings that a single sampling site on a confluence-to-confluence stream reach can
represent general conditions throughout the reach (Wang et al. 2006), the heterogeneous nature

54

of stream habitat conditions seems to require a greater intensity of sampling effort than four
transects. Furthermore, the lower model fit and much higher degree of scatter for the transectscale observations, relative to the reach-scale observations, for both habitat types reflects the
large amount of variation within reaches and the resulting difficulty in effectively capturing
reach-scale dynamics. Temporal stochasticity and disconnects may have played a significant
role in influencing model and predictive performance. Although utilized predictor data were
ideally summarized within a similar time frame as the initial larval habitat data collection (i.e.,
2004-2008), localized land use disturbances on short-time scales, characteristic in stream
systems with a strong anthropogenic influence, may have added further unexplained variation in
the dataset; these anthropogenic disturbances can temporarily overwhelm natural catchment
characteristics in stream influence (Johnson et al. 1997; Allan 2004). Historical, large-scale land
use changes and climatic conditions, again uncaptured by my chosen variables, may also be
exerting remnant influences on localized stream habitat (Langbein and Schumm 1958). Finally,
the greater temporal mismatch between the predictors and the larval habitat data collected in
2016 may help explain the reduced predictive fit, as I did not account for the potential change in
predictive features like land cover between the collection of the two response datasets.
Despite these limitations, there remains room for improvement in future modeling
endeavors. Modeling systems without well-established mechanistic relationships can result in
purely correlative habitat relationships (Miller et al. 2004), and subsequently produce spurious
predictions. Given the observed low predictive ability of the evaluated models, especially for
Type I habitat, it is certainly possible this occurred in these analyses. This source of error
translates directly into model improvement, as simpler models with proven mechanisticallysupported landscape covariates may have resulted in improved fit. Further empirical work

55

should be done to more closely evaluate driving factors of stream habitat in Michigan streams,
with a focus on untangling linkages among spatial and temporal scales. Exploring additional
variable reduction and nonlinear modeling techniques also has the potential to improve modeling
utility in the future.
These modeling efforts to predict upstream larval habitat quantities represented a unique
effort to generate habitat estimates in Lake Michigan tributaries at a scale useful to management
and inform future sea lamprey control efforts. The results have the potential to inform future
sampling by identifying dominant influences on habitat quantity, but have limited utility in
generating accurate predictions. The insights provided from evaluating model predictive
capacity on out-of-sample data also highlight the importance of “ground-truthing” model
performance. Moving forward, the reliability of future habitat predictions can be improved with
a better understanding of landscape-scale drivers and further work in identifying new variable
reduction or nonlinear modeling techniques.

56

APPENDICES

57

APPENDIX A
WinBUGS SCRIPT FOR THE TYPE 1 HABITAT LOGIT REGRESSION
model{
#Introducing Observed Data and the Random Effect of Access Point (AP)
for (i in 1:Num_Transects) {
T1_Logit[i]~dnorm(mu[i],tau)
mu[i]<-AP_effect[AP [i]]
}
#Modeling the Random Effect of Access Point
for (i in 1:Num_Access_Points){
AP_effect[i]~dnorm(mean_AP[i],tau_AP)
mean_AP[i]<-reach_effect[reachID[i]]
}
#Modeling the Random Effect of Reach
for (i in 1:Num_Reaches){
reach_effect[i]~dnorm(mean_reach[i],tau_reach)
mean_reach[i]<-beta_Precip*Precip[i]+beta_NatVeg*NatVeg[i]+
beta_Agr*Agr[i]+beta_Wetl*Wetl[i]+beta_Imperv*Imperv[i]+
beta_CatchmentArea*Catchment_Area[i]+beta_Slope*Slope[i]+
beta_Elev*Elev[i]+beta_GWIndex*GW_Index[i]+
beta_FineLith*FineLith[i]+beta_SoilPerm*SoilPerm[i]+
beta_StrXing*Str_Xing[i]+alpha_Dam[Dam_Presence[i]]
}
#Prior Specification
sigma~dunif(0,20)
sigma_AP~dunif(0,20)
sigma_reach~dunif(0,20)
beta_Precip~dnorm(0,0.001)
beta_NatVeg~dnorm(0,0.001)
beta_Agr~dnorm(0,0.001)
beta_Wetl~dnorm(0,0.001)
beta_Imperv~dnorm(0,0.001)
beta_CatchmentArea~dnorm(0,0.001)
beta_Slope~dnorm(0,0.001)
beta_Elev~dnorm(0,0.001)
beta_GWIndex~dnorm(0,0.001)
beta_FineLith~dnorm(0,0.001)
beta_SoilPerm~dnorm(0,0.001)
beta_StrXing~dnorm(0,0.001)
for(i in 1:2){
alpha_Dam[i]~dnorm(0,0.001)
#Additional Parameters to Track
tau<-1/(sigma * sigma)

58

tau_AP<-1/(sigma_AP * sigma_AP)
tau_reach<-1/(sigma_reach * sigma_reach)

59

APPENDIX B

fpartial

PARTIAL EFFECTS PLOTS FOR DERIVED RESPONSE VARIABLES

Standardized Covariate Value
Figure 1.11: Partial effects plots, generated from the generalized additive model for total larval
habitat, showing the marginal, nonlinear influence of each continuous covariate on the response.

60

fpartial

Standardized Covariate Value
Figure 1.12: Partial effects plots, generated from the generalized additive model for Type II to
Type I habitat ratios, showing the marginal, nonlinear influence of each continuous covariate on
the response.

61

CHAPTER 2: SIMULATING THE SEA LAMPREY RESPONSE TO DAM REMOVALS IN
THE LAKE MICHIGAN DRAINAGE BASIN

62

INTRODUCTION

Although dams are globally pervasive and can provide many societal benefits, they can
be significant barriers to migratory fish. Indeed, dams have been implicated in the declines of
numerous populations of diadromous species (Limburg and Waldman 2009). Thanks to growing
public preference to increase lotic connectivity and benefit aquatic species, dam removal in the
U.S. is accelerating and numerous large-scale dams have been demolished in systems like the
Penobscot, Carmel, and Elwha Rivers in Maine, California, and Washington, respectively. In
numerous regions of the world, however, stream barriers can actually provide an ecosystem
service by blocking fish migration, specifically in impeding the spread and success of invasive
species (Sharov and Liebhold 1998; VĂŠlez-Espino et al. 2011; Rahel 2013).
Sea lampreys have been major parasites within the Laurentian Great Lakes since their
invasion in the early 20th century, and their arrival and subsequent increase in population
abundance are associated with historic declines in commercially important fish populations,
including lake whitefish (Coregonus clupeaformis) and lake trout (Salvelinus namaycush) (Smith
and Tibbles 1980). The parasitic juvenile stage of this species feeds on Great Lakes fish before
maturing and migrating to Great Lakes tributaries to spawn; the resulting larvae live as
burrowing filter-feeders in these streams for several years before metamorphosing and migrating
back to the lakes to begin their parasitic stage (Applegate 1950). Sea lampreys are currently
controlled to more acceptable population levels in the Great Lakes using a combination of
lampricide application and intentional fragmentation (Smith and Tibbles 1980). A generally
fixed budget is spent annually to treat Great Lakes tributaries, typically selected by expected
treatment effect relative to other tributaries, using lampricides to kill the burrowing larval stages

63

before metamorphosis (Christie et al 2003). Stream barriers play a significant complementary
role by preventing migratory spawners from accessing high quality spawning and larval habitat,
and consequently eliminating the need to treat large sections of rivers for larval sea lampreys
(Hunn and Youngs 1980). Indeed, in addition to pre-existing dams, the Great Lakes Fishery
Commission (GLFC) Sea Lamprey Control Program (SLCP) has actively constructed stream
barriers with the express purpose of blocking sea lamprey migration (Lavis et al. 2003).
Due to the complementary roles of lampricide application and barrier presence, barrier
removals have the potential to drastically reduce the effectiveness of sea lamprey control. In the
Lake Michigan basin alone, dams like the Sixth Street Dam, Union Street Dam, and Calkins
Bridge Dam currently block hundreds of miles of viable spawning and larval habitat in the Grand
River, Boardman River, and Kalamazoo River, respectively, and consequently eliminate the need
for costly control efforts in the form of lampricide applications. If any of these dams were
removed without the construction of a replacement lamprey barrier or an increase in the
lampricide control budget, there would be two options available to control agents: 1) ignore the
newly available habitat and do not allocate any treatment efforts, or 2) allocate lampricide
application efforts to the newly available habitat as needed, at the possible expense of neglecting
other river systems. The first option is unlikely to be considered for large systems like the Grand
River, especially in light of a real-life scenario in Lake Michigan’s Manistique River; an initially
unrecognized barrier failure in this system has been linked to unexpected increases in the lake’s
total sea lamprey population in the late 1990’s (Schleen and Klar 1999; Klar and Young 2004).
Assuming a fixed lake-wide control budget for lampricide applications, the second option
requires a shift in control effort from existing stream systems to the newly available habitat,
potentially decreasing treatment effectiveness across the basin as a whole.

64

Although the general risks of barrier removal in the Great Lakes are accepted by fishery
management agencies, there is a need to better understand the actual magnitude of the sea
lamprey response to barrier removal, as well as implications for the long-term effectiveness of
sea lamprey control. Dam removals are often driven by political and stakeholder considerations
or advocated under the assumption that removals are inherently positive, despite a historical
dearth of studies documenting the effect of removals (Born et al. 1998; Doyle et al. 2003). In
response to the divisive and sometimes arbitrary nature of past dam removal decisions,
researchers are increasingly arguing for more careful, comprehensive consideration of the
potential ecological consequences and an increased role for scientists in providing data on these
consequences (Johnson and Graber 2002; Doyle et al. 2003). Barrier managers are also faced
with many competing objectives and pressures, in addition to ecological ones, that include the
maintenance of public safety and enhancement of recreational opportunities. The development
of formal criteria, supported by the necessary scientific and social information, is one solution
for managing these trade-offs (Pejchar and Warner 2001). In the case of barrier removals in the
Great Lakes, research that equips managers with a more explicit understanding of the effect of
barrier removals on sea lamprey control can help formalize the balancing of trade-offs inherent
in decision-making; this, in turn, can lead to more informed barrier removal decisions and
increasingly positive consequences for society and the environment.
Management strategy evaluation (MSE) modeling, using known information about sea
lamprey life history and control in the Great Lakes, represents a feasible, realistic means to
capture the expected effects of barrier removals on the long-term effectiveness of sea lamprey
control. Management strategy evaluation models are powerful tools for research and
management because they tie together biological, observational, and management (including

65

control) processes, account for sources of uncertainty in each of these processes, and allow
researchers to formally compare the ability of competing management strategies to achieve
specified management objectives (Smith et al. 1999; Harwood and Stokes 2003). Researchers
have already developed a MSE model for sea lamprey, specific to the Great Lakes, incorporating
biological, observational, and management models, each with simulated uncertainty, to assess
the effect of alternative management strategies (Jones et al. 2009). Moreover, each model
component is informed by years or decades of data on life history processes, sampling accuracy,
and control effectiveness. This Sea Lamprey Management Strategy Evaluation (SLaMSE)
model has been used to determine optimal control budgets to achieve target economic injury
levels (Irwin et al. 2012) and to explicitly compare the effectiveness of alternative management
strategies at a basin-wide scale (Dawson et al. 2016). Most importantly, this model can be
modified to simulate barrier removal by creating new artificial habitat units, each representing
newly opened sections of streams, and forecasting the effect of the increased habitat availability
on the equilibrium, lake-wide sea lamprey abundances.
Given the established strengths of MSE models and the breadth of information
incorporated into the existing model for Great Lakes sea lamprey, I used a modified version of
the SLaMSE model to evaluate the effects of barrier removals on the response of the Lake
Michigan sea lamprey population. Lake Michigan was selected as the area of focus for this work
due to the detailed understanding of sea lamprey dynamics in this region (Dawson et al. 2016). I
first assessed the system’s general response to increasing habitat availability through the
systematic, incremental addition of discrete habitat units. This approach allowed me to evaluate
basic trends in the sea lamprey population response under varying habitat perturbations and
characterize the drivers of changing population abundances. I also modeled a specific Lake

66

Michigan barrier removal scenario, using relevant input data and management scenarios
prioritized by sea lamprey control agents, to provide useful information for a contentious,
contemporary barrier removal decision. Both approaches helped explain how a complex, heavily
controlled biological system responds to fundamental changes in habitat availability and
illustrated the ability of an MSE approach to provide applicable predictions of ecological
responses to barrier removals. The success of these modeling efforts further supports the use of
the SLaMSE model in informing future barrier removal scenarios in the Great Lakes.

67

METHODS

Modifying the SLaMSE Model
In modifying the SLaMSE model to accommodate habitat additions and unique
management scenarios, I started with the model used by Dawson et al. (2016). In addition to the
original SLaMSE model structure (Jones et al. 2009), the updated model accounted for historical
trapping effort on 16 Lake Michigan tributaries and assigned streams as regular or irregular
producers, based on expert professional judgment, with different corresponding Ricker
recruitment curves (Dawson and Jones 2009). Regular producers are attributed a more
productive Ricker curve (i.e., a larger Ricker Îą value) than irregular producers.
New model modifications were also included to match recent analyses of historical data.
These modifications included the following: 1) allocating 52% and 48% of all Lake Michigan
spawners to northern and southern tributaries, respectively, prior to assigning spawners to
individual streams based on drainage area and larval abundance, and 2) increasing the influence
of drainage area, relative to larval abundance, in determining spawner allocation to individual
tributaries. These changes were made to match simulated spawner numbers with observed adult
distributions in Lake Michigan rivers that received previous spawner assessments (Mullett et al.
2003; H. Dawson and M.L. Jones, Michigan State University, East Lansing, Michigan,
unpublished analysis).
I made two further modifications to the SLaMSE model to accomplish barrier removal
simulation objectives: 1) I enabled the flexible addition of new habitat units depending on
editable data inputs, and 2) I allowed sea lamprey to infest treatment units that were not treated

68

with lampricides. The second modification allows for estimation of maximum potential parasitic
sea lamprey production, for a given number of allocated spawners, from selected treatment units.

Assessing Population Responses to Systematically Increasing Habitat Availability
I first characterized the general response of the Lake Michigan sea lamprey population to
barrier removals by systematically adding standardized habitat blocks. Each block was assigned
identical attributes, including areas of suitable larval sea lamprey habitat types as defined by the
GLFC (i.e., Type I and Type II; Slade et al. 2003), drainage area, treatment cost, and
miscellaneous larval growth and mortality parameters; these are all attributes of existing
treatment units within the original SLaMSE model. Block attributes were calculated as averages
of all existing treatment units in Lake Michigan (Table 2.1). These habitat additions were
intended to simulate the effect of opening new river systems to sea lamprey (e.g., removing
barriers at the river mouths).

Table 2.1: Base model inputs for a single standardized habitat block in the systematic habitat
additions.
Type I Habitat (m2)
72,872

Type II Habitat
(m2)
313,403

Drainage
Area (km2)
842.8

69

Treatment
Cost ($)
127,864

The systematic addition of treatment units was conducted in two ways: 1) combine new
habitat blocks into an ever larger single treatment unit or 2) add habitat blocks as multiple,
discrete treatment units. These two approaches were intended to contrast the effect of opening a
single large river with the effect of opening numerous small tributaries, with the same overall
increase in total habitat area. When additional habitat blocks were combined to form a single
treatment unit, the total habitat area, drainage area, and treatment cost were correspondingly
increased in a 1:1 relationship; a treatment unit composed of six habitat blocks would therefore
have twice the drainage area, treatment cost, and habitat area as one composed of three such
blocks. I systematically assessed the effect of increased habitat availability by adding three
habitat blocks at a time. This was a convenient scale of analysis because nine additional habitat
units represent a 10% increase in total habitat availability across Lake Michigan. In the end, I
chose to evaluate increasing habitat availability up to an additional 18 habitat units, representing
a 20% increase in habitat availability across the basin.
The influence of two categorical treatment unit attributes, namely producer status and
geographic region, on the sea lamprey response to barrier removals were formally evaluated by
running increasing habitat addition simulations for each possible combination of attributes. New
habitat areas were either assigned the status of regular or irregular producers and the geographic
location of northern or southern Lake Michigan, with the corresponding Ricker recruitment
curves and spawner allocations, respectively. The range of outcomes observed across each of
these attribute scenarios reflected the range of variation expected with future dam removals,
especially as new river systems cannot be readily identified as regular or irregular producers a
priori.

70

For each dam removal scenario, I ran the SLaMSE model for 5000 simulations, with a
100 year time horizon for each simulation; this was intended to capture the full range of
stochasticity in model results. For every simulation, the mean number of total lake-wide adult
spawners across the last ten years was calculated to represent expected equilibrium conditions.
This mean system response was summarized by calculating the percent change in mean
abundance from status quo mean abundance using the equation below, in which the original
value refers to mean status quo abundance unless otherwise stated:
(𝑁𝑒𝑤 𝑉𝑎𝑙𝑢𝑒 − 𝑂𝑟𝑖𝑔𝑖𝑛𝑎𝑙 𝑉𝑎𝑙𝑢𝑒)
× 100
𝑂𝑟𝑖𝑔𝑖𝑛𝑎𝑙 𝑉𝑎𝑙𝑢𝑒
The observed range of variation for each model run represented variability among the
simulation-specific 10-year averages. I also took advantage of the stochastic nature of the
simulations to calculate the proportion of the 5000 simulations, for each habitat addition
scenario, exceeding a high threshold relative to average status quo spawner abundance; I selected
an abundance of 152,266 based on the 90th percentile of simulated lamprey abundances under
status quo conditions. This simulated threshold abundance is similar to the maximum observed
Lake Michigan adult abundance of 141,730 over a recent 10-year period (2005-2014). Finally,
to calibrate the model at the current Lake Michigan control budget of $2.42 million, larval
survival was adjusted until the base model (i.e., no habitat additions) successfully projected the
recent observed average adult abundance of 72,200 (M. Siefkes, Great Lakes Fishery
Commission, Ann Arbor, Michigan, personal communication, 2016).
To explain the observed trends in sea lamprey abundance with increasing habitat
availability, I ran additional simulations to characterize trends in the following model
components: stream-specific parasite production, control budget allocation among the newly
added and original treatment units, and lampricide treatment frequency. Parasitic sea lamprey

71

production reflected the total number of metamorphosed sea lampreys leaving streams in each
year and iteration. Tracking stream-specific parasite production facilitated comparison of the
relative contribution of the new and original treatment units to lake-wide adult abundances.
Additionally, looking at both control budget allocation and treatment frequency helped to explain
why the relative contributions of parasite production from new and original treatment units might
change with increasing habitat availability.
I ran these additional simulations 1000 times over the same 100 year timespan; consistent
with other simulations, only the last ten years of data in each simulation were used to
characterize trends. Simulations were run only for increasing habitat availability in which
regular producing streams were added to northern Lake Michigan, as these attributes produced
the strongest trends in sea lamprey abundance and were therefore more amenable for elucidating
population drivers. I ran these simulations across the full range of increasing habitat availability
and for both the single large and multiple small river additions.

Modeling Barrier Removal on the Grand River
I selected the potential removal of the Sixth Street Dam to demonstrate the utility of
SLaMSE in informing a relevant, potentially high impact barrier removal scenario. The Sixth
Street Dam is located in downtown Grand Rapids, MI, and has served as an incidental lamprey
barrier on the Grand River, Michigan’s longest river system. This barrier currently blocks a
large extent of potential lamprey habitat from access by migrating spawners. Approximately 96
km lies between the Sixth Street Dam and the Webber Dam, the next upstream barrier on the
mainstem, and numerous large tributaries, including the Thornapple, Maple, and Rogue Rivers
drain into the Grand River between the two dams, in addition to many smaller streams (Fig. 2.1).

72

a)

Sixth
Street Webber
Dam Dam
North Lansing
Dam

b)
16

15

1
2

3 4

5
Sixth
Street
Dam

11 1314
6 8 10
7
9 12
17

North Lansing Dam

Figure 2.1: Map of the Grand River mainstem (a) and the modeled Grand River system between
the Sixth Street Dam and North Lansing Dam (b). Numbers in (b) correspond to numbered
treatment unit names in Table 2.2.

73

Recently, there has been pressure by stakeholders to remove this dam, with the primary
goal of recreating the historical rapids and establishing new recreational boating opportunities
(Adair and Sullivan 2015). Thanks in large part to the current relevance and extent of protected
upstream habitat, the Sixth Street Dam removal scenario was listed a high priority for modeling
by SLCP managers (P. Hrodey and M. Siefkes, Great Lakes Fishery Commission, Ann Arbor,
Michigan, personal communication, 2015). Furthermore, this system can also be modeled with
some degree of accuracy given the quantity of compiled data; SLCP surveys for larval habitat
quantities and native lamprey densities were conducted in 2014 and 2015 (Adair and Sullivan
2015), in addition to the recent development of treatment cost estimates for the area.
To simulate the removal of the Sixth Street Dam, I incorporated sixteen new treatment
units between the Sixth Street Dam and Webber Dam, each representing distinct Grand River
tributaries, to the SLaMSE database. The mainstem of the Grand River was deemed unlikely to
host significant larval quantities or require treatment (Fig. 2.1; J. Tews, U.S. Fish and Wildlife
Service, Ludington, MI, personal communication, 2015). Each included treatment unit was
known to contain viable habitat for spawning and larval sea lamprey, and had a uniquely
estimable treatment cost.
Additional attributes of the new treatment units were then estimated using all available
data on the Grand River. The following inputs were estimated for each treatment unit using
available data sources from the SLCP: total stream length, average stream width, the proportion
of Type I and Type II habitat, and treatment cost (Table 2.2). Total stream length was calculated
using the Sea Lamprey Control Map, a decision-support tool developed by the SLCP to help
quantify and visualize the extent of potential sea lamprey habitat blocked by barriers in the Great
Lakes (Great Lakes Fishery Commission 2016). Average stream width and the proportion of

74

Type I and Type II habitat for each treatment unit were estimated by calculating the mean of all
transect data within the system (D. Keffer, U.S. Fish and Wildlife Service, Ludington, MI,
unpublished data). If transect data were absent, the mean width from a similarly-sized, nearby
treatment unit and average habitat proportions across all other treatment units served as
replacements; this was required for six of the sixteen units. Treatment unit length, width, and
habitat proportion data were used in concert to calculate total Type I and Type II habitat area
values. Treatment cost estimates were obtained for all sixteen tributary units from SLCP agents.
Drainage area estimates for each treatment unit, obtained from the National Hydrography
Dataset Plus (1:100,000 scale; http://www.horizon-systems.com/nhdplus), were calculated then
scaled to match the total drainage area of the Grand River. In the case of independent treatment
units within a single watershed, like the Grand River, treatment units are assigned a proportional
fraction of the entire watershed’s drainage area, based on the relative size of their drainage areas
compared to that of all other units. The drainage areas of the new treatment units, in addition to
existing units downstream of the Sixth Street Dam, therefore were appropriately adjusted based
on their individual drainage area estimates. Other inputs, including river-specific features
describing sea lamprey growth and mortality rates, could not be estimated using external
databases and were assigned the same values as existing Grand River treatment units
downstream of the Sixth Street Dam.

75

Table 2.2: Base model inputs for the new treatment units in the Grand River, upstream of the
Sixth Street Dam. Treatment units are ordered top to bottom by increasing upstream distance
from the Sixth Street Dam.
Treatment Unit
Tributary
1. Rogue River
2. Bear Creek
3. Egypt Creek
4. Honey Creek
5. Thornapple River
6. Lowell Creek
7. Flat River
8. Toles Creek
9. Lake Creek
10. Red Creek
11. Timberlin Creek
12. Sessions Creek
13. Bellamy Creek
14. Trestle Creek
15. Prairie Creek
16. Maple River
17. Looking Glass
River*

Length
(m)
4122
1631
7093
5659
479
4172
412
4544
3359
6438
2842
1962
13536
3473
29846
133946
93979

Width
(m)
25.79
6.86
2.73
4.31
16.70
3.21
16.70
2.93
5.25
2.99
2.85
5.25
6.86
5.25
11.80
18.56
20.86

Type I Habitat
(Proportion)
0.0633
0.1675
0.3013
0.2329
0.1675
0.2407
0.1675
0.1324
0.1675
0.1675
0.2277
0.0763
0.1096
0.1675
0.1320
0.1526
0.1195

76

Type II Habitat
(Proportion)
0.1302
0.3503
0.2989
0.3566
0.3503
0.4824
0.3503
0.4532
0.3503
0.3503
0.4010
0.1529
0.3468
0.3503
0.1658
0.3838
0.1306

Treatment
Cost ($)
131125
49230
90440
97110
387840
73170
245870
171800
83870
196700
64485
73385
164700
68085
196700
503161
471201

I also accounted for the potential influence of numerous small tributaries along the Grand
River mainstem, deemed too small or otherwise insignificant to merit consideration as an
independent treatment unit, on future sea lamprey production. Using estimates of the expected
influence of these small tributaries in attracting sea lamprey spawners, based on their observed
drainage areas, I calculated the probable, proportional increase in available habitat in the
upstream Grand River that might result if lampreys used these small tributaries. This increase in
available habitat was then spread equally among the sixteen new treatment units by artificially
increasing stream lengths; a similar adjustment was made for relative drainage area. The total
cost of treating these small tributaries for larval sea lamprey, again estimated by control agents,
was also spread equally among the sixteen explicit treatment units. Across all new treatment
units, these adjustments resulted in a 2.6% inflation in stream area and drainage area and a
20.0% increase in total treatment cost.
Relevant dam removal scenarios and assumptions for the Sixth Street Dam case study
were identified following a meeting with SLCP and Michigan Department of Natural Resources
managers. Three primary management decisions were selected as the focus for modeling work:
the decision to modify the Webber Dam to block sea lamprey, the decision to treat or ignore the
newly available habitat upstream of the Sixth Street Dam, and the decision to maintain or
increase the current lake-wide control budget (Fig. 2.2). Because the Webber Dam currently has
the potential to pass sea lamprey but can be appropriately modified with sufficient justification, I
evaluated the decision to modify the dam by simulating the effect of including the Looking Glass
River; the Looking Glass River is the only major tributary system between the Webber Dam and
the next mainstem barrier. To include this system, I estimated the necessary inputs using the
steps outlined for the sixteen new units downstream of Webber Dam (Table 2.2). The decision

77

to treat or ignore habitat upstream of the Sixth Street Dam is intended to compare the effect of
potentially pulling treatment effort away from other Lake Michigan tributaries with the effect of
allowing uninhibited lamprey production above the Sixth Street Dam, respectively. Finally, if
the system is in fact treated, I both evaluated the effect of treating the system under the current
budget of $2.42 million and estimated the necessary budget increase to prevent a lake-wide
increase in sea lamprey abundance from status quo. I only evaluated the necessary budget
increase for the scenario in which the Webber Dam was modified to block sea lamprey.

Figure 2.2: Visual representation of the modeling scenario combinations considered in
simulating the Sixth Street Dam removal.

78

I also formally assessed the influence of the assumed degree to which sea lampreys
utilize the newly available larval habitat upstream of the Sixth Street Dam. Among all inputs,
larval habitat quantity is especially important to evaluate given its observed role in influencing
recruitment success (Jones et al. 2003) and explicit incorporation into the SLaMSE model (Jones
et al. 2009). I therefore assessed the response of sea lampreys to two levels of assumed habitat
use within added tributaries for each of the control scenario combinations: 10% and 50% habitat
use. The 10% habitat use represents a reasonable approximation of expected lamprey use of total
river length based on professional judgment (A. Jubar, U.S. Fish and Wildlife Service,
Ludington, Michigan, personal communication, 2016) and preliminary analyses indicating that
the observed lengths of existing Grand River treatment units (obtained from the SLCP’s
database) averaged just 10% of the total tributary lengths calculated from the Sea Lamprey
Control Map (A. Jensen, Michigan State University, East Lansing, Michigan, unpublished
analysis). Expected use of total river length is as low as 10% because linear referencing, in
which even marginal lotic habitats unsuitable for larval sea lamprey (e.g., drainage ditches,
ephemeral headwater creeks) are digitized to form stream GIS datasets, can produce
overestimates of total river lengths. I chose to assess the influence of 50% habitat use on the sea
lamprey response in order to evaluate a presumed worst-case scenario for extent of habitat use.
The SLaMSE model was run and summarized in the same manner as for the systematic
habitat additions (i.e., 5000 simulations, 100 year time horizon, ten year averages) for every
scenario and assumption, and the proportions of simulation results above the same status quo
threshold were again calculated.

79

RESULTS

Systematic Habitat Increases
There was a nonlinear increase in the adult sea lamprey population in response to
systematically increasing habitat availability (Figs. 2.3, 2.4). The smallest percent increase in
mean abundance from status quo conditions with a 20% increase in habitat availability was
161%; the greatest increase exceeded 800%. This nonlinear response with increasing habitat
appeared to be adequately characterized by an exponential function across all combinations of
types of habitat addition, producer status, and geographic location (Figs. 2.3, 2.4).
The type of barrier removal (i.e., whether there is one large-scale barrier removal or
multiple small-scale events) influenced the sea lamprey population’s response to barrier removal,
with the addition of a single large stream having the greater effect. The largest percent increase
in abundance for the single stream addition was 885%, in comparison to 452% for multiple
stream additions. This difference in abundance between the types of habitat addition held true
across all combinations of producer status and geographic location. Corresponding with the
different trends in mean abundance, the proportion of simulations with forecasted abundances
greater than the high threshold relative to status quo abundance (152,266) also approached one
more rapidly, relative to the amount of added habitat, when additions were conducted as a single
large river.

80

Geographic Location

Proportion of Simulations > Status Quo 90th Percentile

South

Adult Abundance (Hundreds of Thousands)

Irregular

Producer Status

Regular

North

Habitat Increase (%)
Figure 2.3: Adult sea lamprey abundance trends with increasing habitat availability, assuming
habitat units are added within a single stream unit. Boxes, whisker bars, and open circles
represent the 25th and 75th, 10th and 90th, and 5th and 95th percentiles of simulated adult
abundances, respectively. Solid horizontal lines and black circles represent corresponding
median and mean values, respectively, and the dashed line represents the exponential model fit to
the mean abundance values. The asterisk indicates mean lamprey abundance from Scenario #1
of the Grand River case study, with an assumed 10% habitat use (see Figure 2.7). Gray squares
indicate proportions of simulations with abundances greater than the status quo 90th percentile.

81

Geographic Location

Proportion of Simulations > Status Quo 90th Percentile

South

Adult Abundance (Hundreds of Thousands)

Irregular

Producer Status

Regular

North

Habitat Increase (%)
Figure 2.4: Adult sea lamprey abundance trends with increasing habitat availability, assuming
habitat units are added as independent stream units. Boxes, whisker bars, and open circles
represent the 25th and 75th, 10th and 90th, and 5th and 95th percentiles of simulated adult
abundances, respectively. Solid horizontal lines and black circles represent corresponding
median and mean values, respectively, and the dashed line represents the exponential model fit to
the mean abundance values. The asterisk indicates mean lamprey abundance from Scenario #1
of the Grand River case study, with an assumed 10% habitat use (see Figure 2.7). Gray squares
indicate proportions of simulations with abundances greater than the status quo 90th percentile.

82

The status of opened river systems as regular or irregular producers, as well as their
geographical location, also had implications for the effectiveness of sea lamprey control under
dam removal scenarios. Habitat additions to northern Lake Michigan resulted in higher
abundances than habitat additions in southern Lake Michigan, and new treatment units assigned
as regular producers produced greater overall abundances than habitat additions assigned as
irregular producers (Figs. 2.3, 2.4). Among the two categorical factors of producer status and
geographic location, producer status of the newly added habitat had a slightly greater effect on
resulting adult sea lamprey abundances. With a 20% increase in habitat availability and the
combination of geographic region and type of habitat addition held constant, regular producer
habitat additions resulted in 38.2% to 115% greater mean adult sea lamprey abundances relative
to abundances arising from habitat additions assigned as irregular producers. With the same
20% increase in habitat availability, habitat added to northern Lake Michigan resulted in mean
abundances 23.3% to 92.2% greater than those achieved with habitat added to southern Lake
Michigan.
A combination of novel parasite production from newly added habitat and exponentially
increasing parasite production from the original treatment units, caused in part by a shifting
allocation of treatment effort away from original units to new ones, appeared to drive the
exponential response of adult sea lamprey abundance to habitat increases. There was a trend of
increasing average contribution of new treatment units to basin-wide parasite production with
increasing habitat availability (Fig. 2.5a). This increasing relative contribution from new habitat
constituted novel production of parasites due to barrier removals, but did not completely explain
the magnitude of the sea lamprey response to habitat change, especially if we assume the
production of parasites from the original Lake Michigan treatment units remained unchanged.

83

Perhaps more importantly, increasing habitat availability also caused a steep, concurrent increase
in parasite production within the original treatment units (Fig. 2.5b); the nature of the response
was consistent across both types of habitat addition. This response may be explained in part by
the reduced overall annual treatment frequency among original treatment units with increasing
habitat additions (Fig 2.5c). Concurrently, the average annual allocation of the control budget to
original treatment units declined from $2.42 million to a median of $2.07 and $1.79 million for
the single and multiple treatment unit additions, respectively, when 18 habitat units were added
to the Lake Michigan basin (Fig. 2.6).

84

Relative Contribution of
New Treatment Units to
Total Parasite Production
Total Parasite
Production From
Original Treatment Units
(Hundreds of Thousands)
Annual Treatment
Frequency Across
Original Treatment Units

a)

b)

c)

Habitat Increase (%)
Figure 2.5: Changing model characteristics with increasing habitat. Lines and polygons
represent the median and 10th and 90th percentiles, respectively, across all simulations (a, b) or
treatment units (c). The dotted line and lighter polygon illustrate the effect of adding a single,
large unit, and the solid line and darker polygon illustrate the addition of habitat units as
multiple, discrete treatment units, respectively.

85

Budget Expenditure on Original
Treatment Units ($ Million)

Habitat Increase (%)
Figure 2.6: Average annual budget expenditure on the original treatment units with increasing
habitat availability. The dotted and solid lines illustrate the response when habitat is added as a
single, ever-larger system and multiple, discrete treatment units, respectively.

86

Case Study: Barrier Removal on the Grand River
Treatment unit attributes were successfully estimated for all 17 potential new treatment
units (Table 2.2). The length measurements represent values expected for 10% habitat
utilization; stream length was subsequently adjusted upwards to accommodate the 50% habitat
use assumption. Drainage area values are not listed as they varied depending on the simulation
scenario; the decision to exclude or include the Looking Glass River necessarily altered drainage
area values for all other Grand River treatment units, as they are expressed relative to the Grand
River’s total drainage area.
All management scenarios pertaining to the dam removal forecasted large increases in
adult sea lamprey abundance in Lake Michigan, assuming the control budget remains unchanged
(Fig. 2.7). Among the simulations, the lowest mean percent increase in adult abundance, over
that observed under status quo conditions, was 52%. It occurred when the Webber Dam was
modified to block sea lamprey, new habitat units were treated, and sea lampreys used 10% of
available habitat. For the same scenario, just over 24% of simulations produced abundances
exceeding the status quo 90th percentile. The largest mean increase of 269% occurred when an
unmodified Webber Dam allowed sea lamprey to infest the Looking Glass River, none of the
new habitat units were treated, and sea lampreys used 50% of potentially available habitat.
Approximately 87% of simulations resulted in spawner abundances exceeding the threshold for
this worst case scenario.

87

Adult Abundance (Hundreds of Thousands)

0.87

0.77

0.58
0.44

0.61

0.50
0.34

0.25
0.10

Scenario

Figure 2.7: Expected sea lamprey abundances for each of the management scenarios. Scenarios
#1 and #2 exclude the Looking Glass River, while Scenarios #3 and #4 account for its influence.
New treatment units are treated by the SLCP in Scenarios #1 and #3, and ignored in Scenarios #2
and #4. Boxes, whisker bars, and open circles represent the 25th and 75th, 10th and 90th, and 5th
and 95th percentiles of simulated adult abundances, respectively. The solid horizontal lines and
black circles represent median and mean values, respectively. Numbers above the upper whisker
bars indicate the proportion of simulations greater than the status quo 90th percentile.

88

The decision to modify the Webber Dam, the decision to treat the upstream Grand River,
and the assumed degree of habitat use each had substantial effects on equilibrium sea lamprey
abundances, but the relative magnitude of effects differed. When the decision to treat and
assumed habitat use were otherwise held constant among scenarios, the percent difference in
mean lake-wide sea lamprey abundance between simulations including and excluding the
Looking Glass River ranged between 13.1% and 19.6%, with higher simulated abundances for
scenarios including the Looking Glass River. The decision whether or not to treat the upstream
Grand River system had a larger effect on sea lamprey numbers than the decision to modify
Webber Dam. Compared to simulation scenarios with the new habitat units considered for
treatment, again assuming that all other factors are held constant, the decision to not treat these
units resulted in a 40.4% to 52.1% increase in average adult abundance. Assuming increased
habitat utilization in the new treatment units had a similarly dramatic effect on equilibrium sea
lamprey abundances. With management decisions held constant, the increase in assumed habitat
utilization from 10% to 50% caused the percent differences in abundances between these habitat
use assumptions to range from 34.7% to 49.1%.
For the dam removal scenario in which the Webber Dam is modified to block sea
lamprey, substantial increases in the annual Lake Michigan control budget were needed to
restore mean sea lamprey abundances to levels at or below status quo under the two assumptions
of habitat use. Simulations suggested an annual control budget of $2.62 million per year,
representing a $200,000 increase from the current $2.42 million budget, was needed to decrease
mean abundances to levels at or just below status quo when assumed habitat use was 10% (Fig.
2.8). In comparison, a control budget of $2.78 million was required to accomplish the same

89

objective when assumed habitat use was 50%, representing an annual budget increase of

Adult Abundance
(Hundreds of Thousands)

$360,000.

Current Control Budget ($ Million)
Figure 2.8: Expected sea lamprey abundances when the Sixth Street Dam is removed, the
Webber Dam is modified to block sea lamprey, lampreys are assumed to use 10% (a) or 50% (b)
of maximum potential river length, and the new treatment units are allocated control efforts with
a steadily increasing Lake Michigan control budget. Boxes, whisker bars, and open circles
represent the 25th and 75th, 10th and 90th, and 5th and 95th percentiles of simulated adult
abundances, respectively. The solid horizontal lines and black circles represent median and
mean values, respectively.

90

DISCUSSION

This analysis quantitatively characterized the probable response of the Lake Michigan sea
lamprey population to increasing habitat availability associated with future barrier removals.
Simulations of systematic habitat additions resulted in an exponential population response to
changing habitat, revealed causative drivers of this change within the tightly-controlled
biological system, and illustrated sensitivity to habitat quality attributes. The simulated influence
of alternative management actions and expected cost of maintaining status quo sea lamprey
abundance after removal of the Sixth Street Dam can inform future conservation actions within
the Grand River watershed. More broadly, modeling results revealed complex dynamic behavior
of an intensively managed biological system.
The systematic habitat addition simulations showed that a heavily-controlled invasive
species, like sea lamprey, responds to the localized easing of key management-imposed
constraints in a complex, disproportionate manner. The primary constraints on sea lamprey
population growth in the Great Lakes are habitat limitations created by barriers in large river
systems and lampricide treatment-induced mortality at the larval stage (Christie et al. 2003;
Lavis et al. 2003). When these two constraints were diminished by the addition of habitat and
the subsequent shifting of treatment efforts to these new habitat patches, sea lamprey production
increased in both the new and original river systems and Lake Michigan adult abundances
exhibited a sharp, exponential increase. The population responded to the new habitat through the
immediate contribution of additional parasite production from the new treatment unit(s).
Additionally, the simulated shifts in control effort allocation, following established management
procedures, seemed to drive a more complex response from the original treatment units. The

91

necessity to treat new habitat in a zero-sum control system caused an overall decrease in
treatment frequency across original treatment units, which contributed to increased production
from the these units and helped drive the overall population response. The inter-connected
relaxation of the two critical constraints of habitat and lampricide control determined the sea
lamprey response to increasing habitat availability. Similarly complex responses in population
abundance to changing top-down controls have been observed in outbreaks of mesopredators
(“mesopredator release”), in which small reductions in the abundance of apex predators trigger
disproportion increases in mesopredator abundance (Ritchie and Johnson 2009).
Compounding these shifts in treatment allocation and parasite production is the likely
presence of positive feedback. Positive feedback effects are often observed in complex fisheries
systems under changing predation pressure or habitat conditions (Rose and Cowan Jr 2000;
Walters and Kitchell 2001; Kirby et al. 2009; Audzijonyte et al. 2013). In the Lake Michigan
system, positive feedback refers to increased recruitment and subsequent parasite production
from previously spawner-limited river systems due to increased lake-wide production of
parasites from newly added habitat. Great Lakes sea lampreys seem to exhibit density-dependent
recruitment with a large degree of density-independent variation (Jones et al. 2003; Dawson and
Jones 2009). Great Lakes sea lamprey populations are currently suppressed to abundances less
than 10% of historical abundances by the combination of lampricides and barriers, and this
suppression of adult numbers likely results in spawner-limited tributaries. Increased spawner
allocations to each Lake Michigan tributary, driven by elevated spawner abundances across Lake
Michigan from the addition of habitat, should produce greater future recruitment in these rivers,
and contribute to increasing lamprey production across the basin. Of course, the magnitude of

92

this effect will eventually be limited by compensation (i.e., spawner over-saturation), but this
likely does not occur until after some positive feedback effects are observed.
In total, the combined new production from new habitats, increased production from old
habitats due to shifted control efforts, and positive feedback drove a large response of sea
lamprey abundance from a comparatively small increase in habitat. There is evidence for such a
population response from historical barrier failures in Lake Michigan, as unrestricted
colonization of 220 km of the Manistique River in the late 1990s and early 2000s was associated
with approximately a 100% increase in the estimated Lake Michigan sea lamprey abundance
(Klar and Young 2004). Although the example represents uncontrolled production from a river
system, it suggests that increased habitat availability for a controlled species can have
disproportionate consequences.
The simulated effects of positive feedbacks and shifting treatment efforts in driving the
Lake Michigan sea lamprey population response underscore the importance of explicitly
considering linked impacts when evaluating barrier removal decisions. For species that exhibit
natal homing for spawning, increased spawning and larval habitat in a single stream will
typically translate into increased productivity for that stream. Sea lampreys, however, do not
exhibit natal homing, and the allocation of spawners to rivers in the SLaMSE model reflects this
reality (Bergstedt and Seelye 1995; Jones et al. 2009). The effects of increased production from
a single river for a non-homing species, like sea lamprey, are spread throughout the basin in the
form of broad, but more modest, increases in spawners in multiple streams, with a corresponding
increase in average recruitment success. Furthermore, the need to treat new regions of habitat
made accessible by barrier removals disrupts the finely tuned balance of control effort among the
linked river systems and leads to escalating production from pre-existing treatment units. These

93

results suggest that evaluations of barrier removals focusing on potential fish responses must
consider broader scale, linked implications of increased fish production, especially for systems in
which species do not exhibit strict natal homing and control effort is necessarily balanced among
many streams.
My analysis also revealed that barrier removal decisions need to account for factors in
addition to habitat quantity to accurately assess the effects of barrier removal. Specifically, the
decision to remove a dam from a single large river system versus multiple small rivers, while
providing access to the same overall amount of additional habitat, had a large influence on
lamprey response, as adding ever-larger single river systems produced much greater increases in
lake-wide lamprey abundance than adding multiple discrete units. The difference is due to the
challenge of incorporating increasingly expensive single-system treatments into an annual stream
ranking system; in the current scheme, streams in a given year are treated in descending order
based on the projected cost per expected larva killed until the available budget is exhausted
(Slade et al. 2003). As a system becomes more expensive to treat, the probability that there
would be sufficient budget remaining by the time it ranks for treatment decreases; if there is
insufficient budget remaining at the point a treatment unit is selected, the unit will be passed over
and lower ranked, less expensive systems will instead be treated. There is evidence for this
occurring in the model results: trends in budget expenditure and treatment frequency among
original treatment units appeared to flatten and even rise slightly as increasing habitat availability
in a single large river approaches the final value of 20% (Figs. 2.5c, 2.6). Furthermore,
consistent with the differing treatment frequencies for the types of habitat addition, parasite
production from new habitat for the single river dam removal increased more steeply with
increasing habitat additions than that for multiple dam removals (Fig. 2.5a). Within the existing

94

model framework, an increasingly large single river system is simply more difficult to treat than
multiple, smaller river systems, and ends up producing more parasites than the same amount of
added habitat spread over several smaller rivers.
The geographic location and production potential (i.e., regular or irregular producers) of
affected areas also played important roles in mediating the sea lamprey response to increasing
barrier removals. The geographic location of barrier removals regulates the sea lamprey
response by dictating spawner allocation to new habitat units. Given that a fixed 52% of lakewide spawners are allocated to northern Lake Michigan tributaries, with smaller overall drainage
areas relative to the southern tributaries, habitat additions in northern Lake Michigan will result
in relatively greater spawner allocations to these new habitat units, which in turn increases the
production potential and required treatment frequency in the new habitat units. In comparison,
inherent production potential regulates the response by affecting recruitment dynamics; regular
producers produce more recruits per spawner, on average, than irregular ones. This has the same
effect of increasing production potential and required treatment effort in the new habitat areas.
In a similar vein, the influence of dams and fragmentation on indicator fish species has been
shown to be mediated by river attributes including stream size and thermal regime (Cooper et al.
2016). The variability in projected abundances among these scenarios, dependent on the
attributes of the newly available habitat, therefore should be considered when making barrier
removal decisions.
Furthermore, the high degree of variability across simulations for each of the barrier
removal scenarios reflects very real uncertainty in our understanding of sea lamprey dynamics
and should be explicitly recognized in decision making. One of the strengths of the MSE
approach is the incorporation of multiple sources of uncertainty (Bunnefeld et al. 2011); for the

95

SLaMSE model, these sources included stochasticity in biological processes, larval abundance
assessments, and control efforts. The resulting variability in model results indicates that the
calculated mean responses in abundance are by no means guaranteed outcomes. Instead, the
results indicate expected trends in the sea lamprey response to barrier removals, with a
correspondingly wide range of plausible alternative outcomes. Reporting results as proportions
of simulations with values above some threshold value, selected in this paper as the upper limit
of expected status quo conditions, is one useful means of capturing uncertainty in the simulated
outcomes; decision-makers can then use information like this to assess the risk of key decisions.
In some cases, precautionary approaches are used to deal with uncertainty over management
actions (Kell et al. 1999); in others, uncertainty can be formally accounted for using a decision
analysis process (Irwin et al. 2008). Modeling approaches like that of SLaMSE are not intended
to make decisions (Sharov and Liebhold 1998). Instead, these models provide information on
the best current understanding of the influence of management actions on biological processes,
and the uncertainty around this understanding is a critical component of the simulation results
whenever decisions involve risks.
Simulating the removal of the Sixth Street dam can inform decision making efforts on the
Grand River. The difference between the treated and untreated habitat scenarios points to the
value of treating the upstream Grand River in the case of barrier removal. The decision to ignore
the upstream habitat resulted in higher sea lamprey abundances, despite the dilution of basinwide treatment effort that would have occurred if upstream habitat were to be treated. The
decision to modify the Webber Dam also appeared to be secondary in importance to that
regarding treatment, as the decision to include upstream habitat in treatment efforts had the
relatively bigger effect on the sea lamprey response. Finally, the modest control budget

96

increases required to suppress sea lamprey to status quo abundances, relative to the percent
increases in abundance with no change in budget, seem to justify the recommended increases in
annual control spending. Similar to this application of management strategy evaluation to the
Sixth Street Dam removal, MSE has guided other fisheries-based decisions in the face of
biological and management uncertainty (Kell et al. 1999; Smith et al. 1999; Irwin et al. 2012;
Dawson et al. 2016). Additional SLaMSE-based analyses around the Great Lakes could be used
to inform management decisions for future barrier removal scenarios.
The responsiveness of model results to uncertainty about habitat use by colonizing sea
lampreys underscores the importance of identifying accurate input values. Moving from a 10%
to 50% habitat use assumption, for each of the four evaluated management scenarios,
meaningfully increased expected mean lamprey abundances. Although the 10% assumption can
be considered a reasonable estimate based on professional judgment and preliminary analyses, it
remains a rough approximation, as does the worst case scenario of 50% habitat use. Identifying
reliable habitat area estimates in future modeling endeavors will require both more detailed GIS
data integrating length and width information along streams and improved empirical
understanding of habitat use by spawning sea lamprey within tributary systems. Although not
utilized in these simulation efforts, results from habitat modeling for sea lamprey-specific
attributes like larval habitat quality also have the potential to improve future population
dynamics work (see chapter 1).
An important, implicit assumption of all case study simulations is my assumption that
migrating sea lampreys would utilize all selected upstream tributary systems after the dam
removal. This assumption has been largely supported by dam removal studies on coastal river
systems smaller than the Grand River, in which lampreys were observed to quickly re-colonize

97

previously blocked upstream habitat (Hogg et al. 2013; Lasne et al. 2014). Additionally, sea
lampreys appeared to rapidly colonize upstream reaches of the Manistique River in northern
Michigan, a river section over 220 km in length (Klar and Young 2004). Upstream migration
rates for land-locked sea lamprey have been observed to exceed 3 km/d (Wigley 1959) and
appear to increase with increasing stream size (Kelso and Gardner 2000). There is also evidence
that lampreys can access and utilize upstream tributaries on larger river systems. Tagged and
released anadromous sea lamprey were observed to travel over 18 km to access upstream
tributary systems in Portugal’s River Mondego, a river system draining a watershed slightly less
than half the size of Michigan’s Grand River (Almeida et al. 2002).
I also did not simulate the effect of a seasonal barrier on sea lamprey populations.
Stakeholder groups have proposed the construction of a seasonally-adjusted structure, in place of
the Sixth Street Dam, to operate as a barrier only during sea lamprey migrations (Adair and
Sullivan 2015). I chose not to account for this possibility in simulating the dam removal due to
the uncertainty surrounding its actual installation and potential success at blocking sea lamprey.
If the goals of such a barrier are blocking sea lamprey and allowing passage of other nonjumping, migratory species, the overlapping migration phenologies of Great Lakes fish have
been demonstrated to prevent the successful balancing of such objectives without the installation
of an effective fishway (VĂŠlez-Espino et al. 2011). While certainly not inclusive of all
possibilities, the simulation scenarios and assumptions selected for model analysis reflected best
judgment based on the current understanding of the Grand River system and lamprey migratory
capacities.
As with most modeling approaches, there are key limitations with the SLaMSE model
structure that deserve mentioning. First, the SLaMSE model is concerned with a single species,

98

and does not take any multi-species interactions within the Great Lakes into account. Related to
this, our models assume that sea lampreys in Lake Michigan do not become limited by the
availability of hosts, even at projected abundances exceeding observed historical levels in the
1950’s and 1960’s (Smith and Tibbles 1980). It is reasonable to expect that prey limitation
would occur at some of the higher forecasted abundances and restrict further production (LaRue
1980). The current SLaMSE model also lacks the ability to flexibly treat large treatment units,
with correspondingly high production potential, outside of the annual stream ranking system.
The SLCP sometimes selects high-profile river systems for preferential or consecutive treatments
based on professional judgment, independent of the stream ranking system (Adair and Sullivan
2015). It is possible that a more flexible treatment strategy might reduce the effect of adding
single large habitat units, in both the systematic habitat additions and Grand River case study, by
reducing parasite production from these larger systems. Both modeling limitations have the
potential to forecast future sea lamprey abundances above those which might be observed
following a real-world barrier removal.
Although other modeling-based approaches have been used to inform barrier removal
decisions and predict fish response to changing habitat availability, few, if any, have matched
both the extent and resolution of the reported MSE modeling efforts for sea lamprey populations.
At the broadest extent, barrier removal prioritization efforts, using information like existing
barrier passability, blocked stream length, and estimated cost of removal, are being advocated to
optimize barrier removals across varying spatial extents like the coastal rivers of the Pacific
Northwest (Kemp and O’Hanley 2010). These models help synthesize multiple sources of
information and coordinate decisions at scales similar to the Great Lakes, but often make
simplifying assumptions in relating passability, stream length, and habitat quality to future fish

99

production (Kuby et al. 2005; O’Hanley and Tomberlin 2005; Zheng et al. 2009). At a smaller
spatial extent, landscape models are increasingly used to predict some aspect of fish response to
barrier removal , including projected steelhead (Onchorhynchus mykiss) redd density in the
Pacific Northwest (Steel et al. 2004) and northern pike (Esox lucius) incidence in northern
Sweden (Spens et al. 2007). This type of work also has the potential to inform decision-making
at a broad scale using indirect measures of productivity, but does not directly inform fish
abundance. Finer resolution modeling has occurred to predict fish response to individual barrier
removals using species-specific, population dynamics models for American shad (Alosa
sapidissima), walleye (Sander vitreus), and American eels (Anguilla rostrata) (McCleave 2001;
Cheng et al. 2006; Harris and Hightower 2012). While these models successfully predicted fish
production from systems under different barrier removal scenarios, they did not connect streamspecific production to larger population effects or extrapolate findings to other river systems. To
our knowledge, no previous studies have assessed fish responses to barrier removal using a
detailed, species-specific management strategy evaluation approach across a spatial extent
comparable to Lake Michigan, nor have they explicitly considered the implications of barrier
removals in a coupled management system with trade-offs. Results from these analyses
therefore represent unique, spatially explicit, and comprehensive population predictions for a
suppressed invasive species.
The MSE modeling efforts for the Lake Michigan sea lamprey population represent a
useful quantification of the sea lamprey response to barrier removal and can be used to guide
future barrier decisions. Modeling results revealed new insights into interactions and feedbacks
among the different system components, and highlighted the high degree of uncertainty in the
expected responses to habitat additions. This type of MSE approach has the potential to be

100

replicated in other biological systems similar to that of Great Lakes sea lamprey, in which a
detailed understanding of a controlled species’ population dynamics is already established and
there is a readily manipulated environmental constraint. Future SLaMSE-based analyses of
barrier removals can be improved by implementing some form of prey-based limitation in
population growth, continuing to increase the quality of input data, and clarifying our
understanding of sea lamprey migratory capacity in the Great Lakes.

101

MANAGEMENT RECOMMENDATIONS

The limited success in predicting the quantity of Type I and Type II data in stream
reaches in the Lake Michigan basin, due in part due to sampling methods and inconsistent
sample sizes among reaches, reveals an opportunity to improve the utility of future habitat data
collection for research use. Designing future sampling efforts to capture habitat dynamics at the
scale of natural river features like stream reaches and linking collected data to a common spatial
framework (e.g., National Hydrography Dataset) can facilitate modeling efforts similar to those
reported in chapter 1. Furthermore, changing the sampling protocol to identify more detailed
measures of substrate type (i.e., classify habitat by objective grain size) in place of the current
habitat types can lead to more effective modeling, as Type I and II habitat refer to observed
functional habitat types for sea lamprey rather than sharply defined measures of substrate
composition. Sampling protocols similar to those used in the U.S. EPA’s Wadeable Stream
Assessment program may serve as a useful guide for such changes (USEPA 2004; USEPA
2006). Finally, in the case of sampling new potential habitats for sea lamprey (e.g., river
sections upstream of dams), protocols should incorporate sampling the upstream extent of
spawning habitat and ultimately frame larval habitat quantity relative to this distribution. Given
the inability of larvae to disperse upstream of nest sites, actual potential larval habitat only exists
downstream of the furthest upstream spawning habitat patches.
The most important, and simplest, recommendation from the population simulations is
that an increase in the annual control budget is necessary to avoid a disproportionate sea lamprey
population response with Great Lakes barrier removals. The increases in budget required for
barrier removals like the Grand River represent a small proportional increase in cost relative to

102

the proportional increase in sea lamprey abundance that would otherwise result. Furthermore,
decision makers should account for river attributes other than sheer available habitat when
considering possible effects of barrier removals on sea lamprey populations. The different
spawner allocations and simulated effect of adding habitat in northern and southern Lake
Michigan support the recommendation that expected implications of barrier removals between
these two regions should vary. Additionally, the implications of opening multiple small rivers or
a single large river system for the effectiveness of sea lamprey control need to be carefully
considered on a case-by-case basis. Modeling results suggest addition of a single, expensive
river system can be more difficult to accommodate in the existing ranking system than several
small additions. On the other hand, widely spaced small rivers may be more logistically
challenging to sample and treat, given required travel times between sites.
In addition to the projected costs of barrier removals associated with sea lamprey control,
all relevant objectives of barrier removals should be considered with a similar level of scrutiny.
Sea lamprey control comprises just one aspect of the numerous biological objectives associated
with barrier removal, which range from aquatic ecosystem restoration to the enhancement of
desirable diadromous fish populations. Beyond these biological objectives are equally pressing
social and economic concerns, including public safety, recreation, and power generation.
Researchers are increasingly urging that all of these objectives are formally accounted for and
compared to formal criteria (Pejchar and Warner 2001; Johnson and Graber 2002).
To this end, improving our understanding of the population-level processes influencing
the abundance of desirable fish species represents a major next step to achieve a formal
comparison of biological objectives with dam removals. Simulations of the sea lamprey
population response to increasing habitat were facilitated by decades of persistent research on

103

life history processes for this species. Equivalent levels of understanding regarding population
dynamics processes, or similarly developed simulation-based models, are lacking for desirable
fish species like lake sturgeon (Acipenser fulvescens) and walleye (Sander vitreus).
Incorporating an improved understanding of desirable fish dynamics into comparable models is
critical if we want to explicitly compare the competing objectives of helping certain species and
inhibiting the success of others.
Further research needs include developing an improved mechanistic understanding of
drivers on sea lamprey habitat and incorporating predator-prey dynamics into our understanding
of sea lamprey population processes. Purely modeling-based attempts to predict fine-scale river
habitat rarely explain more than 50% of variation (Wang et al. 2013); in situ observations and
manipulative experiments on fine substrate habitat distribution in the Great Lakes basin can
improve our mechanistic understanding of landscape processes and improve predictive
performance. Finally, time-series analyses that explicitly link the abundances of sea lamprey and
prominent host species should be conducted to inform future population analyses. An improved
understanding of population responses at the simulated extremes will improve the utility of
generated model predictions, particularly in the case of large-scale barrier removals.

104

LITERATURE CITED

105

LITERATURE CITED

Adair, R., and Sullivan, P. 2015. Sea lamprey control in the Great Lakes 2014. Annual Report to
the Great Lakes Fishery Commission.
Allan, J.D. 2004. Landscapes and riverscapes: the influence of land use on stream ecosystems.
Annual Review of Ecology, Evolution, and Systematics 35: 257-284.
Almeida, P.R., Quintella, B.R., and Dias, N.M. 2002. Movement of radio-tagged anadromous
sea lamprey during the spawning migration in the River Mondego (Portugal).
Hydrobiologia 483: 1-8.
Applegate, V. C. 1950. Natural history of the sea lamprey, Petromyzon marinus in Michigan.
U.S. Fish and Wildlife Service Special Scientific Report 55:1-60.
ASDSO (Association of State Dam Safety Officials). 2016. Fact or fiction? [online].
Available from www.damsafety.org [accessed 28 November 2016].
Audzijonyte, A., Kuparinen, A., Gorton, R., and Fulton, E.A. 2013. Ecological consequences of
body size decline in harvested fish species: positive feedback loops in trophic interactions
amplify human impact. Biology Letters 9: 20121103.
Avenetti, L.D., Robinson, A.T., and Cantrell, C.J. 2006. Short-term effectiveness of constructed
barriers at protecting Apache trout. North American Journal of Fisheries Management
26: 213-216.
Bagnold, R.A. 1980. An empirical correlation of bedload transport rates in flumes and natural
rivers. Proceedings of the Royal Society of London 372A: 453-473.
Beamish, F.W.H. 1980. Biology of the North American anadromous sea lamprey, Petromyzon
marinus. Canadian Journal of Fisheries and Aquatic Sciences 37: 1924-1943.
Bergstedt, R.A., and Seelye, J.G. 1995. Evidence for lack of homing by sea lampreys.
Transactions of the American Fisheries Society 124: 235-239.
Born, S.M., Genskow, K.D., Filbert, T.L., Hernandez-Mora, N., Keefer, M.L., and White, K.A.
1998. Socioeconomic and institutional dimensions of dam removals: the Wisconsin
experience. Environmental Management 22: 359-370.
Brenden, T.O., Clark, R.D., Cooper, A.R., Seelbach, P.W., Wang, L., Aichele, S.S., Bissell,
E.G., and Stewart, J.S. 2006. A GIS framework for collecting, managing, and analyzing
multi-scale landscape variables across large regions for river conservation and
management. In Landscape influences on stream habitats and biological assemblages.

106

Edited by R.M. Hughes, L. Wang and P.W. Seelbach. American Fisheries Society,
Bethesda, Maryland.
Brenden, T.O., Wang, L., Clark, R.D., Jr., and Seelbach, P.W. 2007. Comparison between
model-predicted and field-measured stream habitat features for evaluating fish
assemblage-habitat relationships. Transactions of the American Fisheries Society 136:
580-592.
Bunnefeld, N., Hoshino, E., and Milner-Gulland, E.J. 2011. Management strategy evaluation: a
powerful tool for conservation. Trends in Ecology and Evolution 26: 441-447.
Charlton, Ro. 2008. Fundamentals of Fluvial Geomorphology. Routledge, Abingdon, Oxon
ISBN10: 0-415-334453-5.
Cheng, F., Zika, U., Banachowski, K., Gillenwater, D., and Granata, T. 2006. Modelling the
effects of dam removal on migratory walleye (Sander vitreus) early life-history stages.
River Research and Applications 22: 837-851.
Christie, G., Adams, J.V., Steeves, T.B., Slade, J.W., Cuddy, D.W., Fodale, M.F., Young, R.J.,
Kuc, M., and Jones, M.L. 2003. Selecting Great Lakes streams for lampricide treatment
based on larval sea lamprey surveys. Journal of Great Lakes Research 29: 152-160.
Cooper, A.R., Infante, D.M., Wehrly, K.E., Wang, L., and Brenden, T.O. 2016. Identifying
indicators and quantifying large-scale effects of dams on fishes. Ecological Indicators
61: 646-657.
Cowles, M.K., and Carlin, B.P. 1996. Markov Chain Monte Carlo convergence diagnostics: A
comparative review. Journal of the American Statistical Association 91: 883-904.
Curry, R.R. 1972. Rivers – a geomorphic and chemical overview. In River ecology and man.
Edited by R.T. Oglesby, C.A. Carlson and J.A. McCann. Academic Press, New York.
pp. 9-32
D’Angelo, V.S., Muhlfeld, C.C., Miller, B.J., Fredenberg, C.R. 2010. 2010 Summary Report:
Preservation of Native Bull Trout in Glacier National park: Continued Experimental
Suppression of Nonnative Lake Trout in Quartz Lake, Glacier National Park. US
Geological Survey, Rocky Mountain Science Center.
Dawson, H.A., and Jones, M.L. 2009. Factors affecting recruitment dynamics of Great Lakes
sea lamprey (Petromyzon marinus) populations. Journal of Great Lakes Research 35:
353-360.
Dawson, H.A., Jones, M.L., Irwin, B.J., Johnson, N.S., Wagner, M.C., and Szymanski, M.D.
2016. Management strategy evaluation of pheromone-baited trapping techniques to
improve management of invasive sea lamprey. Natural Resource Modeling 29: 448-469.

107

Davies, N.M., Norris, R.M., and Thoms, M.C. 2000. Prediction and assessment of local stream
habitat features using large-scale catchment characteristics. Freshwater Biology 45: 343369.
Derosier, A.L., Jones, M.L., and Scribner, K.T. 2007. Dispersal of sea lamprey larvae during
early life: relevance for recruitment dynamics. Environmental Biology of Fishes 78: 271284.
Doyle, M.W., Harbor, J.M., and Stanley, E.H. 2003. Toward policies and decision-making for
dam removal. Environmental Management 31: 453-465.
Ferrari, S., and Cribari-Neto, F. 2004. Beta regression for modelling rates and proportions.
Journal of Applied Statistics 31: 799-815.
Frappier, B., and Eckert, R.T. 2007. A new index of habitat alteration and a comparison of
approaches to predict stream habitat conditions. Freshwater Biology 52: 2009-2020.
Friedman, J.H. 2001. Greedy function approximation: A gradient boosting machine. The Annals
of Statistics 29: 1189-1232.
Fox, J., and Monette, G. 1992. Generalized collinearity diagnostics. Journal of the American
Statistical Association 87: 178-183.
Gelman, A. 2006a. Multilevel (hierarchical) modeling: What it can and cannot do.
Technometrics 48: 432-435.
Gelman, A. 2006b. Prior distributions for variance parameters in hierarchical models. Bayesian
Analysis 1: 515-534.
Geweke, J. 1992. Evaluating the accuracy of sampling-based approaches to the calculation of
posterior moments. In Bayesian statistics 4. Edited by J. M. Fernando, J.O. Berger, A.P.
David, and A.F.M. Smith. Oxford University Press, Oxford, UK. pp. 169-193.
Great Lakes Fishery Commission (GLFC). 2011. Strategic vision of the Great Lakes Fishery
Commission 2011-2020. Ann Arbor, MI.
Great Lakes Fishery Commission. 2016. Sea Lamprey Control Map [online]. Available from
http://data.glfc.org/ [accessed 30 August 16]
Hansen, G.J. A., and Jones, M.L. 2008. A rapid assessment approach to prioritizing streams for
control of Great Lakes sea lampreys (Petromyzon marinus): a case study in adaptive
management. Canadian Journal of Fisheries and Aquatic Sciences 65: 2471-2484.
Harris, J.E., and Hightower, J.E. 2012. Demographic population model for American shad: will
access to additional habitat upstream of dams increase population sizes? Marine and
Coastal Fisheries: Dynamics, Management, and Ecosystem Science 4: 262-283.

108

Harwood, J., and Stokes, K. 2003. Coping with uncertainty in ecological advice: lessons from
fisheries. TRENDS in Ecology and Evolution 18: 617-622.
Hastie, T., and Tibshirani, R. 1986. Generalized additive models. Statistical Science 1: 297-318.
Hofner, B., Mayr, A., Robinzonov, N., and Schmid, M. 2014. Model-based boosting in R - A
hands-on tutorial using the R package mboost. Computational Statistics 29: 3-35.
Hogg, R., Coghlan, S.M., Jr., and Zydlewski, J. 2013. Anadromous sea lampreys recolonize a
Maine coastal river tributary after dam removal. Transactions of the American Fisheries
Society 142: 1381-1394.
Howell, J.H., Lech, J.J., and Allen, J.L. 1980. Development of sea lamprey (Petromyzon
marinus) larvicides. Canadian Journal of Fisheries and Aquatic Sciences 37: 2103-2107.
Hubbs, C.L., and Pope, T.E.B. 1937. The spread of the sea lamprey through the Great Lakes.
Transactions of the American Fisheries Society 66: 172-176.
Hunn, J.B., and Youngs, W.D. 1980. Role of physical barriers in the control of sea lamprey
(Petromyzon marinus). Canadian Journal of Fisheries and Aquatic Sciences 37: 21182122.
Hynes, H.B.N. 1975. The stream and its valley. Internationale Vereinigung fur Theoretische
und Angewandte Limnologie 19: 1-15.
Irwin, B.J., Wilberg, M.J., Bence, J.R., and Jones, M.L. 2008. Evaluating alternative harvest
policies for yellow perch in southern Lake Michigan. Fisheries Research 94: 267-281.
Irwin, B.J., Liu, W., Bence, J.R., and Jones, M.L. 2012. Defining economic injury levels for sea
lamprey control in the Great Lakes basin. North American Journal of Fisheries
Management 32: 760-771.
Jeffers, J.N.R. 1998. Characterization of river habitats and prediction of habitat features using
ordination techniques. Aquatic Conservation: Marine and Freshwater Ecosystems 8: 529540.
Johnson, L.B., Richards, C., Host, G.E., and Arthur, J.W. 1997. Landscape influences on water
chemistry in Midwestern stream ecosystems. Freshwater Biology 37: 193-208.
Johnson, S.E., and Graber, B.E. 2002. Enlisting the social sciences in decisions about dam
removal. BioScience 52: 731-738.
Jones, M.L., Bergstedt, R.A., Twohey, M.B., Fodale, M.F., Cuddy, D.W., and Slade, J.W.
2003. Compensatory mechanisms in Great Lakes sea lamprey populations: implications
for alternative controls. Journal of Great Lakes Research 29: 113-29.

109

Jones, M.L., Irwin, B.J., Hansen, G.J.A., Dawson, H.A., Treble, A.J., Liu, W., Dai, W., and
Bence, J.R. 2009. An operating model for the integrated pest management of Great Lakes
sea lampreys. The Open Fish Science Journal 2: 59-73.
Kell, L.T., O’Brien, C.M., Smith, M.T., Stokes, T.K., and Rackham, B.D. 1999. An evaluation
of management procedures for implementing a precautionary approach in the ICES
context for North Sea plaice (Pleuronectes platessa L.). ICES Journal of Marine Science
56: 834-845.
Kelso, J.R.M., and Gardner, W.M. 2000. Emigration, upstream movement, and habitat use by
sterile and fertile sea lampreys in three Lake Superior tributaries. North American Journal
of Fisheries Management 20: 144-153.
Kemp, P.S., and O’Hanley, J.R. 2010. Procedures for evaluating and prioritising the removal of
fish passage barriers: a synthesis. Fisheries Management and Ecology 17: 297-322.
KĂŠry, M. 2010. Introduction to WinBUGS for ecologists: A Bayesian approach to regression,
ANOVA, mixed models and related analyses. Academic Press, San Diego, California
Kirby, R.R., Beaugrand, G., and Lindley, J.A. 2009. Synergistic effects of climate and fishing in
a marine ecosystem. Ecosystems 12: 548-561.
Klar, G.T., and Young, R.J. 2004. Integrated management of sea lampreys in the Great Lakes
2004. Annual Report to the Great Lakes Fishery Commission.
Kuby, M.J., Fagan, W.F., ReVelle, C.S., and Graf, W.L. 2005. A multiobjective optimization
model for dam removal: an example of trading off salmon passage with hydropower and
water storage in the Willamette basin. Advances in Water Resources 28: 845-855.
Langbein,W.B., and Schumm, S.A. 1958. Yield of sediment in relation to mean annual
precipitation. Transactions, American Geophysical Union 39: 1076–1084.
LaRue, W. 1980. Lake trout (Salvelinus namaycush) and sea lamprey (Petromyzon marinus)
populations in Lake Michigan, 1971-78. Canadian Journal of Fisheries and Aquatic
Sciences 37: 2047-2051.
The Larval Assessment Work Group of the Sea Lamprey Control Program. 3 April 2013. Larval
assessment sampling protocol using the AbP-2 backpack electrofisher in Great Lakes
streams.
Lasne, E., SabatiĂŠ, M.-R., Jeannot, N., and Cucherousset, J. 2014. The effects of dam removal on
river colonization by sea lamprey Petromyzon marinus. River Research and Applications
31: 904-911.

110

Lavis, D.S., Hallett, A., Koon, E.M., and McAuley, T.C. 2003. History of and advances in
barriers as an alternative method to suppress sea lampreys in the Great Lakes. Journal of
Great Lakes Research 29: 362-372.
Limburg, K.E., and Waldman, J.R. 2009. Dramatic declines in North Atlantic diadromous
fishes. BioScience 59: 955-965.
Liu, F., and Kong, Y. 2015. zoib: An R package for Bayesian inference for beta regression and
zero/one inflated beta regression. The R Journal 7: 34-51.
Lunn, D.J., Thomas, A., Best, N., and Spiegelhalter, D. 2000. WinBUGS -- a Bayesian
modelling framework: concepts, structure, and extensibility. Statistics and Computing 10:
L325-337.
Manion, P.J., and Hanson, L.H. 1980. Spawning behavior and fecundity of lampreys from the
upper three Great Lakes. Canadian Journal of Fisheries and Aquatic Sciences 37: 16351640.
McCleave, J.D. 2001. Simulation of the impact of dams and fishing weirs on reproductive
potential of silver-phase American eels in the Kennebec River basin, Maine. North
American Journal of Fisheries Management 21: 592-605.
Meckley, T.D., Wagner, C.M., and Gurarie, E. 2014. Coastal movements of migrating sea
lamprey (Petromyzon marinus) in response to a partial pheromone added to river water:
implications for management of invasive populations. Canadian Journal of Fisheries and
Aquatic Sciences 71: 533-544.
Michigan Department of Natural Resources. 2016. Dam Management Grants [online].
www.michigan.gov/dnr/ [accessed 28 November 2016]
Midway, S.R., Wagner, T., and Tracy, B.H. 2014. A hierarchical community occurrence model
for North Carolina stream fish. Transactions of the American Fisheries Society 143:
1348-1357.
Midway, S.R., Wagner, T., Tracy, B.H., Hogue, G.M., and Starnes, W.C. 2015. Evaluating
changes in stream fish species richness over a 50-year time-period within a landscape
context. Environmental Biology of Fishes 98: 1295-1309.
Miller, J.R., Turner, M.G., Smithwick, E.A.H., Dent, C.L., and Stanley, E.H. 2004. Spatial
extrapolation: the science of predicting ecological patterns and processes. BioScience
54: 310-320.
Moore, J.D., and Lychwick, T.J. 1980. Changes in mortality of Lake Trout (Salvelinus
namaycush) in relation to increase sea lamprey (Petromyzon marinus) abundance in
Green Bay, 1974-78. Canadian Journal of Fisheries and Aquatic Sciences 37: 2052-2056.

111

Moore, H.H., and Schleen, L.P. 1980. Changes in spawning runs of sea lamprey (Petromyzon
marinus) in selected streams of Lake Superior after chemical control. Canadian Journal
of Fisheries and Aquatic Sciences 37: 1851-1860.
Morman, R.H., Cuddy, D.W., and Rugen, P.C. 1980. Factors influencing the distribution of sea
lamprey (Petromyzon marinus) in the Great Lakes. Canadian Journal of Fisheries and
Aquatic Sciences 37: 1811-1826.
Mugodo, J., Kennard, M., Liston, P., Nichols, S., Linke, S., Norris, R.H., and Lintermans, M.
2006. Local stream habitat variables predicted from catchment scale characteristics are
useful for predicting fish distribution. Hydrobiologia 572: 59-70.
Mullett, K.M. 1997. Agreement of observers classifying larval sea lamprey (Petromyzon
marinus) habitat. M.S. thesis, Northern Michigan University, Marquette, MI.
Mullett, K.M., Heinrich, J.W., Adams, J.V., Young, R.J., Henson, M.P., McDonald, R.B., and
Fodale, M.F. 2003. Estimating lake-wide abundance of spawning-phase sea lampreys
(Petromyzon marinus) in the Great Lakes: extrapolating from sampled streams using
regression models. Journal of Great Lakes Research 29: 240-252.
Multi-Resolution Land Characteristics Consortium. 2015. Multi-Resolution Land Characteristics
Consortium: National Land Cover Database [online]. Available from
http://www.mrlc.gov/index.php [accessed 9 April 2015].
Neeson, T.M., Koonce, J.F., and Whiting, P.J. 2007. Predicting sea lamprey (Petryomyzon
marinus) ammocoete habitat using Geographic Information Systems. Journal of Great
Lakes Research 33: 546-553.
Neeson, T.M., Adlerstein, S.A., and Wiley, M.J. 2012. Towards a process domain-sensitive
substrate habitat model for sea lampreys in Michigan rivers. Transactions of the
American Fisheries Society 141: 313-326.
NRC (National Research Council). 1992. Restoration of aquatic ecosystems: Science,
technology, and public policy. Washington DC, National Academy Press, 552 pp.
O’Hanley, J.R., and Tomberlin, D. 2005. Optimizing the removal of small fish passage barriers.
Environmental Modeling and Assessment 10: 85-98.
Pejchar, L., and Warner, K. 2001. A river might run through it again: criteria for consideration of
dam removal and interim lessons from California. Environmental Management 28: 561575.
Phillips, J.D. 2003. Sources of nonlinearity and complexity in geomorphic systems. Progress in
Physical Geography 27: 1-23.

112

Phillips, J.D. 2006. Evolutionary geomorphology: thresholds and nonlinearity in landform
response to environmental change. Hydrology and Earth System Sciences 10: 731-742.
PiĂąeiro, G., Perelman, S., Guerschman, J.P., and Paruelo, J.M. 2008. How to evaluate models:
observed vs. predicted or predicted vs. observed? Ecological Modelling 216: 316-322.
Pister, E.P. 2010. California golden trout: perspectives on restoration and management.
Fisheries 35: 550-553.
Plummer, M. 2014. JAGS. http://mcmc-jags.sourceforge.net/
Potter, I.C. 1980. Ecology of larval and metamorphosing lampreys. Canadian Journal of
Fisheries and Aquatic Sciences 37: 1641-1657.
R Core Team. 2016. R: a language and environment for statistical computing. R Foundation for
Statistical Computing, Vienna, Austria.
Rahel, F.J. 2013. Intentional fragmentation as a management strategy in aquatic systems.
BioScience 63: 362-372.
Resh, V.H.A., Brown, A.V., Covich, A.P., Gurtz, M.E., Li, H.W., Minshall, G.W., Reice, S.R.,
Sheldon, A.L., Wallace, J.B., and Wissmar, R.C. 1988. The role of disturbance in stream
ecology. Journal of the North American Benthological Society 7: 433-455.
Ritchie, E.G., and Johnson, C.N. 2009. Predator interactions, mesopredator release and
biodiversity conservation. Ecology Letters 12: 982-998.
Robinson, A.P., Duursma, R.A., and Marshall, J.D. 2005. A regression-based equivalence test
for model validation: shifting the burden of proof. Tree Physiology 25: 903-913.
Rose, K.A., and Cowan, J.H., Jr. 2000. Predicting fish population dynamics: compensation and
the importance of site-specific considerations. Environmental Science and Policy 3:
S433-S443.
Schleen, L.P., and Klar, G.T. 1999. Integrated management of sea lampreys in the Great Lakes
1999. Annual Report to the Great Lakes Fishery Commission.
Schmid, M., and Hothorn, T. 2008. Boosting additive models using component-wise P-Splines.
Computational Statistics and Data Analysis 53: 298-311.
Sharov, A.A., and Liebhold, A.M. 1998. Bioeconomics of managing the spread of exotic species
with barrier zones. Ecological Applications 8: 833-845.
Shreve, R.L. 1979. Models for prediction in fluvial geomorphology. Mathematical Geology
11: 165-173.

113

Slade, J.W., Adams, J.V., Christie, G.C., Cuddy, D.W., Fodale, M.F., Heinrich, J.W.,
Quinlan, H.R., Weise, J.G., Weisser, J.W., and Young, R.J. 2003. Techniques and
methods for estimating abundance of larval and metamorphosed sea lampreys in Great
Lakes tributaries, 1995 to 2001. Journal of Great Lakes Research 29: 137-151.
Smith, S.H. 1968. Species succession and fishery exploitation in the Great Lakes. Journal of
Fisheries Research Board of Canada 25: 667-693.
Smith, B.R. 1972. Sea lampreys in the Great Lakes of North America. Ch. V. In Biology of the
lamprey Vol. 1. Edited by M.W. Hardisty and I.C. Potter. Academic Press Inc., London
and New York.
Smith, B.R., and Tibbles, J.J. 1980. Sea lamprey (Petromyzon marinus) in Lakes Huron,
Michigan, and Superior: history of invasion and control, 1936-78. Canadian Journal of
Fisheries and Aquatic Sciences 37: 1780-1801.
Smith, A.D.M., Sainsbury, K.J., and Stevens, R.A. 1999. Implementing effective fisheriesmanagement systems - management strategy evaluation and the Australian partnership
approach. ICES Journal of Marine Science 56: 967-979.
Sorensen, P.W., and Vrieze, L.A. 2003. The chemical ecology and potential application of the
sea lamprey migratory pheromone. Journal of Great Lakes Research 29: 66-84.
Spangler, G.R., Robson, D.S., and Regier, H.A. 1980. Estimates of lamprey-induced mortality
in whitefish, Coregonus clupeaformis. Canadian Journal of Fisheries and Aquatic
Sciences 37: 2146-2150.
Spens, J., Englund, G., and Lundqvist, H. 2007. Network connectivity and dispersal barriers:
using geographical information system (GIS) tools to predict landscape scale distribution
of a key predator (Esox lucius) among lakes. Journal of Applied Ecology 44: 1127-1137.
Steel, E.A., Feist, B.E., Jensen, D.W., Pess, G.R., Sheer, M.B., Brauner, J.B., and Bilby, R.E.
2004. Landscape models to understand steelhead (Onchorhyncus mykiss) distribution and
help prioritize barrier removals in the Willamette basin, Oregon, USA. Canadian Journal
of Fisheries and Aquatic Sciences 61: 999-1011.
Thomas, D.L., Johnson, D., and Griffith, B. 2006. A Bayesian random effects discrete-choice
model for resource selection: Population-level selection inference. The Journal of
Wildlife Management 70: 404-412.
USEPA. 2004. Wadeable streams assessment: field operations manual. EPA 841-B-04-004. U.S.
Environmental Protection Agency, Office of Water and Office of Research and
Development, Washington, DC.

114

USEPA. 2006. Wadeable Streams Assessment: a collaborative survey of the nation’s streams.
EPA 841-B-06-002. U.S. Environmental Protection Agency, Office of Water and Office
of Environmental Information, Washington, D.C.
Wagner, T., Hayes, D.B., and Bremigan, M.T. 2006. Accounting for multilevel data structures
in fisheries data using mixed models. Fisheries 31: 180-187.
Wagner, C.M., Twohey, M.B., and Fine, J.R. 2009. Conspecific cueing in the sea lamprey: do
reproductive migrations consistently follow the most intense larval odour? Animal
Behaviour 78: 593-599.
Walters, C., and Kitchell, J.F. 2001. Cultivation/depensation effects on juvenile survival and
recruitment: implications for the theory of fishing. Canadian Journal of Fisheries and
Aquatic Sciences 58: 39-50.
Wang, L., Seelbach, P.W., and Hughes, R.M. 2006. Introduction to landscape influences on
stream habitats and biological assemblages. In Landscape influences on stream habitats
and biological assemblages. Edited by R.M. Hughes, L. Wang, and P.W. Seelbach
2006. American Fisheries Society, Bethesda, Maryland.
Wang, L., Infante, D., Esselman, P., Cooper, A., Wu, D., Taylor, W., Beard, D., Whelan, G., and
Ostroff, A. 2011. A hierarchical spatial framework and database for the National River
Fish Habitat Condition Assessment. Fisheries 36: 436-449.
Wang L., Brenden, T., Lyons, J., and Infante, D. 2013. Predictability of in-stream physical
habitat for Wisconsin and northern Michigan wadeable streams using GIS-derived
landscape data. Riparian Ecology and Conservation 1: 11-24.
Warton, D.I., and Hui, F.K.C. 2011. The arcsine is asinine: the analysis of proportions in
ecology. Ecology 92: 3-10.
Wells, L. 1980. Lake trout (Salvelinus namaycush) and sea lamprey (Petromyzon marinus)
populations in Lake Michigan, 1971-78. Canadian Journal of Fisheries and Aquatic
Sciences 37: 2047-2051.
Whiting, P.J., Stamm, J.F., Moog, D.B., and Orndorff, R.L. 1999. Sediment-transporting flows
in headwater streams. Geological Society of America Bulletin 111: 450-466.
Wigley, R.L. 1959. Life history of the sea lamprey of Cayuga Lake, New York. U.S. Fish and
Wildlife Service Fishery Bulletin 59: 559-617.
VĂŠlez-Espino, L.A., McLaughlin, R.L., Jones, M.L., and Pratt, T.C. 2011. Demographic
analysis of trade-offs with deliberate fragmentation of streams: control of invasive
species versus protection of native species. Biological Conservation 144: 1068-1080.

115

Youson, J.H., and Potter, I.C. 1979. A description of the stages in the metamorphosis of the
anadromous sea lamprey, Petromyzon marinus L. Canadian Journal of Zoology 57: 18081817.
Zheng, P.Q., Hobbs, B.F., and Koonce, J.F. 2009. Optimizing multiple dam removals under
multiple objectives: linking tributary habitat and the Lake Erie ecosystem. Water
Resources Research 45: W12417.
Zuur, A.F., Ieno, E.N., and Elphick, C.S. 2010. A protocol for data exploration to avoid
common statistical problems. Methods in Ecology and Evolution 1: 3-14.

116