Ecological models for estimating breakpoints and prediction intervals

Abstract The relationships between an environmental variable and an ecological response are usually estimated by models fitted through the conditional mean of the response given environmental stress. For example, nonparametric loess and parametric piecewise linear regression model (PLRM) are often used to represent simple to complex nonlinear relationships. In contrast, piecewise linear quantile regression models (PQRM) fitted across various quantiles of the response can reveal nonlinearities in its range of variation across the explanatory variable. We assess the number and positions of candidate breakpoints using loess and compare the relative efficiencies of PLRM and PQRM to quantitatively determine the breakpoints' location and precision. We propose a nonparametric method to generate bootstrap confidence intervals for breakpoints using PQRM and prediction bands for loess and PQRM. We illustrated the applications using data from two aquatic studies suspected to exhibit multiple environmental breakpoints: relating a fish multimetric index of community health (MMI) to agricultural activity in wetlands' adjacent drainage basins; and relating cyanobacterial biomass to total phosphorus concentration in Canadian lakes. Two statistically significant breakpoints were detected in each dataset, demarcating boundaries of three linear segments, each with markedly different slopes. PQRM generated less biased, more accurate, and narrower confidence intervals for the breakpoints and narrower prediction bands than PLRM, especially for small samples and large error variability. In both applications, the relationship between the response and environmental variables was weak/nonsignificant below the lower threshold, strong through the midrange of the environmental gradient, and weak/nonsignificant beyond the upper threshold. We describe several advantages of PQRM over PLRM in characterizing environmental relationships where the scatter of points represents natural environmental variation rather than measurement error. The proposed methodology will be useful for detecting multiple breakpoints in ecological applications where the limits of variation are as important as the conditional mean of a function.


