17.3 Estimation Methods for Non-Linear Regression

Non-linear regression models capture complex relationships between variables that can't be described by straight lines. These models use curved functions like exponentials or logarithms to fit data more accurately in many real-world scenarios.

Because these models are non-linear in their parameters, you can't solve for the best-fit values with simple algebra the way you can in ordinary linear regression. Instead, you need specialized estimation methods that search for optimal parameters through iteration. This topic covers how those methods work, when they succeed, and how to interpret what they produce.

Least squares estimation for non-linear models

Concept and application of least squares in non-linear regression

The core idea is the same as in linear regression: find parameter values that minimize the sum of squared residuals (SSR) between observed and predicted values. The difference is that in non-linear models, the SSR surface is no longer a simple bowl shape with a single, analytically solvable minimum.

• The objective function is still $S(\boldsymbol{\theta}) = \sum_{i=1}^{n} \left[ y_i - f(x_i, \boldsymbol{\theta}) \right]^2$, where $f(x_i, \boldsymbol{\theta})$ is a non-linear function of the parameter vector $\boldsymbol{\theta}$
• Because $f$ is non-linear in $\boldsymbol{\theta}$, taking derivatives and setting them to zero doesn't yield a closed-form solution. You have to minimize iteratively.
• Non-linear model forms include exponential growth ($y = \alpha e^{\beta x}$), logistic curves, power functions, and trigonometric functions, among others
• Example: In a population growth model $P(t) = \frac{K}{1 + e^{-r(t - t_0)}}$, least squares estimation adjusts $K$, $r$, and $t_0$ iteratively to minimize the squared differences between observed and predicted population sizes

Role of initial parameter values and iterative optimization

Unlike linear least squares, where the solution is unique and directly computable, non-linear least squares depends heavily on where you start.

• The SSR surface for a non-linear model can have multiple local minima. If your starting values land you near a local (but not global) minimum, the algorithm may converge there and miss the true best fit.
• Choosing good initial values matters a lot. Common strategies include:
  • Using domain knowledge (e.g., a rough estimate of carrying capacity from the data)
  • Fitting a linearized version of the model first to get ballpark estimates
  • Running the algorithm from several different starting points and comparing results
• Example: In a logistic growth model, if you set the initial carrying capacity far below the observed data range, the algorithm may converge to a poor solution. Starting with $K$ near the maximum observed value is a much better choice.

Iterative methods for non-linear estimation

Gauss-Newton and Levenberg-Marquardt algorithms

These are the two most widely used algorithms for non-linear least squares. Both work by repeatedly linearizing the model around the current parameter estimates and solving the resulting linear problem.

Gauss-Newton method:

1. Start with initial parameter estimates $\boldsymbol{\theta}^{(0)}$
2. At each iteration $k$, compute the Jacobian matrix $\mathbf{J}$ (partial derivatives of the model with respect to each parameter, evaluated at the current estimates)
3. Approximate the Hessian as $\mathbf{J}^T \mathbf{J}$ (this avoids computing second derivatives)
4. Solve the update equation: $\Delta \boldsymbol{\theta} = (\mathbf{J}^T \mathbf{J})^{-1} \mathbf{J}^T \mathbf{r}$, where $\mathbf{r}$ is the vector of residuals
5. Update: $\boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} + \Delta \boldsymbol{\theta}$
6. Repeat until convergence

The Gauss-Newton method converges quickly when the model is mildly non-linear and the residuals are small, but it can become unstable when the approximation $\mathbf{J}^T \mathbf{J} \approx \mathbf{H}$ is poor.

Levenberg-Marquardt method:

This extends Gauss-Newton by adding a damping parameter $\lambda$ to the update equation:

$\Delta \boldsymbol{\theta} = (\mathbf{J}^T \mathbf{J} + \lambda \mathbf{I})^{-1} \mathbf{J}^T \mathbf{r}$

• When $\lambda$ is large, the update behaves like gradient descent (small, cautious steps)
• When $\lambda$ is small, the update behaves like Gauss-Newton (larger, faster steps)
• The algorithm adjusts $\lambda$ adaptively: it increases $\lambda$ if a step fails to reduce SSR, and decreases it when steps are successful

This makes Levenberg-Marquardt more robust than pure Gauss-Newton, which is why it's the default in most non-linear regression software.

Jacobian matrix and damping factor

• The Jacobian matrix $\mathbf{J}$ has dimensions $n \times p$ (observations by parameters). Each entry $J_{ij} = \frac{\partial f(x_i, \boldsymbol{\theta})}{\partial \theta_j}$ tells you how sensitive the model prediction for observation $i$ is to a small change in parameter $j$.
• Both algorithms require recalculating $\mathbf{J}$ at every iteration, since the partial derivatives change as $\boldsymbol{\theta}$ changes. For models without closed-form derivatives, numerical differentiation (finite differences) is used instead.
• The damping factor $\lambda$ in Levenberg-Marquardt controls the trade-off between convergence speed and stability. Too small and you risk overshooting; too large and convergence slows to a crawl. Most implementations handle this automatically, but understanding the trade-off helps when diagnosing convergence problems.

Convergence assessment and termination criteria

An iterative algorithm needs clear rules for when to stop. The three standard criteria are:

1. Parameter change criterion: Stop when the relative change in parameter estimates between iterations falls below a tolerance, e.g., $\frac{||\boldsymbol{\theta}^{(k+1)} - \boldsymbol{\theta}^{(k)}||}{||\boldsymbol{\theta}^{(k)}||} < 10^{-6}$

