Introduction

Survival is one of the determinants of population viability. Its accurate estimation and manipulation is essential for the management of endangered species, invasive species, plagues and commercially exploited species. Survival is also one of the fundamental drivers of fitness and life-cycle evolution. It has been hypothesized as one of the vital rates through which sympatric kelp species differentiate their adaptation to the environment and consequently, their niche occupation [1]. In isomorphic biphasic life-cycles, dissimilar survival rates between the haploid and diploid generations cause their uneven field abundances [2–6] and are necessary for the prevalence of this type of life-cycle [7]. Furthermore, when these species follow a life strategy dominated by investment in survival, it takes only small differences between haploid and diploid survival for niche differentiation [8–11]. Population density is one of the factors most conspicuously affecting survival. Its most generalized dynamics is decreasing survival with increasing densities due to competition for resources. It is the basic driver of the self-thinning of monospecific stands of plants [12–14], algae [15–18] and animals [19–20]. However, survival (as well as other vital rates) can also decrease at very low densities, a dynamics known as the Allee effects that occurs in animal, plant and algae populations by a multitude of factors [21,22].

The integration of both competition and Allee effects inevitably leads to a non-linear density-dependent survival function. To facilitate parameter estimation, this function can be transposed into a 3^rd degree polynomial in order to density (x). Although it is non-linear, if each of the xⁿ is considered as an independent variable, it becomes a linear-in-the-parameters polynomial of which coefficients can be estimated using Ordinary Least Squares (OLS) or other model I regression algorithms minimizing the vertical residuals. However, this class of methods is only valid when the error in the x estimates is minimal when compared to the error in the y estimates. The appropriate alternative is a Model II regression minimizing the residuals perpendicular to the regression line [23–27]. Unfortunately for the common user (i.e, non-specialized in numerical analysis), this class of methods requires complex calculus, particularly when applied to non-linear models. Furthermore, their numerical stability demands for sample sizes that are unrealistically large for most ecological applications. Hence, model II regressions are not frequently used in ecology and rarely so when non-linear models are involved.

In this study, we derived a density-dependent survival function for Gracilaria chilensis, a red seaweed with a classical isomorphic biphasic life cycle, and tested for differences among its haploid males, haploid females and diploids stages. We developed a non-linear model II regression with an improved line-fit even with small sample sizes. It is an iterative, two-step, Linear-in-the-parameters, Oblique Least Squares (LOLS) method, that is numerically stable in situations where other methods are usually unstable. To compare the parameters obtained from the several algorithms tested, we developed a user-friendly Matlab-based software package that includes a simple and short tutorial (S1 File).

Gracilaria chilensis

This red alga is heavily cultivated and harvested for agar in the intertidal along the Chilean shore. Its life-cycle alternates between free living isomorphic haploid (gametophytes) and diploid (tetrasporophytes) phases. The gametophytes are either male or female [4,7,28]. Diploids (D), haploid males (M) and haploid females (F) were monitored in 5 intertidal rock-pools within 2 sites (Corral 39°52’27”S / 73°24’02” W and Niebla 39°55’47”S / 73°23’57”W) along the margins of the Valdivia river estuary, from October 2009 to February 2011 at 4 month intervals. No specific permissions were required for these locations and activities. The sampling sites were not in protected areas, G. chilensis is not an endangered or protected species, and the sampling method was non-destructive. All individuals within each rock-pool at each time were mapped. M, F and D fronds were identified using first the observation of reproductive structures in a tissue sample under a binocular microscope and the methodology developed in Guillemin et al. [29] for the remaining vegetative individuals. Frond length and diameter were recorded for each individual observed at each census in each rock pool. Individual biomass was represented by the volume (v_i) of a cylinder of equal length and diameter as the frond, significantly correlated with dry weight (r² = 0.877;P < 0.0001; n = 281). Every individual absent after 4 months was re-checked after 8 month for confirmation of death.

The density-dependent survival (s)

