Correlation analysis provides a single summary statistic – the correlation coefficient – describing the strength of the association between two variables. The most common types of association which are investigated by correlation analysis are linear relationships, and there are a number of forms of linear correlation coefficients for use with different types of data.

The (Pearson) product-moment correlation coefficients measure a linear relationship, while Kendall's tau and Spearman's rank order correlation coefficients measure monotonicity only. All three coefficients range from  to . A coefficient of zero always indicates that no linear relationship exists; a  coefficient implies a "perfect" positive relationship (i.e., an increase in one variable is always associated with a corresponding increase in the other variable); and a coefficient of  indicates a "perfect" negative relationship (i.e., an increase in one variable is always associated with a corresponding decrease in the other variable).

Consider the bivariate scattergrams in Figure 1: (a) and (b) show strictly linear functions for which the values of the product-moment correlation coefficient, and (since a linear function is also monotonic) both Kendall's tau and Spearman's rank order coefficients, would be  and  respectively. However, though the relationships in figures (c) and (d) are respectively monotonically increasing and monotonically decreasing, for which both Kendall's and Spearman's non-parametric coefficients would be  (in (c)) and  (in (d)), the functions are nonlinear so that the product-moment coefficients would not take such "perfect" extreme values. There is no obvious relationship between the variables in figure (e), so all three coefficients would assume values close to zero, while in figure (f) though there is an obvious parabolic relationship between the two variables, it would not be detected by any of the correlation coefficients which would again take values near to zero; it is important therefore to examine scattergrams as well as the correlation coefficients.

In order to decide which type of correlation is the most appropriate, it is necessary to appreciate the different groups into which variables may be classified. Variables are generally divided into four types of scales: the nominal scale, the ordinal scale, the interval scale, and the ratio scale. The nominal scale is used only to categorise data; for each category a name, perhaps numeric, is assigned so that two different categories will be identified by distinct names. The ordinal scale, as well as categorising the observations, orders the categories. Each category is assigned a distinct identifying symbol, in such a way that the order of the symbols corresponds to the order of the categories. (The most common system for ordinal variables is to assign numerical identifiers to the categories, though if they have previously been assigned alphabetic characters, these may be transformed to a numerical system by any convenient method which preserves the ordering of the categories.) The interval scale not only categorises and orders the observations, but also quantifies the comparison between categories; this necessitates a common unit of measurement and an arbitrary zero-point. Finally, the ratio scale is similar to the interval scale, except that it has an absolute (as opposed to arbitrary) zero-point.

For a more complete discussion of these four types of scales, and some examples, the user is referred to Churchman and Ratoosh (1959) and Hays (1970).

Product-moment correlation coefficients are used with variables which are interval (or ratio) scales; these coefficients measure the amount of spread about the linear least-squares equation. For a product-moment correlation coefficient, , based on  pairs of observations, testing against the null hypothesis that there is no correlation between the two variables, the statistic

has a Student's -distribution with  degrees of freedom; its significance can be tested accordingly.

Ranked and ordinal scale data are generally analysed by non-parametric methods – usually either Spearman's or Kendall's tau rank-order correlation coefficients, which, as their names suggest, operate solely on the ranks, or relative orders, of the data values. Interval or ratio scale variables may also be validly analysed by non-parametric methods, but such techniques are statistically less powerful than a product-moment method. For a Spearman rank-order correlation coefficient, , based on  pairs of observations, testing against the null hypothesis that there is no correlation between the two variables, for large samples the statistic

has approximately a Student's -distribution with  degrees of freedom, and may be treated accordingly. (This is similar to the product-moment correlation coefficient, , see above.) Kendall's tau coefficient, based on  pairs of observations, has, for large samples, an approximately Normal distribution with mean zero and standard deviation

when tested against the null hypothesis that there is no correlation between the two variables; the coefficient should therefore be divided by this standard deviation and tested against the standard Normal distribution, .