2. Objective function criterion: Stop when the reduction in SSR between iterations is negligibly small, e.g., less than $10^{-12}$

3. Safeguard criteria: Stop after a maximum number of iterations (e.g., 100 or 500) or a time limit to prevent runaway computation

If the algorithm hits the maximum iteration limit without satisfying the convergence tolerance, that's a warning sign. It may indicate poor starting values, an ill-conditioned problem, or a model that doesn't fit the data well.

Convergence and stability of non-linear methods

Factors influencing convergence rate

Several factors determine how quickly (or whether) an algorithm converges:

• Initial values: Starting closer to the true solution generally means fewer iterations. Multiple local minima make this especially important.
• Model complexity: Highly curved or overparameterized models create more difficult optimization landscapes. A model with many parameters relative to the data size is harder to fit.
• Data quality: Noisy data or datasets with outliers can distort the SSR surface, slowing convergence or pulling estimates toward misleading solutions.
• Parameter scaling: If parameters differ by orders of magnitude (e.g., one is around 0.001 and another around 10,000), the algorithm may struggle. Rescaling parameters to similar magnitudes often helps.

Stability and ill-conditioning

Stability means the algorithm consistently converges to the same (or very similar) solution despite small changes in starting values or data.

• An unstable method might give you very different parameter estimates when you drop a few data points or shift your starting values slightly. That's a red flag.
• Ill-conditioning of the $\mathbf{J}^T \mathbf{J}$ matrix is a common source of instability. This happens when columns of $\mathbf{J}$ are nearly linearly dependent, meaning two or more parameters have very similar effects on the model predictions. The matrix becomes nearly singular, and small numerical errors get amplified.
• Remedies for ill-conditioning include:
  • Reparameterization: Rewrite the model so parameters are more orthogonal (less correlated in their effects)
  • Tikhonov regularization: Add a small constant $\lambda$ to the diagonal of $\mathbf{J}^T \mathbf{J}$ (note that Levenberg-Marquardt already does this implicitly)
  • Centering or scaling the independent variables

Diagnostic tools for assessing convergence and stability

• Convergence plots show parameter estimates or SSR values across iterations. A smooth, monotonically decreasing SSR plot indicates healthy convergence. Oscillating or diverging traces suggest instability or a poorly chosen algorithm/starting point.
• Residual analysis after convergence works the same way as in linear regression: plot residuals vs. fitted values and check for randomness. A non-random pattern (e.g., a curve) suggests the non-linear model form itself may be wrong. Fanning patterns suggest heteroscedasticity.
• Profile likelihood plots can reveal whether the SSR surface is well-behaved (a clear single minimum) or problematic (flat regions, ridges, or multiple minima) for each parameter.

These diagnostics aren't optional extras. They're how you verify that your estimation actually worked.

Parameter interpretation in non-linear models

Meaning and interpretation of parameter estimates

Unlike linear models where each coefficient has a straightforward "unit change" interpretation, non-linear model parameters often have context-specific meanings that depend on the functional form.

• In a logistic growth model $P(t) = \frac{K}{1 + e^{-r(t-t_0)}}$, $K$ is the carrying capacity (the upper asymptote), $r$ is the intrinsic growth rate, and $t_0$ is the inflection point. Each parameter controls a different aspect of the curve's shape.
• The effect of changing one parameter typically depends on the values of the other parameters. This interdependence is a key difference from linear models.
• Always interpret estimates in the context of the specific model and the units of your variables.

Standard errors, confidence intervals, and hypothesis tests

Inference for non-linear models relies on large-sample (asymptotic) approximations, so these results are less exact than in linear regression, especially with small samples.

• Standard errors are derived from the variance-covariance matrix, estimated as $\hat{\sigma}^2 (\mathbf{J}^T \mathbf{J})^{-1}$, where $\hat{\sigma}^2$ is the residual mean square and $\mathbf{J}$ is evaluated at the converged estimates. These standard errors are approximate.
• Confidence intervals use the form $\hat{\theta}_j \pm t_{\alpha/2, \, n-p} \cdot SE(\hat{\theta}_j)$. A 95% confidence interval that excludes zero suggests the parameter is significantly different from zero at the 0.05 level.
• Hypothesis tests use the t-statistic $t = \frac{\hat{\theta}_j}{SE(\hat{\theta}_j)}$, compared against the t-distribution with $n - p$ degrees of freedom. If $|t|$ exceeds the critical value (e.g., approximately 1.96 for large samples at the 0.05 level), you reject the null that the parameter equals zero.

Statistical significance and variable importance

• A statistically significant parameter estimate indicates that the corresponding term contributes meaningfully to the model's fit. A non-significant estimate suggests the variable (or that particular non-linear term) may not be needed.
• Comparing raw parameter magnitudes across variables is misleading when variables are on different scales. Standardized estimates, computed by scaling each estimate by the ratio of the standard deviation of the corresponding independent variable to the standard deviation of the dependent variable, allow fairer comparisons.
• Example: In a non-linear crop yield model, if the standardized estimate for soil moisture is 0.45 and for fertilizer application is 0.20, soil moisture has a stronger relative influence on yield.
• Be cautious about dropping non-significant terms from non-linear models without careful thought. Removing a parameter can fundamentally change the shape of the fitted function, not just shift a line as in linear regression.