The biomass density of each pool was approximated by total frond volume per rock-pool area (V_p = ∑V_i), rescaled to units of L∙m^-2. The survival (s) dependency on V_p followed the shape of a quadratic function (Fig 1), with a peak survival at intermediate densities and increased mortality at lower and higher densities, corresponding to Allee effects and competition, respectively. The quadratic function s = s_max-b₀(V_p-V_opt)² defined the maximum survival (s_max) attained at an optimal biomass density (V_opt). The parameter b₀ represented the sensitivity of s to V_p, with higher b₀ yielding steeper curves. However, the sensitivity of s was observed to be asymmetrical relative to lower vs higher V_p, meaning Allee effects and competition did not act proportionally. Therefore, b₀ was replaced by b₀(1+b₁(V_p-V_opt)), resulting in a density-dependent survival function also including the sensitivity asymmetry parameter b₁ (Eq 1). If survival changed more with higher densities (i.e, Allee effects < competition), then b₁>0. On the other hand, if survival changed more with lower densities (i.e, Allee effects > competition), then b₁<0.

(1)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Linear-in-the-parameters regressions to estimate the survival (s) dependency on population density (V).

A: all stages pooled together for an estimation by (OLS) Ordinary Least Squares, (WLS) Weigthed Least Squares, (IRLS) Iterative Re-weighted Least Squares, (N-RM) maximum likelihood estimation iterated by Newton-Raphson Method (grp) with data grouped, (LOLS) or Linear-in-the-parameters Oblique Least Squares. B: estimation by LOLS independently for (F) females, (M) males, and (D) diploids.

https://doi.org/10.1371/journal.pone.0167418.g001

The procedure to estimate the parameters s_max, V_opt, b₀ and b₁ was to re-formulate Eq (1) as a 3^rd degree polynomial in order to x = V_p i.e, s = β₀+β₁x+β₂x²+β₃x³, and determine the β_i coefficients from a regression algorithm. The several algorithms tested are presented below in increasing order of complexity. Afterwards, the solutions to the survival parameters were obtained from Eq (2a). However, with V_p ranging from 0.1 to 12 and s ranging between 0.2 and 0.8, several methods had unstable calculus leading to unreliable results. Hence, both variables were forced to vary within {-1 1} by transforming x = (2V-^maxV-^minV)/(^maxV-^minV) and y = (2s-^maxs-^mins)/(^maxs-^mins), with the left superscripts ^max and ^min corresponding to the maximum and minimum observed values, respectively. While it solved the numerical instability, the estimation of the model parameters needed Eq (2b) proceeding Eq (2a).

(2a)

(2b)

Ordinary Least Squares (OLS)

Polynomials of orders higher than 1 are non-linear models. However, considering the xⁱ as independent variables, polynomials are linear-in-the-parameters and these can be estimated by OLS. Calculus was simpler transposing the model to matrix algebra (Eq 3). Using matrix notation: y = Xβ+ε, where y was the (n×1) vector with the observations of y, X was the (n×4) matrix with the observations of xⁱ (i = 0,…,3), β was the (4×1) vector with the β_i coefficients and ε was the (n×1) vector with the error in the τ_p observations. Then, β was obtained from β = (X^TX)^-1X^Ty, where the superscript ^T denotes the matrix transpose.

(3)

Iterative Re-weigthed Least Squares (IRLS)

The IRLS estimated β minimizing the weighted y-residuals (Eq 4) while re-estimating their weights (w_p) at each iteration (Eq 5). This procedure demanded for initial estimates to be provided. The interval 0<ϕ<2 guaranteed that the observations with smaller residuals i.e, closer to the central tendency, had iteratively increasing weight on the parameter estimation.

(4)(5)

Maximum Likelihood Estimation (MLE)

A preliminary manual fit suggested the error in the s estimates (ε_p = y_obs-y_est) followed a normal distribution N(0,σ²). Therefore, its associated likelihood function (ℒ) used the normal distribution formula (Eq 6). However, it was preferable to use its log-likelihood function (ℓ = logℒ) because a posterior step required the partial derivatives and these were easier to estimate from a sum (∑) than from a product (∏) and an exponent. The logarithmic function being monotonic, the maximum of ℒ and ℓ coincided. Hence, the estimates of β maximized ℓ(0,σ²|ε_p) from Eqs (7) and (8) knowing that ε_p = f(y_obs, y_est) and y_est = g(β,x).

(6)(7)(8)

