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 $-1.0$ to $+1.0$. A coefficient of zero always indicates that no linear relationship exists; a $+1.0$ 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 $-1.0$ 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 $+1.0$ and $-1.0$ 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 nonparametric coefficients would be $+1.0$ (in (c)) and $-1.0$ (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, you are 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, $r$, based on $n$ pairs of observations, testing against the null hypothesis that there is no correlation between the two variables, the statistic has a Student's $t$-distribution with $n-2$ degrees of freedom; its significance can be tested accordingly.

Ranked and ordinal scale data are generally analysed by nonparametric 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 nonparametric methods, but such techniques are statistically less powerful than a product-moment method. For a Spearman rank order correlation coefficient, $R$, based on $n$ 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 $t$-distribution with $n-2$ degrees of freedom, and may be treated accordingly. (This is similar to the product-moment correlation coefficient, $r$, see above.) Kendall's tau coefficient, based on $n$ 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, $\mathrm{N}\left(0,1\right)$.

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 you.

The correlation coefficients described above measure the association between two variables ignoring any other variables in the system. Suppose there are three variables $X,Y$ and $Z$ as shown in the path diagram below. The association between $Y$ and $Z$ is made up of the direct association between $Y$ and $Z$ and the association caused by the path through $X$, that is the association of both $Y$ and $Z$ with the third variable $X$. For example if $Z$ and $Y$ were cholesterol level and blood pressure and $X$ 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 known as the partial correlation. If ${\rho }_{zy}$, ${\rho }_{zx}$ and ${\rho }_{xy}$ represent the correlations between $x$, $y$ and $z$ then the partial correlation between $Z$ and $Y$ given $X$ 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 $Y$ and $X$ with $n_y$ variables in $Y$ and $n_x$ variables in $X$ and let the variance-covariance matrix of all $n_y+n_x$ variables be partitioned into Then the variance-covariance of $Y$ conditional on fixed values of the $X$ variables is given by The partial correlation matrix is then computed by standardizing ${\Sigma }_{y\mid x}$.

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, $C$, can be written as where ${\tau }^2$ is a correction factor to give an unbiased estimator if the data is Normal and $A$ is a lower triangular matrix. Let $x_i$ be the vector of values for the $i$th observation and let $z_i=A\left(x_i-\theta \right)$, $\theta$ being a robust estimate of location, then $\theta$ and $A$ are found as solutions to and where $w\left(t\right)$, $u\left(t\right)$ and $v\left(t\right)$ are functions such that they return a value of $1$ for reasonable values of $t$ and decreasing values for large $t$. The correlation matrix can then be calculated from the variance-covariance matrix. If $w$, $u$, and $v$ returned $1$ for all values then the product-moment correlation coefficient would be calculated.

When there are missing values in the data these may be handled in one of two ways. Firstly, if a case contains a missing observation for any variable, then that case is omitted in its entirety from all calculations; this may be termed casewise treatment of missing data. Secondly, if a case contains a missing observation for any variable, then the case is omitted from only those calculations involving the variable for which the value is missing; this may be called pairwise treatment of missing data. Pairwise deletion of missing data has the advantage of using as much of the data as possible in the computation of each coefficient. In extreme circumstances, however, it can have the disadvantage of producing coefficients which are based on a different number of cases, and even on different selections of cases or samples; furthermore, the ‘correlation’ matrices formed in this way need not necessarily be positive semidefinite, a requirement for a correlation matrix. Casewise deletion of missing data generally causes fewer cases to be used in the calculation of the coefficients than does pairwise deletion. How great this difference is will obviously depend on the distribution of the missing data, both among cases and among variables.

Pairwise treatment does therefore use more information from the sample, but should not be used without careful consideration of the location of the missing observations in the data matrix, and the consequent effect of processing the missing data in that fashion.

A correlation matrix is, by definition, a symmetric, positive semidefinite matrix with unit diagonals and all elements in the range $\left[-1,1\right]$.

In practice, rather than having a true correlation matrix, you may find that you have a matrix of pairwise correlations. This usually occurs in the presence of missing values, when the missing values are treated in a pairwise fashion as discussed in [Missing values]. Matrices constructed in this way may not be not positive semidefinite, and therefore are not a valid correlation matrix. However, a valid correlation matrix can be calculated that is in some sense ‘close’ to the original.

Given an $n\times n$ matrix, $G$, there are a number of available ways of computing the ‘nearest’ correlation matrix, $\Sigma$ to $G$:

(a)

Frobenius Norm Find $\Sigma$ such that $$\sum _{i=1}^n\sum _{j=1}^n{\left(s_{ij}-{\sigma }_{ij}\right)}^2$$ is minimized. Where $S$ is the symmetric matrix defined as $S=\frac{1}{2}\left(G+G^{\mathrm{T}}\right)$ and $s_{ij}$ and ${\sigma }_{ij}$ denotes the elements of $S$ and $\Sigma$ respectively. A weighted Frobenius norm can also be used. The term being summed across therefore becomes $w_i{w}_j{\left(s_{ij}-{\sigma }_{ij}\right)}^2$ if row and column weights are being used or $w_{ij}{\left(s_{ij}-{\sigma }_{ij}\right)}^2$ when element-wise weights are used.

(b)

Factor Loading Method This method is similar to (a) in that it finds a $\Sigma$ that is closest to $S$ in the Frobenius norm. However, it also ensures that $\Sigma$ has a $k$-factor structure, that is $\Sigma$ can be written as $$\Sigma =X{X}^{\mathrm{T}}+\mathrm{diag}\left(I-X{X}^{\mathrm{T}}\right)$$ where $I$ is the identity matrix and $X$ has $n$ rows and $k$ columns. $X$ is often referred to as the factor loading matrix. This problem primarily arises when a factor model $\xi =X\eta +D\epsilon$ is used to describe a multivariate time series or collateralized debt obligations. In this model $\eta \in {ℝ}^k$ and $\xi \in {ℝ}^n$ are vectors of independent random variables having zero mean and unit variance, with $\eta$ and $\epsilon$ independent of each other, and $X\in {ℝ}^{n\times k}$ with $D\in {ℝ}^{n\times n}$ diagonal. In the case of modelling debt obligations $\xi$ can, for example, model the equity returns of $n$ different companies of a portfolio where $\eta$ describes $k$ factors influencing all companies, in contrast to the elements of $\epsilon$ having only an effect on the equity of the corresponding company. With this model the complex behaviour of a portfolio, with potentially thousands of equities, is captured by looking at the major factors driving the behaviour. The number of factors usually chosen is a lot smaller than $n$, perhaps between $1$ and $10$, yielding a large reduction in the complexity. The number of the factors, $k$, which yields a matrix $X$ such that ${‖G-X{X}^{\mathrm{T}}+\mathrm{diag}\left(I-X{X}^{\mathrm{T}}\right)‖}_F$ is within a required tolerance can also be determined, by experimenting with the input $k$ and comparing the norms.

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:

(i)	to obtain best estimates of the unknown regression parameters;
(ii)	to test hypotheses about these parameters;
(iii)	to determine the adequacy of the assumed model; and
(iv)	to verify the set of relevant assumptions.

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 $x_j$, for $\mathit{j}=1,2,3,4$, such that: along with a variable for the mean $x_0$: If there were $7$ observations the data would be:

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 $x_4$ cannot be distinguished from the combined effect of $x_0,x_1,x_2$ and $x_3$. 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 dataset 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 $x_1$ term would also be aliased. In general intrinsic aliasing may be overcome by changing the model, e.g., remove $x_0$ or $x_1$ from the model, or by introducing constraints on the parameters, e.g., ${\beta }_1+{\beta }_2+{\beta }_3+{\beta }_4=0$.

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 (${\beta }_0+{\beta }_1$) and (${\beta }_2-{\beta }_3$) 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.

(a)

All possible regressions In this case all the possible combinations of independent variables are fitted and the one considered the best selected. To choose the best, two conflicting criteria have to be balanced. One is the fit of the model which will improve as more variables are added to the model. The second criterion is the desire to have a model with a small number of significant terms. Depending on how the model is fit, statistics such as $R^2$, which gives the proportion of variation explained by the model, and $C_p$, which tries to balance the size of the residual sum of squares against the number of terms in the model, can be used to aid in the choice of model.

(b)

Stepwise model building In stepwise model building the regression model is constructed recursively, adding or deleting the independent variables one at a time. When the model is built up the procedure is known as forward selection. The first step is to choose the single variable which is the best predictor. The second independent variable to be added to the regression equation is that which provides the best fit in conjunction with the first variable. Further variables are then added in this recursive fashion, adding at each step the optimum variable, given the other variables already in the equation. Alternatively, backward elimination can be used. This is when all variables are added and then the variables dropped one at a time, the variable dropped being the one which has the least effect on the fit of the model at that stage. There are also hybrid techniques which combine forward selection with backward elimination.

One 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 $x_0=1$ for all observations, is to be fitted to the $n$ data points such that where $e_i$ are unknown independent random errors with $E\left(e_i\right)=0$ and $\mathrm{var}\left(e_i\right)={\sigma }^2$, ${\sigma }^2$ being a constant, then the method used is to calculate the estimates of the regression parameters ${\beta }_0,{\beta }_1,{\beta }_2,\dots ,{\beta }_p$ by minimizing 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).

