IMSL C Math Library
regression
Fits a multiple linear regression model using least squares.
Synopsis
#include <imsl.h>
float *imsl_f_regression (int n_observations, int n_independent, float x[], float y[], …, 0)
The type double function is imsl_d_regression.
Required Arguments
int n_observations (Input)
The number of observations.
int n_independent (Input)
The number of independent (explanatory) variables.
float x[] (Input)
Array of size n_observations × n_independent containing the matrix of independent (explanatory) variables.
float y[] (Input)
Array of length n_observations containing the dependent (response) variable.
Return Value
If the optional argument IMSL_NO_INTERCEPT is not used, imsl_f_regression returns a pointer to an array of length n_independent + 1 containing a least-squares solution for the regression coefficients. The estimated intercept is the initial component of the array.
Synopsis with Optional Arguments
#include <imsl.h>
float *imsl_f_regression (int n_observations, int n_independent, float x[], float y[],
IMSL_X_COL_DIM, int x_col_dim,
IMSL_NO_INTERCEPT,
IMSL_TOLERANCE, float tolerance,
IMSL_RANK, int *rank,
IMSL_COEF_COVARIANCES, float **p_coef_covariances,
IMSL_COEF_COVARIANCES_USER, float coef_covariances[],
IMSL_COV_COL_DIM, int cov_col_dim,
IMSL_X_MEAN, float **p_x_mean,
IMSL_X_MEAN_USER, float x_mean[],
IMSL_RESIDUAL, float **p_residual,
IMSL_RESIDUAL_USER, float residual[],
IMSL_ANOVA_TABLE, float **p_anova_table,
IMSL_ANOVA_TABLE_USER, float anova_table[],
IMSL_RETURN_USER, float coefficients[],
0)
Optional Arguments
IMSL_X_COL_DIM, int x_col_dim (Input)
The column dimension of x.
Default: x_col_dim = n_independent
IMSL_NO_INTERCEPT
By default, the fitted value for observation i is
where k = n_independent. If IMSL_NO_INTERCEPT is specified, the intercept term
is omitted from the model.
IMSL_TOLERANCE, float tolerance (Input)
The tolerance used in determining linear dependence. For imsl_f_regression, tolerance = 100 ×imsl_f_machine(4) is the default choice. For imsl_d_regression, tolerance = 100 × imsl_d_machine(4) is the default. See imsl_f_machine.
IMSL_RANK, int *rank (Output)
The rank of the fitted model is returned in *rank.
IMSL_COEF_COVARIANCES, float **p_coef_covariances (Output)
The address of a pointer to the m × m array containing the estimated variances and covariances of the estimated regression coefficients. Here, m is the number of regression coefficients in the model. If IMSL_NO_INTERCEPT is specified, m = n_independent; otherwise, m = n_independent + 1. On return, the pointer is initialized (through a memory allocation request to malloc), and the array is stored there. Typically, float *p_coef_covariances is declared; &p_coef_covariances is used as an argument to this function; and imsl_free(p_coef_covariances) is used to free this array.
IMSL_COEF_COVARIANCES_USER, float coef_covariances[] (Output)
If specified, coef_covariances is an array of length m × m containing the estimated variances and covariances of the estimated coefficients where m is the number of regression coefficients in the model.
IMSL_COV_COL_DIM, int cov_col_dim (Input)
The column dimension of array coef_covariance.
Default: cov_col_dim = m where m is the number of regression coefficients in the model.
IMSL_X_MEAN, float **p_x_mean (Output)
The address of a pointer to the array containing the estimated means of the independent variables. On return, the pointer is initialized (through a memory allocation request to malloc), and the array is stored there. Typically, float *p_x_mean is declared; &p_x_mean is used as an argument to this function; and imsl_free(p_x_mean) is used to free this array.
IMSL_X_MEAN_USER, float x_mean[] (Output)
If specified, x_mean is an array of length n_independent provided by the user. On return, x_mean contains the means of the independent variables.
IMSL_RESIDUAL, float **p_residual (Output)
The address of a pointer to the array containing the residuals. On return, the pointer is initialized (through a memory allocation request to malloc), and the array is stored there. Typically, float *p_residual is declared; &p_residual is used as argument to this function; and imsl_free(p_residual) is used to free this array.
IMSL_RESIDUAL_USER, float residual[] (Output)
If specified, residual is an array of length n_observations provided by the user. On return, residual contains the residuals.
IMSL_ANOVA_TABLE, float **p_anova_table (Output)
The address of a pointer to the array containing the analysis of variance table. On return, the pointer is initialized (through a memory allocation request to malloc), and the array is stored there. Typically, float *p_anova_table is declared; &p_anova_table is used as argument to this function; and imsl_free(p_anova_table) is used to free this array.
The analysis of variance statistics are given as follows:
Element
Analysis of Variance Statistics
0
degrees of freedom for the model
1
degrees of freedom for error
2
total (corrected) degrees of freedom
3
sum of squares for the model
4
sum of squares for error
5
total (corrected) sum of squares
6
model mean square
7
error mean square
8
overall F-statistic
9
p-value
10
R2 (in percent)
11
adjusted R2 (in percent)
12
estimate of the standard deviation
13
overall mean of y
14
coefficient of variation (in percent)
IMSL_ANOVA_TABLE_USER, float anova_table[] (Output)
If specified, the 15 analysis of variance statistics listed above are computed and stored in the array anova_table provided by the user.
IMSL_RETURN_USER, float coefficients[] (Output)
If specified, the least-squares solution for the regression coefficients is stored in array coefficients provided by the user. If IMSL_NO_INTERCEPT is specified, the array requires m = n_independent units of memory; otherwise, the number of units of memory required to store the coefficients is m = n_independent + 1.
Description
The function imsl_f_regression fits a multiple linear regression model with or without an intercept. By default, the multiple linear regression model is
yi =β0 +β1xi1 +β2xi2 +…+βkxik +ɛi i =1, 2, …, n
where the observed values of the yi’s (input in y) are 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 n_independent) independent variables; β0, β1, …, βk are the regression coefficients whose estimated values are to be output by imsl_f_regression; and the ɛi’s are independently distributed normal errors each with mean zero and variance σ2. Here, n is the number of rows in the augmented matrix (x,y), i.e., n equals n_observations. Note that by default, β0 is included in the model.
The function imsl_f_regression 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 as one of the analysis of variance statistics if IMSL_ANOVA_TABLE (or IMSL_ANOVA_TABLE_USER) is specified and is computed as
Another analysis of variance statistic is the total sum of squares. By default, 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. This statistic is computed as
When IMSL_NO_INTERCEPT is specified, the total sum of squares is the sum of squares of yi, the so-called uncorrected total sum of squares. This is computed as
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, imsl_f_regression performs an orthogonal reduction of the matrix of regressors to upper-triangular form. The reduction is based on one pass through the rows of the augmented matrix (x, y) using fast Givens transformations. (See Golub and Van Loan 1983, pp. 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.
By default, the current means of the dependent and independent variables are used to internally center the data for improved accuracy. Let xi be a column vector containing the j-th row of data for the independent variables. Let represent the mean vector for the independent variables given the data for rows 1, 2, …, i. The current mean vector is defined to be
The i-th row of data has subtracted from it and is then weighted by i/(i − 1). Although a crossproduct matrix is not computed, the validity of this centering operation can be seen from the following formula for the sum of squares and crossproducts matrix:
An orthogonal reduction on the centered matrix is computed. When the final computations are performed, the intercept estimate and the first row and column of the estimated covariance matrix of the estimated coefficients are updated (if IMSL_COEF_COVARIANCES or IMSL_COEF_COVARIANCES_USER is specified) to reflect the statistics for the original (uncentered) data. This means that the estimate of the intercept is for the uncentered data.
As part of the final computations, imsl_regression checks for linearly dependent regressors. In particular, linear dependence of the regressors is declared if any of the following three conditions are satisfied:
*A regressor equals zero.
*Two or more regressors are constant.
* is less than or equal to tolerance. 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, 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 coefficient estimate and all elements in the i-th row and i-th column of the estimated variance-covariance matrix of the estimated coefficients (if IMSL_COEF_COVARIANCES or IMSL_COEF_COVARIANCES_USER is specified) are set to zero. Finally, if a linear dependence is declared, an informational (error) message, code IMSL_RANK_DEFICIENT, is issued indicating the model is not full rank.
Examples
Example 1
A regression model
is fitted to data taken from Maindonald (1984, pp. 203-204).
 
#include <imsl.h>
 
#define INTERCEPT 1
#define N_INDEPENDENT 3
#define N_COEFFICIENTS (INTERCEPT + N_INDEPENDENT)
#define N_OBSERVATIONS 9
 
int main()
{
float *coefficients;
float x[][N_INDEPENDENT] = {7.0, 5.0, 6.0,
2.0,-1.0, 6.0,
7.0, 3.0, 5.0,
-3.0, 1.0, 4.0,
2.0,-1.0, 0.0,
2.0, 1.0, 7.0,
-3.0,-1.0, 3.0,
2.0, 1.0, 1.0,
2.0, 1.0, 4.0};
float y[] = {7.0,-5.0, 6.0, 5.0, 5.0, -2.0, 0.0, 8.0, 3.0};
 
coefficients = imsl_f_regression(N_OBSERVATIONS, N_INDEPENDENT,
(float *)x, y, 0);
imsl_f_write_matrix("Least-Squares Coefficients", 1, N_COEFFICIENTS,
coefficients,
IMSL_COL_NUMBER_ZERO,
0);
}
Output
 
Least-Squares Coefficients
0 1 2 3
7.733 -0.200 2.333 -1.667
 
 
Example 2
A weighted least-squares fit is computed using the model
yi =β0xi0 +β1xi1 +β2xi2 +ɛi i =1, 2, …, 4
and weights 1/i2 discussed by Maindonald (1984, pp. 67-68). In order to compute the weighted least-squares fit, using an ordinary least-squares function (imsl_f_regression), the regressors (including the column of ones for the intercept term) and the responses must be transformed prior to invocation of imsl_f_regression. Specifically, the i-th response and regressors are multiplied by a square root of the i-th weight. IMSL_NO_INTERCEPT must be specified since the column of ones corresponding to the intercept term in the untransformed model is transformed by the weights and is regarded as an additional independent variable.
In the example, IMSL_ANOVA_TABLE is specified. The minimum sum of squares for error in terms of the original untransformed regressors and responses for this weighted regression is
where wi = 1/i2. Also, since IMSL_NO_INTERCEPT is specified, the uncorrected total sum-of-squares terms of the original untransformed responses is
 
#include <imsl.h>
#include <math.h>
 
#define N_INDEPENDENT 3
#define N_COEFFICIENTS N_INDEPENDENT
#define N_OBSERVATIONS 4
 
int main()
{
int i, j;
float *coefficients, w, anova_table[15], power;
float x[][N_INDEPENDENT] = {1.0, -2.0, 0.0,
1.0, -1.0, 2.0,
1.0, 2.0, 5.0,
1.0, 7.0, 3.0};
float y[] = {-3.0, 1.0, 2.0, 6.0};
char *anova_row_labels[] = {
"degrees of freedom for regression",
"degrees of freedom for error",
"total (uncorrected) degrees of freedom",
"sum of squares for regression",
"sum of squares for error",
"total (uncorrected) sum of squares",
"regression mean square",
"error mean square", "F-statistic",
"p-value", "R-squared (in percent)",
"adjusted R-squared (in percent)",
"est. standard deviation of model error",
"overall mean of y",
"coefficient of variation (in percent)"};
power = 0.0;
for (i = 0; i < N_OBSERVATIONS; i++) {
power += 1.0;
/* The square root of the weight */
w = sqrt(1.0 / (power*power));
/* Transform response */
y[i] *= w;
/* Transform regressors */
for (j = 0; j < N_INDEPENDENT; j++)
x[i][j] *= w;
}
 
coefficients = imsl_f_regression(N_OBSERVATIONS, N_INDEPENDENT,
(float *)x, y,
IMSL_NO_INTERCEPT,
IMSL_ANOVA_TABLE_USER,
anova_table, 0);
 
imsl_f_write_matrix("Least-Squares Coefficients", 1,
N_COEFFICIENTS, coefficients, 0);
imsl_f_write_matrix("* * * Analysis of Variance * * *\n", 15, 1,
anova_table, IMSL_ROW_LABELS, anova_row_labels,
IMSL_WRITE_FORMAT, "%10.2f", 0);
}
Output
 
Least-Squares Coefficients
1 2 3
-1.431 0.658 0.748
 
* * * Analysis of Variance * * *
 
degrees of freedom for regression 3.00
degrees of freedom for error 1.00
total (uncorrected) degrees of freedom 4.00
sum of squares for regression 10.93
sum of squares for error 1.01
total (uncorrected) sum of squares 11.94
regression mean square 3.64
error mean square 1.01
F-statistic 3.60
p-value 0.37
R-squared (in percent) 91.52
adjusted R-squared (in percent) 66.08
est. standard deviation of model error 1.01
overall mean of y -0.08
coefficient of variation (in percent) -1207.73
Warning Errors
IMSL_RANK_DEFICIENT
The model is not full rank. There is not a unique least-squares solution.