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-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 3rd degree polynomial in order to density (x). Although it is non-linear, if each of the xn 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 (vi) of a cylinder of equal length and diameter as the frond, significantly correlated with dry weight (r2 = 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 (Vp = ∑Vi), rescaled to units of L∙m-2. The survival (s) dependency on V 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-b(V-V)2 defined the maximum survival (s) attained at an optimal biomass density (V). The parameter b represented the sensitivity of s to V, with higher b yielding steeper curves. However, the sensitivity of s was observed to be asymmetrical relative to lower vs higher V, meaning Allee effects and competition did not act proportionally. Therefore, b was replaced by b(1+b(V-V)), 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.

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.

The procedure to estimate the parameters s, V, b and b was to re-formulate Eq (1) as a 3rd degree polynomial in order to x = V i.e, s = β+βx+βx2+βx3, and determine the β 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 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-V-V)/(V-V) and y = (2s-s-s)/(s-s), with the left superscripts and 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 xi 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 xi (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 β = (XTX)-1XTy, 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 β = (XTWX)-1XTWy.

Iterative Re-weigthed Least Squares (IRLS)

The IRLS estimated β minimizing the weighted y-residuals (Eq 4) while re-estimating their weights (wp) 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 (ε = y-y) followed a normal distribution N(0,σ2). 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,σ2|ε) from Eqs (7) and (8) knowing that ε = f(y y) and y = g(β,x).

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, ∂ℓ/∂β1 = 0, ∂ℓ/∂β2 = 0 and ∂ℓ/∂β3 = 0) with 4 unknowns (β0, β1, β2 and β3). However, σ2 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 σ2 (Eq 9). This duality corresponded to a system of 8 equations, where the four ∂σ2/∂β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).

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 μ = ∑xm∙sk, · m = 0,…,6 ∧ k = 0,…,2. Matrix algebra estimated faster μ = {x0, x1, x2, x3, x4, x5, x6}T∙{y0, y1, y2}. It was also defined В = {1, -β, -β, -β, -β}. This way ℓ could be estimated from Eq (13) and σ2 from Eq (14), were n is the sample size and ⊗ is the Hadamard product i.e, matrix element-wise.

Their partial derivatives had the vectorial representation vi = ∂В/∂ci, 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).

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

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).

The calculus was performed using matrix algebra. The vectors with the y components of A and B were defined as A = {y,-β,-βẋ,-βẋ2,-βẋ3} and B = {β,2βẋ,3βẋ2}. The inner product was re-written as:

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 (Eq 22) and B (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.

Step 2: the magnitude of each oblique residual was given by |res| = ((x-ẋ)2+(y-Ẋβ)2)1/2, with the (1×4) design matrix Ẋ = {1, ẋ, ẋ, ẋ}. The minimum sum of their squares i.e, SSE = ∑res2, 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/∂β = -∑(y-Ẋβ)ẋi,∀ 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' = Ẋ(Ẋβ-y). The denominator SSE'' was the (4×4) Jacobian matrix with the second order partial derivatives ∂2SSE/∂βi∂βj = ∑ẋiẋj, ∀ i,j = 0,…,3. Matrix calculus estimated it from SSE'' = (ẊẊ). Merging all yielded the incremental βt+1 = βt-(ẊẊ)-1Ẋ(Ẋβ-y). Only a single iteration within step 1 was performed, after which the perpendicular coordinates of the residuals needed to be re-estimated.

Testing differences between life-cycle stages

We tested whether the parameters s, V, b and b 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 was the only parameter whose estimation was independent, the ANCOVA approach was followed i.e, first we estimated whether V was significantly different among ploidies or sexes. Only in case it was not, did we estimate whether τ, b and b 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-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.

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.

	smax			Vopt			b0			b1
	M	F	D	M	F	D	M	F	D	M	F	D
(value)	.724	.742	.661	4.261	3.898	3.88	.038	.039	.031	-.122	-.118	-.118
(p)	M	F	D	M	F	D	M	F	D	M	F	D
M		.732	.358		.147	.121		.898	.514		.658	.662
F			.283			.932			.484			.992

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) 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 σ2A decreased from 335∙10−4 to 187∙10−4, σ2B increased from 73∙10−4 to 155∙10−4 and σ2C 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-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.

Table 2

The overall survival function with parameters estimated by different methods.

Sample size n = 105.

	smax	Vopt	b0	b1
OLS	0.558	5.004	0.008	-0.059
WLS	0.588	4.432	0.014	-0.108
IRLS	0.561	4.489	0.021	-0.123
N-RM	0.558	5.004	0.008	-0.059
N-RM gp	0.578	4.786	0.011	-0.087
LOLS	0.701	4.009	0.035	-0.119

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 -1fits 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.

S1 software and data.zip.

Software package and Gracilaria chilensis data.

(ZIP)

Click here for additional data file.