Regression and Curve Fitting

Regression and Curve Fitting

Regression and Curve Fitting

Least squares

A regression model states that each observation is y_i = f(x_i; β) + ε_i, a systematic part with unknown coefficients β and an error ε_i. Least squares chooses β to minimize the sum of squared residuals S(β) = Σ w_i (y_i − f(x_i; β))², with weights w_i = 1 when every observation is equally precise. For a model that is linear in the coefficients (a straight line, a polynomial, y = b₀ + b₁ ln x, several predictors) the minimum is the solution of the normal equations XᵀWXβ = XᵀWy, with X the design matrix whose rows are the regressors of each observation. The tool does not form these equations: it factors √W X by Householder QR after scaling each column to unit length, which keeps the full working precision even on badly conditioned data such as the Longley set, and it fits polynomials on the centred and scaled predictor u = (x − x̄)/s_x before transforming the coefficients and their covariance back to the original x. Weights are used as 1/variance when the precision of each observation is known, or as replicate counts: a weight of 3 gives the same fit as the observation entered three times.

The inference that follows rests on four assumptions: the model form is right (the mean of y is f), the errors have constant variance (or a variance inversely proportional to the entered weights), they are independent, and they are normally distributed. The residual plots are the check on the first three; the Q–Q plot on the fourth. The results page gives no robust or autocorrelation-corrected alternative: when the diagnostics fail, the remedies are a transformation of y, a different model, weights, or a time-series method.

Standard errors, t tests and confidence intervals

The residual standard error s = √(SS_E / (n − p)) estimates the error standard deviation (for unit weights), with n observations and p coefficients. The covariance of the coefficients is s² (XᵀWX)⁻¹, and the square roots of its diagonal are the standard errors. The statistic t = b / SE tests whether a coefficient could be zero; its p-value is the two-sided tail probability of Student's t with n − p degrees of freedom. The confidence interval of a coefficient is b ± t_{α/2, n−p} SE at the chosen level 1 − α. A small p-value for the slope says that the relation is real; it says nothing about how strong or how useful it is, which R² and the prediction interval do.

R² and adjusted R²

The total sum of squares about the mean, SS_T = Σ w_i (y_i − ȳ)², splits into the part explained by the regression, SS_R, and the residual part, SS_E; R² = 1 − SS_E/SS_T is the explained fraction. It never falls when a coefficient is added, so the adjusted value R²_adj = 1 − (SS_E/(n − p)) / (SS_T/(n − 1)) is the one to compare models of different size on the same data; a term that does not earn its keep lowers it. For a model fitted through the origin the total sum of squares is taken about zero (Σ w_i y_i²), which makes that R² larger and not comparable with the intercept model's. The Akaike information criterion AIC = n ln(2π SS_E/n) + n + 2(p + 1) (the form of R's lm) is another size-penalized comparison; only differences between models fitted to the same y mean anything.

The analysis-of-variance table

The ANOVA table lists SS_R with p − 1 degrees of freedom (p for a model through the origin), SS_E with n − p, and SS_T with n − 1; the mean squares are the sums of squares divided by their degrees of freedom, and F = MS_R/MS_E tests whether the regression explains more than chance would. Its p-value is the tail probability of the F distribution with those degrees of freedom. For a straight line F = t² of the slope and the two tests are the same.

Confidence interval of the mean and prediction interval

At a chosen x₀ the fitted value ŷ₀ = x₀ᵀb estimates the mean of y there. Its standard error is s √(x₀ᵀ (XᵀWX)⁻¹ x₀) (for a straight line s √(1/n + (x₀ − x̄)²/S_xx)), and ŷ₀ ± t_{α/2, n−p} SE_mean is the confidence interval of the mean response. A single new observation at x₀ scatters about that mean with its own standard deviation s as well, so the prediction interval uses s √(1 + x₀ᵀ (XᵀWX)⁻¹ x₀) and is always wider; it does not shrink to zero with more data, because the scatter of one observation is a property of the process, not of the sample. The bands drawn in Figure 1 are these two intervals evaluated along x; both are narrowest at the mean of x and widen with the distance from it, so an extrapolation far outside the data carries a wide interval even when R² is high. A new observation is taken with unit weight; for an observation of weight w₀ the first term under the root would be 1/w₀.

Residual diagnostics

The leverage h_i = w_i x_iᵀ (XᵀWX)⁻¹ x_i is the diagonal of the hat matrix: how much observation i pulls its own fitted value. The leverages sum to p, so a value above 2p/n marks an observation far from the others in x, which the fit cannot afford to be wrong about. The standardized residual r_i = e_i √w_i / (s √(1 − h_i)) is the residual in units of its own standard deviation; under the model about 95 percent of them lie within ±2, and a value beyond ±3 is a candidate outlier to be explained, not deleted. Figure 2 plots r against the fitted values: a structureless horizontal band is what the model predicts; a bow means a missing curvature term, a funnel means a variance that grows with the response (a transformation of y, or weights, is the remedy). Figure 3 is the normal Q–Q plot, the ordered r_i against the normal quantiles of the plotting positions (i − 3/8)/(n + 1/4); a straight line supports the normal-error assumption, a bend at the ends means heavy tails, an S-shape skewness. The Durbin–Watson statistic Σ(e_i − e_i−1)² / Σe_i² in the order the data were entered is near 2 when successive residuals are independent, below about 1.5 when they are positively correlated (common when x is time), above 2.5 when negatively; it is meaningful only when the data order means something.