When the maximum of ℓ was attained, the first-order partial derivatives of ℓ in order to β (i.e, ∂ℓ/∂β_i) became zero. The Newton-Raphson method [30,31] numerically determined these roots from β_t+1 = β_t-f'(ℓ)/f''(ℓ). This was a system of 4 equations (∂ℓ/∂β₀ = 0, ∂ℓ/∂β₁ = 0, ∂ℓ/∂β₂ = 0 and ∂ℓ/∂β₃ = 0) with 4 unknowns (β₀, β₁, β₂ and β₃). However, σ² was also unknown and needed to be estimated apart from ℓ as it was its forcing function. Then, β assumed the values that maximized ℓ (Eq 8) while minimizing σ² (Eq 9). This duality corresponded to a system of 8 equations, where the four ∂σ²/∂β_i = 0 were added to the four ∂ℓ/∂β_i = 0 while preserving the four β_i unknowns. Finding its roots required adapting the Newton-Raphson method to Eq (10), with the (8×1) vector u of first order partial derivatives (Eq 11a), and the (8×4) Jacobian matrix J extracted from u (Eq 11b).

(9)(10)(11a)(11b)

The computation of u and J was easier with a calculus vectorization similar to that by Vieira et al. [32]. It started by defining M (Eq 12) with entries μ_m,k = ∑x^m∙s_obs^k, · m = 0,…,6 ∧ k = 0,…,2. Matrix algebra estimated faster μ = {x⁰, x¹, x², x³, x⁴, x⁵, x⁶}^T∙{y⁰, y¹, y²}. It was also defined В = {1, -β₀, -β₁, -β₂, -β₃}. This way ℓ could be estimated from Eq (13) and σ² from Eq (14), were n is the sample size and ⊗ is the Hadamard product i.e, matrix element-wise.

(12)(13)(14)

Their partial derivatives had the vectorial representation v_i = ∂В/∂c_i, corresponding to v₀ = {0, -1, 0, 0, 0}, v₁ = {0, 0, -1, 0, 0}, v₂ = {0, 0, 0, -1, 0}and v₃ = {0, 0, 0, 0, -1}. Then, u and J were estimated inserting Eqs (15) to (18) into Eqs (11a and 11b).

(15)(16)(17)(18)

The Newton-Raphson method is very sensitive and prone to overshooting the roots estimates, and so requires initial conditions that are not too far from the solution. Hence, we preformed a manual line-fit to estimate a provisional β by solving the systems of Eqs (2a) and (2b) backwards.

The method was improved by grouping the data into three V classes: (A) 0<V<4, (B) 4<V<5 and (C) 5<V<12 L∙m^-2, with class B giving special attention to the line-fit in the area where s_max and V_opt occurred. The MLE was adapted to run simultaneously over the three classes by concatenating {ℓ_A σ²_A ℓ_B σ²_B ℓ_C σ²_C}^T. The system changed to 24 equations with 4 unknowns, with roots determined from Eqs (10) and (11), where u was a (24×1) vector and J a (24×4) matrix. Eqs (13) to (18) were solved independently for groups A, B and C. An analysis to the line-fit evolution during the Newton-Raphson iterations demonstrated that σ²_A had a far greater weight than σ²_B or σ²_C, with its optimization even driving C into a worsening line-fit. To overcome this bias, the optimization of σ² was standardized to its value at the previous iteration, which only required replacing n by nσ² in the denominators of Eqs (16) and (18).

Linear-in-the-parameters Oblique Least Squares (LOLS)

This is a model II regression minimizing the oblique residuals (i.e, perpendicular to the regression line). However, because the line was curved, each iteration comprised a first step to determine the vectors perpendicular to the third degree polynomial and a second step minimizing the sum of their magnitudes.

Step 1: Each observation was defined by vector {x,y} and their estimated values {ẋ,ẏ}. In cartesian coordinates, each residual was defined as vector A = {x-ẋ,y-ẏ}. A perpendicular to the polynomial was determined by making it perpendicular to the polynomial's derivative, also defined in Cartesian coordinates as vector B = {δ,δ∙∂ẏ/∂ẋ}, by preference using small δ. When vectors A and B were perpendicular, their inner product was null i.e, A◦B = 0. Hence, ẋ and ẏ were determined nesting Eq (19).