| INTRODUC TI ON
"Biotic integrity" is defined as the "ability of a habitat to support and maintain a balanced, integrated, adaptive assemblage of organisms having a composition, diversity, and function comparable to that of a natural habitat (Frey, 1977).'' Maintaining biotic integrity of aquatic habitat is one of the primary objectives set forth by the US Clean Water Act of 1972 and by the European Union Water Framework Directive (EU WFD, 2000). Aggressive human activities and rapid urbanization leave only a few areas, if any, that can be considered as "natural habitat." Instead, comparisons are commonly made to reference locations that are subject to a minimum level of anthropogenic stress (Stoddard et al., 2006).
Questions of identifying the degree of disturbance at which biological changes, from reference to nonreference condition, occur across a stress gradient have long been an important consideration (Karr & Chu, 1998;Qian & Miltner, 2015). An ecological threshold ("numerical criterion") is a point of abrupt change of the response variable of an ecological attribute (such as an index of ecological condition) relative to a measure of habitat, such as human-induced disturbance affecting natural habitat (Fahrig, 2001). Such values can serve as guidelines for the protection of environmental condition of sites, and as theoretical conservation or restoration targets (Johnson, 2013;Larned & Schallenberg, 2019).
When the goals are to apply a precautionary principle to environmental management, the measures of biological response should be based on the limits of the relationships of the response given environmental condition (Cade et al., 1999). For example, Karanth et al. (2004) generated prediction bands for tiger densities as a function of their prey densities using standard centrality assumptions for the response variable. In this paper, we describe a method of generating bootstrap prediction bands for a biological response given environmental condition that is applicable to any ecological model without requiring constraining distributional assumptions for the error term.
Methods of identifying environmental disturbance thresholds have been a topic of considerable research. Regression trees (Bunea et al., 1999) and various parametric, nonparametric and Bayesian approaches have been applied to identifying the location of change points in the environmental stress-biological condition relationship (Brenden et al., 2008;Dodds et al., 2010;Qian et al., 2003). The regression tree approach of Bunea et al. (1999), the nonparametric approach of Qian et al. (2003), and the nonparametric deviance reduction method of Brenden et al. (2008) are all similar and based on the concept of classification and regression tree (CART) models.
Bayesian change-point models evaluated by Qian et al. (2003) and Brenden et al. (2008) are the same model. The models reviewed by Brenden et al. (2008) are designed to detect one change point, which is discontinuous at the inflection point, and based on the concept of the nonparametric CART model. Piecewise linear regression model (PLRM) is often used to estimate the location of environmental thresholds-typically a single breakpoint (Ficetola & Denoël, 2009;Shea & Vecchione, 2002;Toms & Lesperance, 2003;Toms & Villard, 2015). Such models may appear to be limited in two perspectives. Firstly, these models often use aggregated community metrics and consider that the regression relationship is based upon the association of two aggregated metrics each drawn from a single population of values. Yet, the sample data are often collected from many different systems representing multiple taxa, each with differing component properties. And the PLRM, which goes through the conditional mean of the metric representing biological response given the other metric representing environmental condition, may not capture the discontinuity in association present in other quantiles of the conditional distribution. This may make it difficult to identify the precise location of the change points (King & Baker, 2011). Secondly, a limiting factor (the least available factor among all factors in the aggregated metrics (Thomson et al., 1996)) may induce an unequal variance pattern in the biological response via interactions among the constituent factors incorporated in the aggregated metric representing environmental condition, and this can alter the relationships near the center of the conditional distribution of the biological response given environmental condition (Cade et al., 1999). Such limiting factors may also cause wedgeshaped relationships in the conditional distribution of the biological response with respect to environmental condition; as a result, the relationships at the edges of the conditional distribution might appear to be more important than the relationship in the center of the distribution (Cade & Noon, 2003;Cade et al., 1999).
Quantile regression has the potential to accommodate these limitations in that it can estimate relationships between variables defined through different quantiles of the conditional distribution of the response variable. As a result, quantile regression models provide a more complete view of the possible relationships between variables than central tendency models (Cade & Noon, 2003). Cade and Noon (2003) gave a general overview of ecological applications of quantile regression, and discussed linear and nonlinear models with both homogeneous and heterogeneous error variances. Using a large simulation study, Cade et al. (2005) showed that the quantile regression model can reveal hidden bias and uncertainty in habitat models. They also showed that the parameters measured at upper ( > 0.5) and lower ( < 0.5) quantiles are less biased than the parameters defined at the mean of the conditional distribution of the response variable given the predictors in the presence of confounding variables. Austin bootstrap, chlorophyll phosphorus, ecological breakpoints, environmental stress, loess, multimetric index, piecewise linear regression, quantile regression growth rates from temperature using a nonlinear quantile regression model. Planque and Buffaz (2008) used a linear quantile regression model to study fish recruitment-environment relationships in marine ecology. Bryce et al. (2010) employed linear quantile regression to predict the maximum decline of vertebrate and macroinvertebrate assemblage responses against streambed sedimentation. Cade et al. (2005) used linear and nonlinear quantile regression models to reveal hidden bias and uncertainty in habitat models in ecology. Brenden et al. (2008) concluded that quantile regression was the most effective means of detecting a single disturbance threshold of the various approaches they investigated.
Nonparametric regression is another frequently used and effective means of studying simple to complex relationships between a biological response and an environmental stress variable. Trexler and Travis (1993) discussed the application of locally estimated scatterplot smoothing (loess) in ecology. Building on the recommendations of Toms and Lesperance (2003), we used loess to subjectively identify the number and positions of candidate ecological breakpoints along an environmental stressor gradient. We then complemented the use of the loess and the piecewise linear regression model (PLRM) approaches with a novel application of a piecewise linear quantile regression model (PQRM) to estimate the location of the environmental thresholds. We propose a method of quantile-based bootstrap confidence interval (CI) for the environmental thresholds using the PQRM and compare estimates with the parametric CI of the breakpoints inferred using the PLRM.
We compared the performances of our methods using simulated data and illustrated the procedures by applying them to two datasets. Our objectives are threefold. The first objective is to display the shape of the relationship between a biological response and an environmental predictor variable. As per the second objective, we identify the locations and precision of the thresholds. Our third objective is to determine the prediction band for the biological response given environmental condition. We then discuss which method (PLRM versus PQRM) provides more precise estimates of environmental thresholds and the prediction bands; and whether PQRM provides more information about the relationships between variables than the PLRM.

