Linear-In-The-Parameters Oblique Least Squares (LOLS) Provides More Accurate Estimates of Density-Dependent Survival

Survival is a fundamental demographic component and the importance of its accurate estimation goes beyond the traditional estimation of life expectancy. The evolutionary stability of isomorphic biphasic life-cycles and the occurrence of its different ploidy phases at uneven abundances are hypothesized to be driven by differences in survival rates between haploids and diploids. We monitored Gracilaria chilensis, a commercially exploited red alga with an isomorphic biphasic life-cycle, having found density-dependent survival with competition and Allee effects. While estimating the linear-in-the-parameters survival function, all model I regression methods (i.e, vertical least squares) provided biased line-fits rendering them inappropriate for studies about ecology, evolution or population management. Hence, we developed an iterative two-step non-linear model II regression (i.e, oblique least squares), which provided improved line-fits and estimates of survival function parameters, while robust to the data aspects that usually turn the regression methods numerically unstable.


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][3][4][5][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][9][10][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][13][14], algae [15][16][17][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 n 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][24][25][26][27]. Unfortunately for the common user (i.e, nonspecialized 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 2 = 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 0 (V p -V opt ) 2 defined the maximum survival (s max ) attained at an optimal biomass density (V opt ). The parameter b 0 represented the sensitivity of s to V p , with higher b 0 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 0 was replaced by b 0 (1+b 1 (V p -V opt )), resulting in a density-dependent survival function also including the sensitivity asymmetry parameter b 1 (Eq 1). If survival changed more with higher densities (i.e, Allee effects < competition), then b 1 >0. On the other hand, if survival changed more with lower densities (i.e, Allee effects > competition), then b 1 <0.
The procedure to estimate the parameters s max , V opt , b 0 and b 1 was to re-formulate Eq (1) as a 3 rd degree polynomial in order to x = V p i.e, s = β 0 +β 1 x+β 2 x 2 +β 3 x 3 , 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-max V-min V)/( max V-min V) and y = (2s-max s-min s)/( max s-min s), 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).

Ordinary Least Squares (OLS)
Polynomials of orders higher than 1 are non-linear models. However, considering the x i 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 (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 T X) -1 X T y, where the superscript T denotes the matrix transpose.   Weigthed Least Squares (WLS) The WLS estimated β weighting the observations by the (p×p) matrix W, the inverse of their variance-covariance matrix. Hence, the regression model was upgraded from the previous OLS to become β = (X T WX) -1 X T Wy.

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.

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,σ 2 ). Therefore, its associated likelihood function (L) used the normal distribution formula (Eq 6). However, it was preferable to use its log-likelihood function (ℓ = logL) 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 L and ℓ coincided. Hence, the estimates of β maximized ℓ(0,σ 2 | ε p ) from Eqs (7) and (8) knowing that ε p = f(y obs, y est ) and y est = g(β,x).
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 σ 2 A ℓ B σ 2 B ℓ C σ 2 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 σ 2 A had a far greater weight than σ 2 B or σ 2 C , with its optimization even driving C into a worsening linefit. To overcome this bias, the optimization of σ 2 was standardized to its value at the previous iteration, which only required replacing n by nσ 2 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 {_ x,_ y}. In cartesian coordinates, each residual was defined as vector A = {x-_ x,y-_ 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 = {δ,δÁ@ _ y/@ _ x}, by preference using small δ. When vectors A and B were perpendicular, their inner product was null i.e, A•B = 0. Hence, _ x and _ y were determined nesting Eq (19).
The calculus was performed using matrix algebra. The vectors with the y components of A and B were defined as A y = {y,-β 0 ,- The inner product was re-written as: The Newton-Raphson method generally took four iterations of _ to converge into the roots of Eq (20) with a change rate<0,1%. The derivative f'(_ x) was given by Eq (21), which required the derivatives of A y (Eq 22) and B y (Eq 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.

Testing differences between life-cycle stages
We tested whether the parameters s max , V opt , b 0 and b 1 changed significantly with ploidy or sex using non-parametric permutation tests with 1000 randomizations. The test accuracy was improved by keeping the identity of the observations relatively to the non-tested properties i.e, observations were permuted within their site and season [17,18]. A suited definition of the survival function required merging sites and seasons. Thus, differences between treatments within these factors could not be tested. The sample size was 35 and one degree of freedom was lost with the estimation of each of the four model parameters. Because V opt was the only parameter whose estimation was independent, the ANCOVA approach was followed i.e, first we estimated whether V opt was significantly different among ploidies or sexes. Only in case it was not, did we estimate whether τ max , b 0 and b 1 were significantly different among ploidies or sexes.

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][16][17][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 densitydependent 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.
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 σ 2 A decreased from 335Á10 −4 to 187Á10 −4 , σ 2 B increased from 73Á10 −4 to 155Á10 −4 and σ 2 C increased from 117Á10 −4 to 123Á10 −4 . 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][26][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][35][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][24][25][26][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.
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 0 ) and bigger asymmetry between competition and Allee effects (b 1 ). 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 S1 File. S1 software and data.zip. Software package and Gracilaria chilensis data. (ZIP)