Transformations and why their R² is not comparable

An exponential law y = b₀ e^b₁x becomes a straight line in ln y, and a power law y = b₀ x^b₁ a straight line in ln y against ln x. Fitting that line is quick and is how the tool finds its starting values, but it minimizes the squared errors of ln y, which weights the small values of y heavily and the large ones lightly, and the R² it reports is the explained fraction of ln y, a different quantity from the R² of y. The tool therefore finishes every exponential and power fit by nonlinear least squares on the original y and reports R² = 1 − SS_E/SS_T on that scale; it prints the log-linear estimate and its R² beside them so the difference can be seen. If the errors are in fact proportional to y (constant percentage scatter) the log-scale fit is the better model, and the right way to obtain it here is the exponential model with weights 1/y_i², or a linear fit of the logged data pasted directly. The logarithmic model y = b₀ + b₁ ln x transforms only x, so its R² is on the y scale and is comparable.

Nonlinear least squares and Levenberg–Marquardt

When the model is not linear in its parameters there is no closed solution; S(θ) is minimized by iteration from starting values. Gauss–Newton linearizes the model about the current θ with the Jacobian J (∂f_i/∂θ_j) and solves the linear least-squares step JᵀWJ δ = JᵀW(y − f). Far from the minimum that step can overshoot; Marquardt's method adds λ diag(JᵀWJ) to the left-hand side, which turns the step toward the steepest descent direction and shortens it as λ grows. The tool accepts a step when it lowers S and then divides λ by 10; otherwise it multiplies λ by 10 and tries again. Iteration stops when the relative decrease of S falls below 10⁻¹² and the step is below 10⁻⁸ of the parameter size, or when no step lowers S any further; the table on the results page shows the path. The exponential and power models use analytic derivatives; a user model is differentiated by central differences with a step of 10⁻⁶ of each parameter, which costs 2p extra evaluations per iteration and gives derivatives good to about eight digits.

At the solution the model is again linearized, and the asymptotic covariance s² (JᵀWJ)⁻¹ gives the standard errors, t statistics and intervals exactly as in the linear case, with s² = S/(n − p). They are approximations whose quality depends on how linear the model is over the width of the intervals; for a well-determined fit (the NIST Misra1a data, on which the tool reproduces the certified standard errors) they are good, for a fit with strongly correlated parameters the true confidence region is banana-shaped and wider than the ± intervals suggest. The confidence and prediction bands are not drawn for nonlinear models for that reason; the prediction at x₀ uses the gradient of the model there in place of the regressor row and is labelled as linearized.

Starting values

Levenberg–Marquardt finds the minimum nearest its start; a poor start can converge to a local minimum or run the parameters to a region where the model cannot be evaluated. Good starts come from the physics (an asymptote, a rate constant estimated from a half-life), from a transformation that makes the model linear (as the tool does for the exponential and power laws), or from a few points solved exactly. A fit that stops at the iteration limit, or whose standard errors are of the size of the estimates, is reporting a problem with the model or the start, not a result. The NIST Statistical Reference Datasets give two certified starts for every nonlinear problem; the user-model default (Misra1a, start 500 and 0.0001) is one of them, and the fit reaches the certified values from both.

Reading the results page

The fitted equation is written out with its coefficients; the coefficient table gives the estimate, standard error, t, p-value and the interval at the chosen level; the goodness-of-fit block gives R², adjusted R², s, F with its p-value, AIC and Durbin–Watson; the ANOVA table follows; the residual table lists every observation with its fitted value, residual, standardized residual and leverage, and flags those beyond the usual limits; the prediction block gives ŷ₀ with both intervals and the equations they come from. For a polynomial the coefficients are reported in the original x, so their standard errors are large and strongly correlated by nature (that is the reason for fitting on the centred predictor); judge the polynomial by its fit and its bands, not by the t of its individual terms.

Limits

No robust fitting (an outlier is flagged, not down-weighted), no correction for autocorrelated errors (Durbin–Watson reports it), no model selection or stepwise procedures, no errors-in-variables (x is taken as known), no inference for weighted fits beyond the stated unit-weight prediction interval, and no confidence bands for nonlinear models. Up to 5,000 observations, polynomials to degree 6, eight parameters in a user model.

Verification

The library is checked by test-scripts/test_regression.php against the NIST Statistical Reference Datasets: Norris (straight line, the page's default data), Pontius (quadratic), Longley (six predictors, the classic ill-conditioned set) and Misra1a (nonlinear, both certified starts), every certified estimate, standard error, residual standard deviation, R² and ANOVA entry agreeing to at least nine significant figures; the interval formulas against the Norris closed forms and the rocket-propellant examples of Montgomery, Peck and Vining; and the t, F and normal quantiles against the standard tables.