| Loess and bootstrap prediction band
Let y and x be the biological response and environmental stress variables, respectively. The nonparametric locally estimated scatterplot smoothing (loess) model (Cleveland, 1979;Cleveland & Devlin, 1988) is defined as where m(x;h) is the smoothed function of interest with smoothing parameter h and ϵ is an independent error term with mean 0 and standard deviation . Loess can capture both linear and nonlinear relationships between variables. Here, the goal of fitting loess is to approximate the number and positions of the breakpoints in a relationship by inspection. We used the R (R Core Team, 2020) statement loess to fit the model (Equation 1).
We generated a bootstrap (Efron & Tibshirani, 1994) prediction band for loess to provide a measure of the variability of the biological response given the environmental condition (Algorithm 1). The algorithm for generating the prediction band is obtained by adapting the methods of Hӓrdle and Bowman (1988) and Davison and Hinkley (1997).
The proposed algorithm assumes that the model is correctly specified and that the residuals are identically and independently distributed. However, the algorithm requires no distributional assumptions for the residuals. Importantly, it can be applied to any method by replacing m(x;h) by the desired model. In this algorithm, steps i-ii capture the sampling variability of the estimated model, and steps iii-iv capture the extra variability due to prediction.
Algorithm 1 Bootstrap resampling method to construct a nonparametric prediction band for the biological response given the predictor.

| Piecewise linear regression model (PLRM)
A PLRM goes through the conditional mean of the response variable and connects two linear segments at each breakpoint. We define to estimate two breakpoints incorporating three linear segments (Seber & Wild, 2003) as where y i and x i are the values for the i th response and predictor variables, respectively, and 1 and 2 are the breakpoints. Here, = ( 0 , 1 , 2 , 3 ) T represents the vector of regression coefficients and = ( 1 , 2 ) T represents the vector of breakpoints. This model assumes that the errors i are iid normal random variable with mean zero and standard deviation . The parameters ( , , ) are estimated using a nonlinear least squares method. We fitted PLRM (Equation 2) using the R package segmented (Muggeo, 2015). To run the nonlinear least squares method, we supplied the initial parameter values estimated from the fitted loess model.