When the number of ordinal categories a variable takes is large, and the number of ties is relatively small, Spearman's rank-order correlation coefficients have advantages over Kendall's tau; conversely, when the number of categories is small, or there are a large number of ties, Kendall's tau is usually preferred. Thus when the ordinal scale is more or less continuous, Spearman's rank-order coefficients are preferred, whereas Kendall's tau is used when the data is grouped into a smaller number of categories; both measures do however include corrections for the occurrence of ties, and the basic concepts underlying the two coefficients are quite similar. The absolute value of Kendall's tau coefficient tends to be slightly smaller than Spearman's coefficient for the same set of data.

There is no authoritative dictum on the selection of correlation coefficients – particularly on the advisability of using correlations with ordinal data. This is a matter of discretion for the user.

The correlation coefficients described above measure the association between two variables ignoring any other variables in the system. Suppose there are three variables  and  as shown in the path diagram below.

 Figure 2

The association between  and  is made up of the direct association between  and  and the association caused by the path through , that is the association of both  and  with the third variable . For example if  and  were cholesterol level and blood pressure and  were age since both blood pressure and cholesterol level may increase with age the correlation between blood pressure and cholesterol level eliminating the effect of age is required.

The correlation between two variables eliminating the effect of a third variable is know as the partial correlation. If ,  and  represent the correlations between ,  and  then the partial correlation between  and  given  is

The partial correlation is then estimated by using product-moment correlation coefficients.

In general, let a set of variables be partitioned into two groups  and  with  variables in  and  variables in  and let the variance-covariance matrix of all  variables be partitioned into

Then the variance-covariance of  conditional on fixed values of the  variables is given by

The partial correlation matrix is then computed by standardising .

The product-moment correlation coefficient can be greatly affected by the presence of a few extreme observations or outliers. There are robust estimation procedures which aim to decrease the effect of extreme values.

Mathematically these methods can be described as follows. A robust estimate of the variance-covariance matrix, , can be written as

where  is a correction factor to give an unbiased estimator if the data is Normal and  is a lower triangular matrix. Let  be the vector of values for the th observation and let ,  being a robust estimate of location, then  and  are found as solutions to

and

where ,  and  are functions such that they return a value of 1 for reasonable values of  and decreasing values for large . The correlation matrix can then be calculated from the variance-covariance matrix. If , , and  returned 1 for all values then the product-moment correlation coefficient would be calculated.

In regression analysis the relationship between one specific random variable, the dependent or response variable, and one or more known variables, called the independent variables or covariates, is studied. This relationship is represented by a mathematical model, or an equation, which associates the dependent variable with the independent variables, together with a set of relevant assumptions. The independent variables are related to the dependent variable by a function, called the regression function, which involves a set of unknown parameters. Values of the parameters which give the best fit for a given set of data are obtained; these values are known as the estimates of the parameters.

The reasons for using a regression model are twofold. The first is to obtain a description of the relationship between the variables as an indicator of possible causality. The second reason is to predict the value of the dependent variable from a set of values of the independent variables. Accordingly, the most usual statistical problems involved in regression analysis are:

When the regression model is linear in the parameters (but not necessarily in the independent variables), then the regression model is said to be linear; otherwise the model is classified as nonlinear.

The most elementary form of regression model is the simple linear regression of the dependent variable, , on a single independent variable, , which takes the form

(1)

where  is the expected or average value of  and  and  are the parameters whose values are to be estimated, or, if the regression is required to pass through the origin (i.e., no constant term),

(2)

where  is the only unknown parameter.

An extension of this is multiple linear regression in which the dependent variable, , is regressed on the  () independent variables, , which takes the form

(3)

where  and  are the unknown parameters.

A special case of multiple linear regression is polynomial linear regression, in which the  independent variables are in fact powers of the same single variable  (i.e., , for ).

In this case, the model defined by (3) becomes

(4)

There are a great variety of nonlinear regression models; one of the most common is exponential regression, in which the equation may take the form

(5)

It should be noted that equation (4) represents a linear regression, since even though the equation is not linear in the independent variable, , it is linear in the parameters , whereas the regression model of equation (5) is nonlinear, as it is nonlinear in the parameters (,  and ).

The method used to determine values for the parameters is, based on a given set of data, to minimize the sums of squares of the differences between the observed values of the dependent variable and the values predicted by the regression equation for that set of data – hence the term least-squares estimation. For example, if a regression model of the type given by equation (3), namely