Let $X$ be the $n$ by $p$ matrix of independent variables and $y$ be the vector of values for the dependent variable. To find the least squares estimates of the vector of parameters, $\hat{\beta }$, the $QR$ decomposition of $X$ is found, i.e., where $R^{*}=\left(\begin{array}{c}R\\ 0\end{array}\right)$, $R$ being a $p$ by $p$ upper triangular matrix, and $Q$ an $n$ by $n$ orthogonal matrix. If $R$ is of full rank then $\hat{\beta }$ is the solution to where $c=Q^{\mathrm{T}}y$ and $c_1$ is the first $p$ rows of $c$. If $R$ is not of full rank, a solution is obtained by means of a singular value decomposition (SVD) of $R$, where $D$ is a $k$ by $k$ diagonal matrix with nonzero diagonal elements, $k$ being the rank of $R$, and $Q_{*}$ and $P$ are $p$ by $p$ orthogonal matrices. This gives the solution

$P_1$ being the first $k$ columns of $P$ and $Q_{{*}_1}$ being the first $k$ columns of $Q_{*}$.

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 $w$ is required then both $X$ and $y$ are premultiplied by $w^{1/2}$.

The method described above will, in general, be more accurate than methods based on forming ($X^{\mathrm{T}}X$), (or a scaled version), and then solving the equations

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 $t$-tests can be calculated to test hypotheses about the regression parameters – for example, whether the value of the parameter, ${\beta }_k$, is significantly different from a specified value, $b_k$ (often zero). If the estimate of ${\beta }_k$ is ${\hat{\beta }}_k$ and its standard error is $\mathrm{se}\left({\hat{\beta }}_k\right)$ then the $t$-statistic is It should be noted that both the tests and the confidence intervals may not be independent. Alternatively $F$-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 ${\mathit{RSS}}_1$ with degrees of freedom ${\nu }_1$, is a sub-model of model $2$, giving residual sum of squares ${\mathit{RSS}}_2$ with degrees of freedom ${\nu }_2$, 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 $F$-statistic may be used. These tests and confidence intervals require the additional assumption that the errors, $e_i$, are Normally distributed.

