Fits a multiple linear regression model using least squares.

Routine RLSE fits a multiple linear regression model with or without an intercept. If INTCEP = 1, the multiple linear regression model is

where the observed values of the yi’s (input in Y) constitute the responses or values of the dependent variable, the xi1’s, xi2’s, …, xik’s (input in X) are the settings of the k (input in NIND) independent variables, β0, β1, …, βk are the regression coefficients whose estimated values are output in B, and the ɛ i’s are independently distributed normal errors each with mean zero and variance σ2. Here, n is the number of valid rows in the augmented matrix (X, Y), i.e. n equals NOBS ‑ NRMISS (the number of rows that do not contain NaN). If INTCEP = 0, β0 is not included in the model.

Routine RLSE computes estimates of the regression coefficients by minimizing the sum of squares of the deviations of the observed response yi from the fitted response

for the n observations. This minimum sum of squares (the error sum of squares) is output and denoted by

In addition, the total sum of squares is output. For the case, INTCEP = 1, the total sum of squares is the sum of squares of the deviations of yi from its mean

— the so-called corrected total sum of squares; it is denoted by

For the case INTCEP = 0, the total sum of squares is the sum of squares of yi — the so-called uncorrected total sum of squares; it is denoted by

See Draper and Smith (1981) for a good general treatment of the multiple linear regression model, its analysis, and many examples.

In order to compute a least-squares solution, RLSE performs an orthogonal reduction of the matrix of regressors to upper triangular form. If the user needs the upper triangular matrix output for subsequent computing, the routine R2SE can be invoked in place of RLSE. (See the description of R in Comment 1). The reduction is based on one pass through the rows of the augmented matrix (X, Y) using fast Givens transformations. (See routines SROTMG and SROTM Golub and Van Loan, 1983, pages 156-162, Gentleman, 1974.) This method has the advantage that the loss of accuracy resulting from forming the crossproduct matrix used in the normal equations is avoided.

An orthogonal reduction on the centered matrix is computed. When the final computations are performed, the first row of R and the first element of B are updated so that they reflect the statistics for the original (uncentered) data. This means that the estimate of the intercept and the R matrix are for the uncentered data.

As part of the final computations, RLSE checks for linearly dependent regressors. If the i-th regressor is a linear combination of the first i ‑ 1 regressors, the i-th diagonal element of R is close to zero (exactly zero if infinite precision arithmetic could be used) prior to the final computations. In particular, linear dependence of the regressors is declared if any of the following three conditions is satisfied:

Here, Ri⋅1, 2, …, i−1 is the multiple correlation coefficient of the i-th independent variable with the first i ‑ 1 independent variables. If no intercept is in the model (INTCEP = 0), the “multiple correlation” coefficient is computed without adjusting for the mean.

On completion of the final computations, if the i-th regressor is declared to be linearly dependent upon the previous i ‑ 1 regressors, then the i-th element of B and all elements in the i-th row of R are set to zero. Finally, if a linear dependence is declared, an informational (error) message, code 1, is issued indicating the model is not full rank.

A regression model

is fitted to data taken from Maindonald (1984, pages 203‑204).

A weighted least-squares fit is computed using the model

and weights 1/i2 discussed by Maindonald (1984, pages 67 - 68). In order to compute the weighted least-squares fit, using an ordinary least squares routine (RLSE), the regressors (including the column of ones for the intercept term as well as the independent variables) and the responses must be transformed prior to invocation of RLSE. The transformed regressors and responses can be computed by using routine SHPROD (IMSL MATH/LIBRARY). For the i-th case the corresponding response and regressors are multiplied by a square root of the i-th weight. Because the column of ones corresponding to the intercept term in the untransformed model, is transformed by the weights, this transformed column of ones must be input to the least squares subroutine as an additional independent variable along with the option INTCEP = 0.

In terms of the original, untransformed regressors and responses, the minimum sum of squares for error output in SSE is

where here the weight wi = 1/i2. Also, since INTCEP = 0, the uncorrected total sum of squares is output in SST. In terms of the original untransformed responses,