where  for all observations, is to be fitted to the  data points

(6)

such that

where  are unknown independent random errors with  and ,  being a constant, then the method used is to calculate the estimates of the regression parameters  by minimizing

(7)

If the errors do not have constant variance, i.e.,

then weighted least-squares estimation is used in which

is minimized. For a more complete discussion of these least-squares regression methods, and details of the mathematical techniques used, see Draper and Smith (1985) or Kendall and Stuart (1973).

One application of regression models is in the analysis of experiments. In this case the model relates the dependent variable to qualitative independent variables known as factors. Factors may take a number of different values known as levels. For example, in an experiment in which one of four different treatments is applied, the model will have one factor with four levels. Each level of the factor can be represented by a dummy variable taking the values 0 or 1. So in the example there are four dummy variables , for  such that:

along with a variable for the mean :

If there were 7 observations the data would be:

Models which include factors are sometimes known as General Linear (Regression) Models. When dummy variables are used it is common for the model not to be of full rank. In the case above, the model would not be of full rank because

This means that the effect of  cannot be distinguished from the combined effect of  and . This is known as aliasing. In this situation, the aliasing can be deduced from the experimental design and as a result the model to be fitted; in such situations it is known as intrinsic aliasing. In the example above no matter how many times each treatment is replicated (other than 0) the aliasing will still be present. If the aliasing is due to a particular data set to which the model is to be fitted then it is known as extrinsic aliasing. If in the example above observation 1 was missing then the  term would also be aliased. In general intrinsic aliasing may be overcome by changing the model, e.g., remove  or  from the model, or by introducing constraints on the parameters, e.g., .

If aliasing is present then there will no longer be a unique set of least-squares estimates for the parameters of the model but the fitted values will still have a unique estimate. Some linear functions of the parameters will also have unique estimates; these are known as estimable functions. In the example given above the functions () and () are both estimable.

In many situations there are several possible independent variables, not all of which may be needed in the model. In order to select a suitable set of independent variables, two basic approaches can be used.

Having fitted a model two questions need to be asked: first, "are all the terms in the model needed?" and second, "is there some systematic lack of fit?". To answer the first question either confidence intervals can be computed for the parameters or -tests can be calculated to test hypotheses about the regression parameters – for example, whether the value of the parameter, , is significantly different from a specified value,  (often zero). If the estimate of  is  and its standard error is  then the -statistic is

It should be noted that both the tests and the confidence intervals may not be independent. Alternatively -tests based on the residual sums of squares for different models can also be used to test the significance of terms in the model. If model 1, giving residual sum of squares  with degrees of freedom , is a sub-model of model 2, giving residual sum of squares  with degrees of freedom , i.e., all terms in model 1 are also in model 2, then to test if the extra terms in model 2 are needed the -statistic

may be used. These tests and confidence intervals require the additional assumption that the errors, , are Normally distributed.

To check for systematic lack of fit the residuals, , where  is the fitted value, should be examined. If the model is correct then they should be random with no discernable pattern. Due to the way they are calculated the residuals do not have constant variance. Now the vector of fitted values can be written as a linear combination of the vector of observations of the dependent variable, , . The variance-covariance matrix of the residuals is then ,  being the identity matrix. The diagonal elements of , , can therefore be used to standardize the residuals. The  are a measure of the effect of the th observation on the fitted model and are sometimes known as leverages.

If the observations were taken serially the residuals may also be used to test the assumption of the independence of the  and hence the independence of the observations.

Let  be the  by  matrix of independent variables and  be the vector of values for the dependent variable. To find the least-squares estimates of the vector of parameters, , the  decomposition of  is found, i.e.,

where ,  being a  by  upper triangular matrix, and  an  by  orthogonal matrix. If  is of full rank then  is the solution to

where  and  is the first  rows of . If  is not of full rank, a solution is obtained by means of a singular value decomposition (SVD) of ,

where  is a  by  diagonal matrix with non-zero diagonal elements,  being the rank of , and  and  are  by  orthogonal matrices. This gives the solution

 being the first  columns of  and  being the first  columns of .

This will be only one of the possible solutions. Other estimates may be obtained by applying constraints to the parameters. If weighted regression with a vector of weights  is required then both  and  are premultiplied by .

