R - Nonlinear Least Square

R - Nonlinear Least Squares

While linear regression is powerful, many real-world phenomena don’t follow a straight line. They are better described by nonlinear models, such as exponential growth curves in biology or power laws in economics. Nonlinear Least Squares (NLS) is a method for fitting a nonlinear model to a set of data. The goal, similar to linear regression, is to find the parameter values that minimize the sum of the squared differences between the observed data and the values predicted by the model.

In R, the standard function for NLS is nls(). The process involves defining a model formula, providing your data, and—most importantly—giving the algorithm reasonable starting estimates for the parameters. The function then iteratively refines these estimates to find the best fit.

Syntax

The basic syntax for nls() is:

nls(formula, data, start)

• formula: A nonlinear model formula, including variables and the parameters you want to estimate (e.g., y ~ a * x^b).
• data: A data frame containing the variables used in the formula.
• start: A named list or vector of initial guesses for the parameters (e.g., list(a = 1, b = 2)).

Example: Fitting a Quadratic Model

Let’s fit the model y = b1*x^2 + b2 to some sample data. A crucial step is choosing good starting values. We can often get a hint by plotting the data.

# Load modern packages for data manipulation and plotting
library(tidyverse)

# 1. Prepare the data
xvalues <- c(1.6, 2.1, 2, 2.23, 3.71, 3.25, 3.4, 3.86, 1.19, 2.21)
yvalues <- c(5.19, 7.43, 6.94, 8.11, 18.75, 14.88, 16.06, 19.12, 3.21, 7.58)

# Create a tibble (a modern data frame)
data_to_fit <- tibble(x = xvalues, y = yvalues)

# 2. Visualize the data to guess starting values
# The curve looks like a parabola opening upwards (x^2 term) and shifted up.
# Let's guess b1 ≈ 1 and the y-intercept b2 ≈ 2 or 3.
ggplot(data_to_fit, aes(x = x, y = y)) +
  geom_point(color = "blue", size = 3) +
  labs(title = "Data and Initial Guesses", x = "X Values", y = "Y Values")

# 3. Define starting values and fit the model
start_values <- list(b1 = 1, b2 = 3)

# Use nls() to find the best-fit parameters
nls_model <- nls(y ~ b1 * x^2 + b2, data = data_to_fit, start = start_values)

# 4. Examine the model summary
# summary() is more informative than just printing the model
print(summary(nls_model))

# 5. Get confidence intervals for the parameters
print(confint(nls_model))

The summary provides detailed information about the model fit:

Formula: y ~ b1 * x^2 + b2

Parameters:
   Estimate Std. Error t value Pr(>|t|)
_
b1  1.19235    0.05203  22.917 1.48e-08 ***
b2  1.99692    0.44858   4.452  0.00193 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3677 on 8 degrees of freedom

Number of iterations to convergence: 4
Achieved convergence tolerance: 4.887e-07

Waiting for profiling to be done...
       2.5%    97.5%
b1 1.072228 1.312471
b2 0.962804 3.031045

Visualizing the Fit with ggplot2

The best way to assess the model is to plot the fitted curve over the original data. We can use the predict() function to generate points from our model.

# Augment the original data with fitted values from the model
fit_data <- broom::augment(nls_model)

# Plot original points and the fitted NLS curve
ggplot(data = data_to_fit, aes(x = x, y = y)) +
  geom_point(color = "blue", size = 3, alpha = 0.7) +
  geom_line(data = fit_data, aes(y = .fitted), color = "red", size = 1) +
  labs(
    title = "Nonlinear Least Squares Fit",
    subtitle = "y = 1.19*x^2 + 2.00",
    x = "X Values",
    y = "Y Values"
  ) +
  theme_minimal()

Common Challenges: Convergence

A common issue with nls() is a failure to converge. This usually means one of two things:

1. Poor Starting Values: The initial guesses for the parameters were too far from the optimal values. Try plotting your data and making more educated guesses.
2. Inappropriate Model: The model formula you’ve chosen may be a fundamentally poor fit for the data’s underlying structure. Consider a different model.