To check for systematic lack of fit the residuals, $r_i=y_i-{\hat{y}}_i$, where ${\hat{y}}_i$ is the fitted value, should be examined. If the model is correct then they should be random with no discernible 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, $y$, $\hat{y}=Hy$. The variance-covariance matrix of the residuals is then $\left(I-H\right){\sigma }^2$, $I$ being the identity matrix. The diagonal elements of $H$, $h_{ii}$, can therefore be used to standardize the residuals. The $h_{ii}$ are a measure of the effect of the $i$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 $e_i$ and hence the independence of the observations.

When data on predictor variables $x$ are multicollinear, ridge regression models provide an alternative to variable selection in the multiple regression model. In the ridge regression case, parameter estimates in the linear model are found by penalised least squares: where the value of the ridge parameter $h$ controls the trade-off between the goodness-of-fit and smoothness of a solution.

Consider the piecewise linear loss function

The minimum of the expectation can be obtained by using the integral rule of Leibnitz to differentiate with respect to $z$ and then setting the result to zero, giving hence $\alpha =F^{-1}\left(\tau \right)$ when the solution is unique. If the solution is not unique then there exists a range of quantiles, each of which is equally valid. Taking the smallest value of such a range ensures that the empirical quantile function is left-continuous. Therefore obtaining the $\tau$th quantile of a distribution $F$ can be achieved by minimizing the expected value of the loss function ${\rho }_{\tau }$.