The method described above will, in general, be more accurate than methods based on forming (), (or a scaled version), and then solving the equations

Least-squares regression can be greatly affected by a small number of unusual, atypical, or extreme observations. To protect against such occurrences, robust regression methods have been developed. These methods aim to give less weight to an observation which seems to be out of line with the rest of the data given the model under consideration. That is to seek to bound the influence. For a discussion of influence in regression, see Hampel et al. (1986) and Huber (1981).

There are two ways in which an observation for a regression model can be considered atypical. The values of the independent variables for the observation may be atypical or the residual from the model may be large.

The first problem of atypical values of the independent variables can be tackled by calculating weights for each observation which reflect how atypical it is, i.e., a strongly atypical observation would have a low weight. There are several ways of finding suitable weights; some are discussed in Hampel et al. (1986).

The second problem is tackled by bounding the contribution of the individual  to the criterion to be minimized. When minimizing (7) a set of linear equations is formed, the solution of which gives the least-squares estimates. The equations are

These equations are replaced by

(8)

where  is the variance of the , and  is a suitable function which down weights large values of the standardized residuals . There are several suggested forms for , one of which is Huber's function,

(9)

The solution to (8) gives the -estimates of the regression coefficients. The weights can be included in (8) to protect against both types of extreme value. The parameter  can be estimated by the median absolute deviations of the residuals or as a solution to, in the unweighted case,

where  is a suitable function and  is a constant chosen to make the estimate unbiased.  is often chosen to be  where  is given in (9). Another form of robust regression is to minimize the sum of absolute deviations, i.e.,

For details of robust regression, see Hampel et al. (1986) and Huber (1981).

Robust regressions using least absolute deviations can be computed using routines in Chapter e02.

Generalized linear models are an extension of the general linear regression model discussed above. They allow a wide range of models to be fitted. These included certain non-linear regression models, logistic and probit regression models for binary data, and log-linear models for contingency tables. A generalized linear model consists of three basic components:

The model is fitted by maximum likelihood; this is equivalent to least-squares in the case of the Normal distribution. The residual sums of squares used in regression models is generalized to the concept of deviance. The deviance is the logarithm of the ratio of the likelihood of the model to the full model in which , where  is the estimated value of . For the Normal distribution the deviance is the residual sum of squares. Except for the case of the Normal distribution with the identity link, the  and -tests based on the deviance are only approximate; also the estimates of the parameters will only be approximately Normally distributed. Thus only approximate - or -tests may be performed on the parameter values and approximate confidence intervals computed.

The estimates are found by using an iterative weighted least-squares procedure. This is equivalent to the Fisher scoring method in which the Hessian matrix used in the Newton–Raphson method is replaced by its expected value. In the case of canonical links the Fisher scoring method and the Newton–Raphson method are identical. Starting values for the iterative procedure are obtained by replacing the  by  in the appropriate equations.

Let  be the sum of squares of deviations from the mean, , for the variable  for a sample of size , i.e.,

and let  be the cross-products of deviations from the means,  and , for the variables  and  for a sample of size , i.e.,

Then the sample covariance of  and  is

and the product-moment correlation coefficient is

g02buc (nag_sum_sqs) computes the sample sums of squares and cross-products deviations from the means (optionally weighted). g02btc (nag_sum_sqs_update) updates the sample sums of squares and cross-products and deviations from the means by the addition/deletion of a (weighted) observation. g02bwc (nag_cov_to_corr) computes the product-moment correlation coefficients from the sample sums of squares and cross-products of deviations from the means. The three routines compute only the upper triangle of the correlation matrix which is stored in a one-dimensional array in packed form. g02bxc (nag_corr_cov) computes both the (optionally weighted) covariance matrix and the (optionally weighted) correlation matrix. These are returned in two-dimensional arrays. (Note that g02btc (nag_sum_sqs_update) and g02buc (nag_sum_sqs) can be used to compute the sums of squares from zero.)

g02brc (nag_ken_spe_corr_coeff) computes Kendall and/or Spearman non-parametric rank correlation coefficients. The function allows for a subset of variables to be selected and for observations to be excluded from the calculations if, for example, they contain missing values.