(19)

The calculus was performed using matrix algebra. The vectors with the y components of A and B were defined as A_y = {y,-β₀,-β₁ẋ,-β₂ẋ²,-β₃ẋ³} and B_y = {β₁,2β₂ẋ,3β₃ẋ²}. The inner product was re-written as: (20)

The Newton-Raphson method generally took four iterations of ẋ_t+1 = ẋ_t-f(ẋ)/f'(ẋ) to converge into the roots of Eq (20) with a change rate<0,1%. The derivative f'(ẋ) was given by Eq (21), which required the derivatives of A_y (Eq 22) and B_y (Eq 23).

(21)(22)(23)

A provisional β was required, providing the initial conditions for the first time step 1 was implemented. This was arbitrarily chosen aiming to be close to the expected final solution. Subsequent re-runs of step 1 used the β provided by step 2.

Step 2: the magnitude of each oblique residual was given by |res| = ((x-ẋ)²+(y-Ẋβ)²)^1/2, with the (1×4) design matrix Ẋ = {1, ẋ, ẋ², ẋ³}. The minimum sum of their squares i.e, SSE = ∑res², was achieved when ∂SSE/∂β = 0. These roots were found numerically by iterating the Newton-Rahpson method under the form β_t+1 = β_t-SSE'/SSE''. The numerator SSE' was the (4×1) vector with the first order partial derivatives ∂SSE/∂β_i = -∑(y-Ẋβ)ẋⁱ,∀ i = 0,…,3. Upgrading Ẋ to the (n×4) design matrix and y to the (n×1) vector y including all n observations, matrix calculus estimated SSE' = Ẋ^T(Ẋβ-y). The denominator SSE'' was the (4×4) Jacobian matrix with the second order partial derivatives ∂²SSE/∂β_i∂β_j = ∑ẋⁱẋ^j, ∀ i,j = 0,…,3. Matrix calculus estimated it from SSE'' = (Ẋ^TẊ). Merging all yielded the incremental β_t+1 = β_t-(Ẋ^TẊ)^-1Ẋ^T(Ẋβ-y). Only a single iteration within step 1 was performed, after which the perpendicular coordinates of the residuals needed to be re-estimated.

Results and Discussion

Gracilaria chilensis exhibited a clear pattern of maximum survival attained at optimal intermediate densities (Fig 1). The decreased survival of the ramets at higher densities has been widely detected in algae and attributed to competition [15–18]. But their decreased survival at lower densities (Allee effects) has hardly been reported. Crowding has been demonstrated to protect from desiccation in intertidal algae stands [33] and we hypothesize that it may also protect against hydrodynamic stress. The LOLS applied to each of the life-cycle stages showed no significant differences between males, females and diploids for any of the survival function parameters (Fig 1B and Table 1). If there were any actual differences between the density-dependent survival of the stages, these differences were masked by the error. Whenever testing the relation between variables x and y, there are two usual sources of error: (i) measurement error, and (ii) the influence of extra variables over x and/or y. In our case, a weakness of the data set constituted a third source of error undermining the parameter estimation: our observations predominantly corresponded to low population densities and thus, the line-fits that described the dynamics at high population densities were relatively weak. Our four monthly survey interval enhanced this issue. In hind side, a shorter time interval would have provide better density-dependent survival data and a better definition of the full function.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. The ploidy specific survival functions with parameters estimated by the LOLS (value) and the significance (p) of their differences among (M) males, (F) females and (D) diploids, as estimated from randomization tests.

Sample size n = 35.

https://doi.org/10.1371/journal.pone.0167418.t001