This idea of obtaining the quantile by solving an optimization problem can be extended to finding the $\tau$th sample quantile. Given a vector of $n$ observed values, $y$, from some distribution the empirical distribution function, $F_n\left(\alpha \right)=n^{-1}{\displaystyle \sum _{\mathit{i}=1}^n}I\left(y_i\le \alpha \right)$ provides an estimate of the unknown distribution function $F$ giving an expected loss of and therefore the problem of finding the $\tau$th sample quantile, $\hat{\alpha }\left(\tau \right)$, can be expressed as finding the solution to the problem effectively replacing the operation of sorting, usually required when obtaining a sample quantile, with an optimization.

Given the vector $y$ it is a well known result that the sample mean, $\hat{y}$, solves the least squares problem

This result leads to least squares regression where, given design matrix $X$ and defining the conditional mean of $y$ as $\mu \left(X\right)=X\beta$, an estimate of $\beta$ is obtained from the solution to

Quantile regression can be derived in a similar manner by specifying the $\tau$th conditional quantile as $Q_y\left(\tau |X\right)=X\beta \left(\tau \right)$ and estimating $\beta \left(\tau \right)$ as the solution to

By introducing $2n$ slack variables, $u=\left\{u_i:i=1,2,\dots ,n\right\}$ and $v=\left\{u_i:i=1,2,\dots ,n\right\}$, the quantile regression minimization problem, (10), can be expressed as a linear programming (LP) problem, with primal and associated dual formulations