g02byc (nag_partial_corr) computes a matrix of partial correlation coefficients from the correlation coefficients or variance-covariance matrix returned by g02bxc (nag_corr_cov).

g02hlc (nag_robust_m_corr_user_fn) and g02hmc (nag_robust_m_corr_user_fn_no_derr) compute robust esimates of the variance-covariance matrix by solving the equations

and

as described in Section [Partial Correlation] for user-supplied functions  and . Two options are available for , either  for all  or .

g02hmc (nag_robust_m_corr_user_fn_no_derr) requires only the function  and  to be supplied while g02hlc (nag_robust_m_corr_user_fn) also requires their derivatives.

In general g02hlc (nag_robust_m_corr_user_fn) will be considerably faster than g02hmc (nag_robust_m_corr_user_fn_no_derr) and should be used if derivatives are available.

g02hkc (nag_robust_corr_estim) computes a robust variance-covariance matrix for the following functions:

and

for constants ,  and .

These functions solve a minimax space problem considered by Huber (1981). The values of ,  and  are calculated from the fraction of gross errors; see Hampel et al. (1986) and Huber (1981).

To compute a correlation matrix from the variance-covariance matrix g02bwc (nag_cov_to_corr) may be used.

Two functions are provided for simple linear regression. The function g02cac (nag_simple_linear_regression) calculates the parameter estimates for a simple linear regression with or without a constant term. The function g02cbc (nag_regress_confid_interval) calculates fitted values, residuals and confidence intervals for both the fitted line and individual observations. This function produces the information required for various regression plots.

Note: g02dec (nag_regsn_mult_linear_add_var) also allows the user to initialise a model building process and then to build up the model by adding variables one at a time.

To aid the selection of a regression model the following routines are available.

For robust regression using -estimates instead of least-squares the routine g02hac (nag_robust_m_regsn_estim) will generally be suitable. g02hac (nag_robust_m_regsn_estim) provides a choice of four -functions (Huber's, Hampel's, Andrew's and Tukey's) plus two different weighting methods and the option not to use weights. If other weights or different -functions are needed the routine g02hdc (nag_robust_m_regsn_user_fn) may be used. g02hdc (nag_robust_m_regsn_user_fn) requires the user to supply weights, if required, and also routines to calculate the -function and, optionally, the -function. g02hbc (nag_robust_m_regsn_wts) can be used in calculating suitable weights. The routine g02hfc (nag_robust_m_regsn_param_var) can be used after a call to g02hdc (nag_robust_m_regsn_user_fn) in order to calculate the variance-covariance estimate of the estimated regression coefficients.

For robust regression, using least absolute deviation, e02gac (nag_lone_fit) can be used.

There are four routines for fitting generalized linear models. The output includes: the deviance, parameter estimates and their standard errors, fitted values, residuals and leverages. The routines are:

While g02gac (nag_glm_normal) can be used to fit linear regression models (i.e., by using an identity link) this is not recomended as g02dac (nag_regsn_mult_linear) will fit these models more efficiently. g02gcc (nag_glm_poisson) can be used to fit log-linear models to contingency tables.

In addition to the routines to fit the models there are two routines to aid the interpretation of the model if a model which is not of full rank has been fitted, i.e., aliasing is present.

No routines are currently provided in this chapter for polynomial regression. Users wishing to perform polynomial regressions do however have three alternatives: they can use the multiple linear regression routines, g02dac (nag_regsn_mult_linear), with a set of independent variables which are in fact simply the same single variable raised to different powers, or they can use the routine g04eac (nag_dummy_vars) to compute orthogonal polynomials which can then be used with g02dac (nag_regsn_mult_linear), or they can use the routines in Chapter e02 (Curve and Surface Fitting) which fit polynomials to sets of data points using the techniques of orthogonal polynomials. This latter course is to be preferred, since it is more efficient and liable to be more accurate, but in some cases more statistical information may be required than is provided by those routines, and it may be necessary to use the routines of this chapter.

More general nonlinear regression models may be fitted using the optimization routines in Chapter e04, which contains routines to minimize the function

where the regression parameters are the variables of the minimization problem.