Finding the best parameter estimates also required solving problems specific to the numerical methods (Fig 1A and Table 2). The OLS line-fit concentrated the positive residuals on the centre, largely underestimating the maximum survival rate (s_max) at optimal densities. The WLS represented only a slight improvement, while the IRLS represent no improvement at all. Besides, when diploids, males and females were pooled separately (n = 3 months×3 years×5 pools = 35) the results varied significantly with the ϕ and the initial estimates provided. Upon smaller sample sizes, particular observations erroneously became anchor-points (or attractors) to the IRLS line-fit. The standard MLE provided a line-fit equal to the OLS. When the data was split into the three V classes, it became evident that the optimization of the fit at lower densities was worsening the fit at intermediate and higher densities. While σ²_A decreased from 335∙10⁻⁴ to 187∙10⁻⁴, σ²_B increased from 73∙10⁻⁴ to 155∙10⁻⁴ and σ²_C increased from 117∙10⁻⁴ to 123∙10⁻⁴. The MLE with these three V classes weighted by the inverse of their variance provided a line-fit that was only slightly different. All these methods fall upon the class of model I regressions, characterized by a minimization of the vertical residuals i.e, of the error in the response variable. Besides method-specific fails, their generalized bias also resulted from the fact that the error in the estimation of the predictor was larger than the error in the estimation of the response. In fact, we had more confidence in the estimation of the survival rates (ramets were either dead or alive) than in the estimation of the individual sizes and consequent population biomass densities. In such cases, model I regressions are inadequate and model II regression are the proper alternative [25–27]. Classical error-in-the-variables algorithms for model II regressions estimate the true parameter values by use of high-order moments (or cumulants). However, besides extremely complicated for non-linear models, the use of high-order moments turns them numerically unstable, requiring extremely large data sets (n>900) for reliable estimates [34–36]. The instrumental variables (an alternative class of error-in-the-variables algorithms) demands for a subset of data with known errors to train the algorithm into the estimation of the true parameter values. Principal Components Analysis (PCA) and Reduced Major Axis (RMA) are model II regression algorithms with a different strategy: they do not aim at estimating the true parameter values but simply at minimizing the oblique residuals i.e., perpendicular to the line-fit [23–27]. Although this simplifies and numerically stabilizes them, working fine with smaller data sets, they are exclusive to linear models, and thus inapplicable to our linear-in-the-parameters third degree polynomial.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. The overall survival function with parameters estimated by different methods.

Sample size n = 105.

https://doi.org/10.1371/journal.pone.0167418.t002

It was based in the PCA and RMA rationale that we developed the Linear-in-the-parameters Oblique Least Squares (LOLS) regression method, minimizing the squared residuals perpendicular to a curve line. This procedure distributed the positive and negative residuals evenly along the x axis, contrasting with the other tested algorithms (Fig 1A). However, as the LOLS had no closed form solution, it required being iterated over two steps: one to update the orientations of the vectors of residuals and another to minimize their magnitudes. This method was still fallible when (-if) working with the original variables. Whenever at least one of the variables ranged beyond |1| or their variances were conspicuously unbalanced, the method became numerically unstable. Often, the vertical (or the horizontal) component of the oblique residuals became meaningless and the LOLS was in practice performing a horizontal (or vertical) regression. This problem was solved by transforming both variables to range within -1<x_j<1. Overall, although computationally more demanding, the LOLS was still fast providing a solution that was much more satisfactory than the solutions provided by the former algorithms (Fig 1A and Table 2). It estimated conspicuously higher maximum survival rate (s_max) at lower optimum density (V_opt) with higher density-dependence (b₀) and bigger asymmetry between competition and Allee effects (b₁). An apparent draw-back was the function increasing back again at the extreme of higher densities. However, this error was only a consequence of the small number of observations in that area of the function, it should not happen with more complete data sets and hence should be negligible for the majority of future model applications. The consistency of the LOLS fits among stages (n = 35) contrasted with the inconsistency demonstrated by other methods, and in particular by the IRLS.

Conclusions

The accurate estimation of survival is essential for ecological and evolutionary studies, as well as for population management. As found in other organisms, the survival of Gracilaria chilensis is sensitive to both competition and Allee effects, and thus is optimized at intermediate densities. Although computationally more demanding, the LOLS algorithm was the method providing the best description of this non-linear function. This method is better suited in situations demanding for non-linear model II regression methods i.e, when the error in the estimation of the predictor x is comparable to the error in the estimation of the response y. Furthermore, its conjugation with the transformations of both predictor and response to vary within {-1, 1} made the LOLS robust in situations where the other methods tested were fallible and numerically unstable. Therefore, the LOLS should have a wide applicability.

Supporting Information