(a)	Primal form $$\underset{\left(u,v,\beta \right)\in {ℝ}_{+}^n\times {ℝ}_{+}^n\times {ℝ}^p}{\mathrm{minimize}}\tau e^{\mathrm{T}}u+\left(1-\tau \right)e^{\mathrm{T}}v\text{ \hspace{1em} subject to \hspace{1em}}y=X\beta +u-v$$ (11) where $e$ is a vector of length $n$, where each element is $1$. If $r_i$ denotes the $i$th residual, $r_i=y_i-x_i^{\mathrm{T}}\beta$, then the slack variables, $\left(u,v\right)$, can be thought as corresponding to the absolute value of the positive and negative residuals respectively with $$\begin{array}{ccc}{u}_i=\left\{\begin{array}{l}{r}_i\text{ if }{r}_i>0\\ 0\text{ otherwise}\end{array}\right.& & v_i=\left\{\begin{array}{l}-r_i\text{ if }{r}_i<0\\ 0\text{ otherwise}\end{array}\right.\end{array}$$
(b)	Dual form The dual formulation of (11) is given by $$\underset{d}{\mathrm{maximize}}{y}^{\mathrm{T}}d\text{ \hspace{1em} subject to \hspace{1em}}{X}^{\mathrm{T}}d=0,d\in {\left[\tau -1,\tau \right]}^n$$ which, on setting $a=d+\left(1-\tau \right)e$, is equivalent to $$\underset{a}{\mathrm{maximize}}{y}^{\mathrm{T}}a\text{ \hspace{1em} subject to \hspace{1em}}{X}^{\mathrm{T}}a=\left(1-\tau \right)X^{\mathrm{T}}e,a\in {\left[0,1\right]}^n$$ (12)
(c)	Canonical form Linear programming problems are often described in a standard way, called the canonical form. The canonical form of an LP problem is $$\underset{z}{\mathrm{minimize}}{c}^{\mathrm{T}}z\text{ \hspace{1em} subject to \hspace{1em}}{l}_l\le \left\{\begin{array}{c}z\\ Az\end{array}\right\}\le l_u\text{.}$$ Letting $0_p$ denote a vector of $p$ zeros $\pm {\infty }_p$ denote a vector of $p$ arbitarily small or large values, $I_{n\times n}$ denote the $n\times n$ identity matrix, $c=\left\{a,b\right\}$ denote the row vector constructed by concatenating the elements of vector $b$ to the elements of vector $a$ and $C=\left[A,B\right]$ denote the matrix constructed by concatenating the columns of matrix $B$ onto the columns of matrix $A$ then setting $$\begin{array}{ccc}{c}^{\mathrm{T}}=\left\{0_p,\tau e^{\mathrm{T}},\left(1-\tau \right)e^{\mathrm{T}}\right\}& & z^{\mathrm{T}}=\left\{{\beta }^{\mathrm{T}},u^{\mathrm{T}},v^{\mathrm{T}}\right\}\\ A=\left[X,I_{n\times n},-I_{n\times n}\right]& & b=y\\ l_u=\left\{{+\infty }_p,{\infty }_n,{\infty }_n,y\right\}& & l_l=\left\{{-\infty }_p,0_n,0_n,y\right\}\end{array}$$ gives the quantile regression LP problem as described in (11).

Once expressed as an LP problem the parameter estimates $\hat{\beta }\left(\tau \right)$ can be obtained in a number of ways, for example via the inertia-controlling method of Gill and Murray (1978) (see e04mf), the simplex method or an interior point method as used by g02qf and g02qg.

Koenker (2005) shows that the limiting covariance matrix of $\sqrt{n}\left(\hat{\beta }\left(\tau \right)-\beta \left(\tau \right)\right)$ is of the form of a Huber Sandwich. Therefore, under the assumption of Normally distributed errors where and $f_i\left(Q_{y_i}\left(\tau |x_i\right)\right)$ denotes the conditional density of the response $y$ evaluated at the $\tau$th conditional quantile.

More generally, the asymptotic covariance matrix for $\hat{\beta }\left({\tau }_1\right),\hat{\beta }\left({\tau }_1\right),\dots ,\hat{\beta }\left({\tau }_n\right)$ has blocks defined by

Under the assumption of independent, identically distributed (iid) errors, (13) simplifies to where $s\left(\tau \right)$ is the sparsity function, given by a similar simplification occurs with (14).

In cases where the assumption of iid errors does not hold, Powell (1991) suggests using a kernel estimator of the form for some bandwidth parameter $c_n$ satisfying ${\displaystyle \underset{n\to \infty }{\lim }}{c}_n\to 0$ and ${\displaystyle \underset{n\to \infty }{\lim }}\sqrt{n}{c}_n\to \infty$ and Hendricks and Koenker (1991) suggest a method based on an extension of the idea of sparsity.

Rather than use an asymptotic estimate of the covariance matrix, it is also possible to use bootstrapping. Roughly speaking the original data is resampled and a set of parameter estimates obtained from each new sample. A sample covariance matrix is then constructed from the resulting matrix of parameter estimates.

Let ${\mathit{SS}}_x$ be the sum of squares of deviations from the mean, $\stackrel{-}{x}$, for the variable $x$ for a sample of size $n$, i.e., and let ${\mathit{SC}}_{xy}$ be the cross-products of deviations from the means, $\stackrel{-}{x}$ and $\stackrel{-}{y}$, for the variables $x$ and $y$ for a sample of size $n$, i.e., Then the sample covariance of $x$ and $y$ is and the product-moment correlation coefficient is

g02ba computes the product-moment correlation coefficients.

g02bt updates the sample sums of squares and cross-products and deviations from the means by the addition/deletion of a (weighted) observation.

g02bu computes the sample sums of squares and cross-products deviations from the means (optionally weighted). The output from multiple calls to g02bu can be combined via a call to (G02BZF not in this release), allowing large datasets to be summarised across multiple processing units.

g02bt updates the sample sums of squares and cross-products and deviations from the means by the addition/deletion of a (weighted) observation.

g02bw computes the product-moment correlation coefficients from the sample sums of squares and cross-products of deviations from the means.

The three methods compute only the upper triangle of the correlation matrix which is stored in a one-dimensional array in packed form.

g02bx computes both the (optionally weighted) covariance matrix and the (optionally weighted) correlation matrix. These are returned in two-dimensional arrays. (Note that g02bt and g02bu can be used to compute the sums of squares from zero.)

g02bg can be used to calculate the correlation coefficients for a subset of variables in the data matrix.

If there are missing values then g02bu and g02bx, as described above, will allow casewise deletion by you giving the observation zero weight (compared with unit weight for an otherwise unweighted computation).

Other methods also handle missing values in the calculation of unweighted product-moment correlation coefficients. Casewise exclusion of missing values is provided by g02bb while pairwise omission of missing values is carried out by g02bc. These two methods calculate a correlation matrix for all the variables in the data matrix; similar output but for only a selected subset of variables is provided by methods g02bh and g02bj respectively. As well as providing the Pearson product-moment correlation coefficients, these methods also calculate the means and standard deviations of the variables, and the matrix of sums of squares and cross-products of deviations from the means. For all four methods you are free to select appropriate values for consideration as missing values, bearing in mind the nature of the data and the possible range of valid values. The missing values for each variable may be either different or alike and it is not necessary to specify missing values for all the variables.

There are five methods which perform nonparametric correlations, each of which is capable of producing both Spearman's rank order and Kendall's tau correlation coefficients. The basic underlying concept of both these methods is to replace each observation by its corresponding rank or order within the observations on that variable, and the correlations are then calculated using these ranks.

It is obviously more convenient to order the observations and calculate the ranks for a particular variable just once, and to store these ranks for subsequent use in calculating all coefficients involving that variable; this does however require an amount of store of the same size as the original data matrix, which in some cases might be excessive. Accordingly, some methods calculate the ranks only once, and replace the input data matrix by the matrix of ranks, which are then also made available to you on exit from the method, while others preserve the data matrix and calculate the ranks a number of times within the method; the ranks of the observations are not provided as output by methods which work in the latter way. If it is possible to arrange the program in such a way that the first technique can be used, then efficiency of timing is achieved with no additional storage, whereas in the second case, it is necessary to have a second matrix of the same size as the data matrix, which may not be acceptable in certain circumstances; in this case it is necessary to reach a compromise between efficiency of time and of storage, and this may well be dependent upon local conditions.

Methods g02bn and g02bq both calculate Kendall's tau and/or Spearman's rank order correlation coefficients taking no account of missing values; g02bn does so by calculating the ranks of each variable only once, and replacing the data matrix by the matrix of ranks, whereas g02bq calculates the ranks of each variable several times. Methods g02bp and g02br provide the same output, but treat missing values in a ‘casewise’ manner (see above); g02bp calculates the ranks of each variable only once, and overwrites the data matrix of ranks, while g02br determines the ranks of each variable several times. For ‘pairwise’ omission of missing data (see above), the method g02bs provides Kendall and/or Spearman coefficients.

Since g02bn and g02bp order the observations and calculate the ranks of each variable only once, then if there are $M$ variables involved, there are only $M$ separate ‘ranking’ operations; this should be contrasted with the method used by methods g02bq and g02br which perform $M\left(M-1\right)/2+1$ similar ranking operations. These ranking operations are by far the most time-consuming parts of these nonparametric methods, so for a matrix of as few as five variables, the time taken by one of the slower methods can be expected to be at least a factor of two slower than the corresponding efficient method; as the number of variables increases, so this relative efficiency factor increases. Only one method, g02bs, is provided for pairwise missing values, and this method carries out $M\left(M-1\right)$ separate rankings; since by the very nature of the pairwise method it is necessary to treat each pair of variables separately and rank them individually, it is impossible to reduce this number of operations, and so no alternative method is provided.

g02by computes a matrix of partial correlation coefficients from the correlation coefficients or variance-covariance matrix returned by g02bx.

g02hl and g02hm compute robust estimates of the variance-covariance matrix by solving the equations and as described in [Robust estimation of correlation coefficients] for user-supplied functions $w$ and $u$. Two options are available for $v$, either $v\left(t\right)=1$ for all $t$ or $v\left(t\right)=u\left(t\right)$.

g02hm requires only the function $w$ and $u$ to be supplied while g02hl also requires their derivatives.

In general g02hl will be considerably faster than g02hm and should be used if derivatives are available.

g02hk computes a robust variance-covariance matrix for the following functions: and for constants $a_u$, $b_u$ and $c_w$.

These functions solve a minimax space problem considered by Huber (1981). The values of $a_u$, $b_u$ and $c_w$ 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 g02bw may be used.

Four methods are provided to calculate a nearest correlation matrix. The choice of method will depend on what definition of ‘nearest’ is required and whether there is any particular structure desired in the resulting correlation matrix.

g02aa computes the nearest correlation matrix in the Frobenius norm, using the method of Qi and Sun (2006).

g02ab uses an extension of the method implemented in g02aa allowing for the row and column weighted Frobenius norm to be used as well as bounds on the eigenvalues of the resulting correlation matrix to be specified.

g02ae computes the factor loading matrix, allowing a correlation matrix with a $k$-factor structure to be computed.

(G02AJF not in this release) again computes the nearest correlation matrix in the Frobenius norm, but allows for element-wise weighting as well as bounds on the eigenvalues.

Four methods are provided for simple linear regressions: g02ca and g02cc perform the simple linear regression with a constant term (equation (1) above), while g02cb and g02cd fit the simple linear regression with no constant term (equation (2) above). Two of these methods, g02cc and g02cd, take account of missing values, which the others do not. In these two methods, an observation is omitted if it contains a missing value for either the dependent or the independent variable; this is equivalent to both the casewise and pairwise methods, since both are identical when there are only two variables involved. Input to these methods consists of the raw data, and output includes the coefficients, their standard errors and $t$ values for testing the significance of the coefficients; the $F$ value for testing the overall significance of the regression is also given.

g02ka calculates a ridge regression, optimizing the ridge parameter according to one of four prediction error criteria.

g02kb calculates ridge regressions for a given set of ridge parameters.

No methods are currently provided in this chapter for polynomial regression. If you wish to perform polynomial regressions you have three alternatives: you can use the multiple linear regression methods, g02da, with a set of independent variables which are in fact simply the same single variable raised to different powers, or you can use the method (G04EAF not in this release) to compute orthogonal polynomials which can then be used with g02da, or you can use the methods in E02 class (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 methods, and it may be necessary to use the methods of this chapter.

More general nonlinear regression models may be fitted using the optimization methods in E04 class, which contains methods to minimize the function where the regression parameters are the variables of the minimization problem.

g02da fits a general linear regression model using the $QR$ method and an SVD if the model is not of full rank. The results returned include: residual sum of squares, parameter estimates, their standard errors and variance-covariance matrix, residuals and leverages. There are also several methods to modify the model fitted by g02da and to aid in the interpretation of the model.

g02dc adds or deletes an observation from the model.

g02dd computes the parameter estimates, and their standard errors and variance-covariance matrix for a model that is modified by g02dc, g02de or g02df.

g02de adds a new variable to a model.

g02df drops a variable from a model.

g02dg fits the regression to a new dependent variable, i.e., keeping the same independent variables.

g02dk calculates the estimates of the parameters for a given set of constraints, (e.g., parameters for the levels of a factor sum to zero) for a model which is not of full rank and the SVD has been used.

g02dn calculates the estimate of an estimable function and its standard error.

Note:  g02de also allows you to initialize a model building process and then to build up the model by adding variables one at a time.

If you wish to use methods based on forming the cross-products/correlation matrix (i.e., ($X^{\mathrm{T}}X$) matrix) rather than the recommended use of g02da then the following methods should be used.

For regression through the origin (i.e., no constant) g02ch preceded by:

• g02bd (no missing values, all variables)
• g02bk (no missing values, subset of variables)
• g02be (casewise missing values, all variables)
• g02bl(casewise missing values, subset of variables)
• g02bf* (pairwise missing values, all variables)
• g02bm* (pairwise missing values, subset of variables)

For regression with intercept (i.e., with constant) g02cg preceded by:

• g02ba (no missing values, all variables)
• g02bg (no missing values, subset of variables)
• g02bb (casewise missing values, all variables)
• g02bh (casewise missing values, subset of variables)
• g02bc* (pairwise missing values, all variables)
• g02bj* (pairwise missing values, subset of variables)

Note that the four methods using pairwise deletion of missing value (marked with $*$) should be used with great caution as the use of this method can lead to misleading results, particularly if a significant proportion of values are missing.

Both g02cg and g02ch require that the correlations/sums of squares involving the dependent variable must appear as the last row/column. Because the layout of the variables in your data array may not be arranged in this way, two methods, g02ce and g02cf, are provided for rearranging the rows and columns of vectors and matrices. g02cf simply reorders the rows and columns while g02ce forms smaller vectors and matrices from larger ones.

Output from g02cg and g02ch consists of the coefficients, their standard errors, $R^2$-values, $t$ and $F$ statistics.

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

(G02EAF not in this release) computes the residual sums of squares for all possible regressions for a given set of dependent variables. The method allows some variables to be forced into all regressions.

(G02ECF not in this release) computes the values of $R^2$ and $C_p$ from the residual sums of squares as provided by (G02EAF not in this release).

(G02EEF not in this release) enables you to fit a model by forward selection. You may call (G02EEF not in this release) a number of times. At each call the method will calculate the changes in the residual sum of squares from adding each of the variables not already included in the model, select the variable which gives the largest change and then if the change in residual sum of squares meets the given criterion will add it to the model.

g02ef uses a full stepwise selection to choose a subset of the explanatory variables. The method repeatedly applies a forward selection step followed by a backward elimination step until neither step updates the current model.

g02fa computes the following standardized residuals and measures of influence for the residuals and leverages produced by g02da:

(i)	Internally studentized residual;
(ii)	Externally studentized residual;
(iii)	Cook's $D$ statistic;
(iv)	Atkinson's $T$ statistic.

g02fc computes the Durbin–Watson test statistic and bounds for its significance to test for serial correlation in the errors, $e_i$.

For robust regression using $M$-estimates instead of least squares the method g02ha will generally be suitable. g02ha provides a choice of four $\psi$-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 $\psi$-functions are needed the method g02hd may be used. g02hd requires you to supply weights, if required, and also methods to calculate the $\psi$-function and, optionally, the $\chi$-function. g02hb can be used in calculating suitable weights. The method g02hf can be used after a call to g02hd in order to calculate the variance-covariance estimate of the estimated regression coefficients.

For robust regression, using least absolute deviation, (E02GAF not in this release) can be used.

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

g02ga Normal distribution.

g02gb binomial distribution.

g02gc Poisson distribution.

g02gd gamma distribution.

While g02ga can be used to fit linear regression models (i.e., by using an identity link) this is not recommended as g02da will fit these models more efficiently. g02gc can be used to fit log-linear models to contingency tables.

In addition to the methods to fit the models there is one method to predict from the fitted model and two methods to aid interpretation when the fitted model is not of full rank, i.e., aliasing is present.

g02gp computes a predicted value and its associated standard error based on a previously fitted generalized linear model.

g02gk computes parameter estimates for a set of constraints, (e.g., sum of effects for a factor is zero), from the SVD solution provided by the fitting method.

g02gn calculates an estimate of an estimable function along with its standard error.

There are four methods for fitting linear mixed effects regression.

g02ja (G02JDF not in this release) uses restricted maximum likelihood (REML) to fit the model.

g02jb (G02JEF not in this release) uses maximum likelihood to fit the model.

For all methods the output includes: either the maximum likelihood or restricted maximum likelihood and the fixed and random parameter estimates, along with their standard errors. Whilst it is possible to fit a hierachical model using g02ja or g02jb, (G02JDF not in this release) (G02JEF not in this release) allow the model to be specified in a more intuitive way. (G02JCF not in this release) must be called prior to calling (G02JDF not in this release) (G02JEF not in this release).

As the estimates of the variance components are found using an iterative procedure initial values must be supplied for each $\sigma$. In all four methods you can either specify these initial values, or allow the method to calculate them from the data using minimum variance quadratic unbiased estimation (MIVQUE0). Setting the maximum number of iterations to zero in any of the methods will return the corresponding likelihood, parameter estimates and standard errors based on these initial values.

Two methods are provided for performing linear quantile regression, g02qf and g02qg. Of these, g02qf provides a simplified interface to g02qg, where many of the input parameters have been given default values and the amount of output available has been reduced.

Prior to calling g02qg the optional parameter array must be initialized by calling (G02ZKF not in this release) with optstr set to initialize. Once these arrays have been initialized (G02ZLF not in this release) can be called to query the value of an optional parameter.

g02la calculates a nonlinear, iterative PLS by using singular value decomposition.

g02lb calculates a nonlinear, iterative PLS by using Wold's method.

g02lc calculates parameter estimates for a given number of PLS factors.

g02ld calculates predictions given a PLS model.