| Piecewise linear quantile regression model (PQRM)
Quantile regression models (Cade & Noon, 2003;Koenker, 2005) are defined through the quantiles of the conditional distribution of the biological response variable. Such models allow one to evaluate relationships among variables through the conditional median of the biological response, as well as the full range of other conditional quantile functions. By supplementing the classical regression model, which is defined at the conditional mean, quantile regression models provide a more complete statistical analysis of the relationships among ecological variables (Mosteller & Tukey, 1977). The PQRM, which is defined at conditional quantiles, provides much richer information in terms of estimating a relationship and breakpoints than the PLRM, which is defined at the conditional mean.
Let m (x; , ) be the th quantile of the conditional distribution of the ecological response given environmental condition as.
Then the PQRM with two breakpoints is defined as: where 1 and 2 are the first and second breakpoints, respectively, defined at the th quantile of the conditional distribution.
Here, = ( 0 , 1 , 2 , 3 ) T represents the vector of regression coefficients and = ( 1 , 2 ) T represents the vector of breakpoints defined at the th quantile. The advantage of quantile regression is that there is no restriction for any distribution of the error term . We used the statement nlrq of the R package quantreg (Koenker et al., 2018) to fit PQRM where the initial values of the parameters are supplied from the fitted loess.
Following Feng et al. (2011), we used wild bootstrap residuals to fit multiple PQRMs defined at the median to calculate the confidence interval (CI) for the breakpoints. The bootstrap CIs for the breakpoints of the PQRM at the median are obtained using Algorithm 2. In this algorithm, f is the kernel density function of the distribution of 3. For b in 1 to B: (i) Generate the weights w i from the two-point mass distribution (ii) Fit a PQRM m(x i ; * , * ) using the bootstrapped observations (x i , y * i ), and calculate the breakpoints (̂

| Relating a wetland fish multimetric index (MMI) to variation in agricultural stress among Laurentian Great Lakes coastal wetlands
The first application relates to estimating threshold effects of a measure of agricultural activity in watersheds draining into the Laurentian Great Lakes on scores of a multimetric index of community composition of fishes in bordering coastal wetlands (Bhagat et al., 2007). Runoff associated with agriculture is a major source of human-induced disturbance affecting natural habitat loss for fishes (Brazner & Beals, 1997;Crosbie & Chow-Fraser, 1999 characterize the risk of degradation of natural habitat using GISbased data. We rescaled the AG (a PCA score) to a 0-1 range with larger numbers reflecting more extensive agricultural activities. The measure of biological condition is a wetland fish multimetric index (MMI), a measure representing the inferred health of the fish assemblage in an ecoregion or watershed. Uzarski et al. (2005) developed and Bhagat et al. (2007) validated the fish multimetric index by assessing fish assemblages in stands of bulrush (Schoenoplectus, spp) in 30 coastal wetland distributed across the US Great Lakes coast (Table A1). Scores vary from 0 to 100, with larger scores representing greater ecological health of the fish assemblage. Traditionally, MMI scores falling in the lowest and highest quintiles are classified as "degraded" and "excellent" conditions, respectively. Bhagat et al. (2007) observed a statistically significant negative linear association between fish IBI and AG scores, but suggested the presence of threshold responses. They did not quantitatively test for the presence of breakpoints.

| Relating cyanobacteria biomass to total phosphorus concentrations among lakes
The second application relates to identifying putative threshold effects of total phosphorus (TP) on the risk of development of harmful algal blooms (dominated by toxigenic Cyanobacteria) in lakes (Beaulieu et al., 2014;Downing et al., 2001;Watson et al., 1992).
TP is a limiting nutrient whose loads to lakes and rivers reflect contributions of sewage from urban centers, agricultural runoff, and other manifestations of human activity (Qian et al., 2003;Reynolds & Walsby, 1975). Cyanobacteria biomass per unit volume (CB) is a standard index of concentration, and often used as a proxy for the risk of toxicity of harmful algal blooms. Cyanobacteria blooms are manifestations of eutrophication whose prevalence is increasing globally (Bullerjahn et al., 2016). CB harbors compounds that can be acutely toxic (Campos & Vasconcelos, 2010;Roegner et al., 2014) and that are linked to diseases such as carcinoma (Labine & Minuk, 2009;Lone et al., 2015). Thus, CB is directly related to risks to human and animal health (Downing et al., 2001;Svendsen et al., 2018).
Opinion on the shape of the relationship between TP and CB is varied. TP is arguably one of the top single predictors of CB (Chlorophyll a), and empirically derived linear models are widely used in lake management (Beaulieu et al., 2014;Dillon & Rigler, 1974Stow & Cha, 2013). However, sigmoidal relationships between TP and CB are also well documented (Chow-Fraser et al., 1994;Downing et al., 2001;Filstrup et al., 2014;Watson et al., 1992Watson et al., , 1997

| Simulation: Evaluating effects of sample size and precision
For the simulation, we generated data from the following model.

| Simulation results
We first present the results in terms of bias, variance, mean-squared error (MSE), coverage, and width for the confidence intervals of the breakpoints identified by the PLRM and PQRM (Table 1). The biases, variances, and mean-squared errors of the point estimates of the breakpoints are smaller for PQRM than for PLRM especially when the sample sizes and degrees of freedoms are small. The differences between the metrics (bias, variance, and MSE) for PLRM and PQRM become smaller as the sample sizes and degrees of freedoms grow larger. However, estimates for the first and second breakpoints are positively and negatively biased, respectively. For the first breakpoint, around which the generated data were more variable, the coverage for PQRM is larger than for PLRM. For the second breakpoint, around which the generated data were less variable, the coverage for PQRM is smaller than the coverage for PLRM. For the first breakpoint, the width of the confidence interval is smaller for PQRM than for PLRM especially for small samples. For the second breakpoint, the width of the confidence interval is smaller for PQRM than for PLRM for both small and large samples.
The mean area within the prediction bands (AWC) and the standard errors against sample size (n) and error degrees of freedom (df) are shown for loess, PLRM, and PQRM (Table 2). For the 80% prediction band, the smallest area is from PQRM followed by loess and PLRM irrespective of whether sample sizes were small or large. For the 95% prediction band, the smallest area is from loess followed by PQRM and PLRM for small sample sizes (n = 30). For the large sample (n = 150), the smallest area is from PLRM followed by PQRM and loess.
In a nutshell, the performance of PQRM is less biased and more accurate than the PLRM in situations with small samples and large error variability reflected by small degrees of freedom.

| Loess and bootstrap prediction band
The relationship between Fish MMI and AG was negative ( Figure 1a

| Piecewise linear regression model
The PLRM identified two breakpoints at AG scores of 0.263 and 0.488, respectively (Table A2)

| Piecewise linear quantile regression model
The PQRM defined at the median ( = 0.50) of the conditional distribution of Fish MMI scores as a function of AG stress estimated the location of two significantly different breakpoints: ̂ 1 = 0.264 (CI of (0.179, 0.347); Table A3) and ̂ 2 = 0.466 (CI of (0.393, 0.553)). with slope ̂ 1 +̂ 2 +̂ 3 = 100.594 and 95% confidence interval (28.830, 193.894). We believe that this counterintuitive result of apparent rising trend in MMI scores in the third segment of the model for AG scores greater than 0.466 is an artifact of the sparse data especially in the vicinity of the estimated breakpoint and small sample size.
Algorithm 1 was used to obtain the prediction bands for the Fish MMI defined at the median (Figure 1c). The 80% and 95% confidence bands covered surface areas of 13.668 and 22.078, respectively. As in the worst-case simulation scenario of small n and df, the prediction bands for PQRM were narrower than those for PLRM.

| Loess and bootstrap prediction band
The fitted loess and the prediction bands indicated that there was a positive relationship between log 10 (CB) and log 10 (TP) (Figure 2a).  The trends of CB against TP estimated by the PLRM (Figure 2b) were consistent in relative magnitude with the patterns subjectively described by the loess. Values for regression coefficients in the first and second segments (but not the third segment) of the regression lines were significantly greater than zero.

| Piecewise linear regression model
The 80% and 95% prediction bands covered surface areas of 3.956 and 6.073 square units, respectively. The PLRM prediction bands were substantially wider than those generated from loess. This is due to large second threshold), the predicted log(CB) was 3.580 with 80% and 95% prediction intervals (2.767, 4.392) and (2.332, 4.827), respectively. The lower limit of log(CB) of 2.767 using the 80% prediction band at log(TP) of 2.015 is higher than the upper limit of log(CB) of 2.443 using the 80% prediction band at log(TP) of 1.012.

| Piecewise linear quantile regression model
The PQRM defined at the median ( = 0.50) of the conditional distribution provided estimates of the two breakpoints ̂ 1 = 1.110 and ̂ 2 = 1.662 (Table A6)

| COMPARISONS
We examined which method provided the narrowest CIs of the breakpoints (Table A8). For the Fish MMI and AG data, PLRM generated a narrower interval for 1 , and PQRM generated a narrower interval for 2 . Similarly, for the CB versus TP data, PLRM generated a narrower interval for 1 , and PQRM generated a narrower interval for 2 .
We also compared the surface areas bounded by the prediction bands (Table A9). For the 80% prediction band of the Fish MMI versus AG data, the surface areas covered by loess, PLRM, and PQRM were 12.089, 17.240, and 13.668 square units, respectively. For the 95% prediction band, the surface areas covered by loess, PLRM and PQRM were 17.730, 27.000, and 22.078 square units, respectively.
The smallest surface area was derived from loess followed by PQRM and PLRM. For the CB versus TP data, the surface areas of 80% prediction bands are 3.445, 3.956, and 3.485 square units, respectively.
The surface areas of 95% prediction bands are 5.948, 6.073, and 6.099 square units, respectively. The smallest surface area was derived from loess. The surface areas from PLRM and PQRM are close to each other.

| Fish multimetric index versus agricultural stress
Loess, PLRM, and PQRM all identified 2 breakpoints in the same locations along the AG stress. Of the two methods from which CIs could be empirically calculated, the 95% CI of the PQRM was 124.44% and 82.47% as wide as those of the PLRM for the lower and upper thresholds, respectively. PQRM produced a narrower CI for the second breakpoint demarcating the sharp transition from second regime to the third.
The PLRM approach identified breakpoints that were significantly different from each other and represented marked discontinuities in the MMI-AG relationship. Fish MMI scores were independent of AG stress range below 0.263, but were a significant negative function of increasing AG from 0.263 to 0.488. Fish MMI was also independent of AG at stress values over 0.488. Similar results were obtained from PQRM. Tests for the slope indicated that fish MMI score was a negative function of AG between the two stress thresholds. The upper limit of fish MMI using the 80% prediction band given AG score of 0.517 (a point above the second threshold) was lower than the lower limit of fish MMI using the 80% prediction band given AG score of 0.210 (a point below the first threshold).
The ecological and environmental management implications (Larned & Schallenberg, 2019) of this interpretation are significant. The results suggest that AG values less than 0.263 have no detectable influence on the fish assemblages, relative to the range of natural variation. In contrast, under high levels of AG (>0.488) management practices that slightly reduce agricultural effects are unlikely to improve fish assemblage condition as expressed in MMI scores. Agricultural changes to watersheds draining into coastal wetlands are only likely to influence fish community condition in wetlands with stress scores between the two breakpoints.
The loess prediction envelope was the smallest: The PQRM enveloped an area of 113.06% and 124.52% for the 80% and 95% prediction bands, respectively, of that produced by the loess. The PQRM produced the narrower prediction bands enveloping an area of 79.28% and 81.77% for the 80% and 95% prediction bands, respectively, than PLRM. These results of smaller prediction band using PQRM than PLRM match the simulation results of the worstcase scenario of small sample size and large error variability.
PQRM revealed important information for the Fish MMI versus AG relationship confirming the presence of discontinuities (thresholds) in the stress-response relationship that were qualitatively suggested by the original researchers (Bhagat et al., 2007). practices that reduce agricultural stresses emanating from watersheds will alter fish community health (Johnson, 2013). If the biological variable is truly controlled by environmental stress, then management actions are likely to be most effective at locations where AG scores fall within the central segment of the stress range. Furthermore, Locations in which MMI scores correspond to lower quantiles ( ≤ 0.60) are likely to be more sensitive to management actions than locations whose fish assemblages have relatively high MMI scores for a particular conditional stress level ( > 0.60).

| Cyanobacterial biomass versus total phosphorus
Loess, PLRM, and PQRM all identified breakpoints in the same loca- Biomass was interpreted to be a monotonically increasing function of TP by Dillon and Rigler (1974) and Beaulieu et al. (2014).
But the identification of 3 significantly different segments of the relationship separated by breakpoints in the nutrient gradient supports the sigmoidal interpretation of the relationship (Filstrup et al., 2014;Watson et al., 1992). The management implications (Larned & Schallenberg, 2019) associated with applying single-slope versus 3-segmented interpretations of the relationship are significant. Use of a linear model to guide management implies that any alteration in TP concentration in a receiving water body can be expected to elicit a cyanobacterial response. In contrast, adherence to a sigmoidal model implies that there are points of inflection beyond which the two variables may behave independently, possibly obviating the need for the control of TP below a particular threshold concentration. Phosphorus has traditionally been regarded as the key nutrient limiting phytoplankton biomass in lakes (Schindler et al., 2008). However, Beaulieu et al. (2014) observed that CB values could be predicted equally well and with similar patterns from concentrations of Total Nitrogen (TN). Yet, a strong correlation between TP and TN, made it difficult for the authors to identify which of the two nutrients was the ultimate predictor of CB. The trisegmented relationship that we observed could reflect colimitation of these two nutrients (Müeller & Mitrovic, 2015). A strength of PQRM is that bivariate relationships can be modeled even when potential confounding factors add variation that reduces the signal to noise ratio in the center of the conditional distribution (Cade & Noon, 2003).
We have proposed methods for quantifying the positions and precision of breakpoints that are subjectively identified by empirical observations, justified by visual analysis of the relationships between a response variable and its predictor using nonparametric loess, which minimized potential observer biases. Loess, PLRM, and PQRM all identified breakpoints in the same stress locations, which were statistically significantly nonoverlapping. As a result, we rule out the possibility that these ecological thresholds are spurious (Daily et al., 2012). However, our observations are consistent with the findings of Daily et al. (2012) that greater precision in estimation is achieved with larger sample sizes and a higher frequency of observations across the environmental stress gradient.
Following Qian (2014), we investigated the goodness-of-fit of our models by inspecting residual plots against environmental stress gradient ( Figure A1). The residuals are distributed symmetrically around the horizontal line at 0 suggesting an absence of bias in the selected models. Again, the larger variability of residuals evident across the center of the environmental gradient than the edges supports the validity of using quantile regression. Cade et al. (1999) showed the applications of linear quantile regression with varying error variances to two ecological applications and pointed out that estimating a range of regression quantiles provides a comprehensive description of biological response patterns for exploratory and inferential analyses. Cade and Noon (2003) explored the applications of both linear and nonlinear quantile regression models and showed how stronger and more useful predictive relationships can be found in other parts of the response distribution than that are observed only in the center. Cade et al. (2005) used linear and nonlinear quantile regression models to habitat data and showed that these models are less biased and uncertain than the classical models defined at the center. Brenden et al. (2008) used a collection of ecological models including the quantile piecewise linear (QPL) with applications in aquatic resource management to model a single breakpoint. Their models were mainly nonparametric in nature and depended heavily on CART. Moreover, the piecewise models were discontinuous in the threshold location between the two linear regimes and did not estimate the precision of the breakpoint. Feng et al. (2011) proposed wild bootstrap for linear quantile regression model via bootstrapping the residuals. In this paper, we have proposed a piecewise linear quantile regression model to detect two thresholds with continuous transition from one linear regime to the adjacent regime. We used wild bootstrap to identify the confidence intervals for the breakpoints and applied the proposed methods to two ecological datasets. The quantile regression estimates for the thresholds are less biased and more accurate than their counterparts of classical piecewise linear regression. Furthermore, the piecewise linear quantile regression model provides the smallest width of the prediction band for the simulated and real data especially for small samples and large error variances.

ACK N OWLED G M ENTS
We thank the Editor, Associate Editor, and two anonymous reviewers for their insightful comments, which substantially strengthened the arguments presented here. We thank Karen Fung for her guidance on an early draft of this manuscript. We also thank Sue Watson for informative discussions on Chlorophyll-phosphorus relationships and Yakuta Bhagat, Donald Uzarski and Lucinda Johnson for discussions about the shape of environmental stress-biological response relationships. We thank two anonymous reviewers for detailed constructive comments on an earlier version of this paper. This research was supported by a grant from the Natural Sciences and Engineering Research Council of Canada to JJHC.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasets used in this paper are available in the Dryad data repository and can be accessed via the following link (https://doi.  Uzarski et al. (2005).