diff --git a/doc/src/Regression/Regression.do.txt b/doc/src/Regression/Regression.do.txt index f57d66c11..216fda678 100644 --- a/doc/src/Regression/Regression.do.txt +++ b/doc/src/Regression/Regression.do.txt @@ -18,6 +18,39 @@ A regression model aims at finding a likelihood function $p(y\vert \hat{x})$, th !eblock + +Consider an experiment in which $p$ characteristics of $n$ samples are +measured. The data from this experiment are denoted $\mathbf{X}$, with +$\mathbf{X}$ as above. The matrix $\mathbf{X}$ is called the *design +matrix*. Additional information of the samples is available in the +form of $\mathbf{Y}$ (also as above). The variable $\mathbf{Y}$ is +generally referred to as the *response variable*. The aim of +regression analysis is to explain $\mathbf{Y}$ in terms of +$\mathbf{X}$ through a functional relationship like $Y_i = +f(\mathbf{X}_{i,\ast})$. When no prior knowledge on the form of +$f(\cdot)$ is available, it is common to assume a linear relationship +between $\mathbf{X}$ and $\mathbf{Y}$. This assumption gives rise to +the *linear regression model*: + +In model (ref{form.linRegressionModel}) $\beta = (\beta_1, \ldots, +\beta_p)^{\top}$ is the *regression parameter*. The parameter +$\beta_j$, $j=1, \ldots, p$, represents the effect size of covariate +$j$ on the response. That is, for each unit change in covariate $j$ +(while keeping the other covariates fixed) the observed change in the +response is equal to $\beta_j$. The second summand on the right-hand +side of the model, $\varepsilon_i$, is referred to as the error. It +represents the part of the response not explained by the functional +part $\mathbf{X}_{i,\ast} \, \beta$ of the model +(ref{form.linRegressionModel}). In contrast to the functional part, +which is considered to be systematic (i.e. non-random), the error is +assumed to be random. Consequently, $Y_{i_1,\ast}$ need not equal +$Y_{i_2,\ast}$ for $i_1 \not= i_2$, even if $\mathbf{X}_{i_1,\ast}= +\mathbf{X}_{i_2,\ast}$. To complete the formulation of model +(ref{form.linRegressionModel}) we need to specify the probability +distribution of $\varepsilon_i$. + + + !split ===== General linear models ===== !bblock @@ -1090,6 +1123,8 @@ probability plot. + + !split ===== The singular value decompostion ===== !bblock @@ -1109,9 +1144,6 @@ We have then !eblock -!split -===== Ridge regression ===== - !split ===== Lasso regression ===== @@ -1208,67 +1240,22 @@ plt.show() -Consider an experiment in which $p$ characteristics of $n$ samples are -measured. The data from this experiment are denoted $\mathbf{X}$, with -$\mathbf{X}$ as above. The matrix $\mathbf{X}$ is called the *design -matrix*. Additional information of the samples is available in the -form of $\mathbf{Y}$ (also as above). The variable $\mathbf{Y}$ is -generally referred to as the *response variable*. The aim of -regression analysis is to explain $\mathbf{Y}$ in terms of -$\mathbf{X}$ through a functional relationship like $Y_i = -f(\mathbf{X}_{i,\ast})$. When no prior knowledge on the form of -$f(\cdot)$ is available, it is common to assume a linear relationship -between $\mathbf{X}$ and $\mathbf{Y}$. This assumption gives rise to -the *linear regression model*: -!bt -\begin{align} -Y_{i} & = \mathbf{X}_{i,\ast} \, \beta + \varepsilon_i label{form.linRegressionModel} -\\ -\nonumber -& = \beta_1 \, X_{i,1} + \ldots + \beta_{p} \, X_{i, p} + \varepsilon_i. -\end{align} -!et - -In model (ref{form.linRegressionModel}) $\beta = (\beta_1, \ldots, -\beta_p)^{\top}$ is the *regression parameter*. The parameter -$\beta_j$, $j=1, \ldots, p$, represents the effect size of covariate -$j$ on the response. That is, for each unit change in covariate $j$ -(while keeping the other covariates fixed) the observed change in the -response is equal to $\beta_j$. The second summand on the right-hand -side of the model, $\varepsilon_i$, is referred to as the error. It -represents the part of the response not explained by the functional -part $\mathbf{X}_{i,\ast} \, \beta$ of the model -(ref{form.linRegressionModel}). In contrast to the functional part, -which is considered to be systematic (i.e. non-random), the error is -assumed to be random. Consequently, $Y_{i_1,\ast}$ need not equal -$Y_{i_2,\ast}$ for $i_1 \not= i_2$, even if $\mathbf{X}_{i_1,\ast}= -\mathbf{X}_{i_2,\ast}$. To complete the formulation of model -(ref{form.linRegressionModel}) we need to specify the probability -distribution of $\varepsilon_i$. - -======= Ridge regression ======= - -When the design matrix is high-dimensional, the covariates (the -columns of $\mathbf{X}$) are super-collinear. Recall *collinearity* in -regression analysis refers to the event of two (or multiple) -covariates being highly linearly related. Consequently, the subspace -spanned by collinear covariates may not be (or close to not being) of -full rank. When the subspace (onto which $\mathbf{Y}$ is projected) -is (close to) rank deficient, it is (almost) impossible to separate -the contribution of the individual covariates. The uncertainty with -respect to the covariate responsible for the variation explained in -$\mathbf{Y}$ is often reflected in the fit of the linear regression -model to data by a large error of the estimates of the regression -parameters corresponding to the collinear covariates. -Consider the design matrix: +!split +===== From standard regression to Ridge regressions ===== +One of the typical problems we encounter with linear regression, in particular +when the matrix $\hat{X}$ (our so-called design matrix) is high-dimensional, +are problems with near singular or singular matrices. The column vectors of $\hat{X}$ +may be linearly dependent, normally referred to as super-collinearity. +This means that the matrix may be rank deficient and it is basically impossible to +to model the data using linear regression. As an example, consider the matrix !bt \begin{align*} -\mathbf{X} & = \left( +\mathbf{X} & = \left[ \begin{array}{rrr} 1 & -1 & 2 \\ @@ -1277,92 +1264,56 @@ Consider the design matrix: 1 & 2 & -1 \\ 1 & 1 & 0 -\end{array} \right) +\end{array} \right] \end{align*} !et -The columns of $\mathbf{X}$ are linearly dependent: the first column -is the row-wise sum of the other two columns. The rank (more correct, -the column rank) of a matrix is the dimension of space spanned by the +The columns of $\hat{X}$ are linearly dependent. We se this easily since the +the first column is the row-wise sum of the other two columns. The rank (more correct, +the column rank) of a matrix is the dimension of the space spanned by the column vectors. Hence, the rank of $\mathbf{X}$ is equal to the number -of linearly independent columns: $\mbox{rank}(\mathbf{X}) = 2$. +of linearly independent columns. In this particular case the matrix has rank 2. - - - -Super-collinearity of an $(n \times p)$-dimensional design matrix $\mathbf{X}$ implies[^footnote1][^footnote1]: If the (column) rank of $\mathbf{X$ is smaller than $p$, there exists a non-trivial $\mathbf{v} \in \mathbb{R}^p$ such that $\mathbf{X} \mathbf{v} = \mathbf{0}_{p}$. Multiplication of this inequality by $\mathbf{X}^{\top}$ yields $\mathbf{X}^{\top} \mathbf{X} \mathbf{v} = \mathbf{0}_{p}$. As $\mathbf{v} \not= \mathbf{0}_{p}$, this implies that $\mathbf{X}^{\top} \mathbf{X}$ is not invertible.} that the rank of the $(p \times p)$-dimensional matrix $\mathbf{X}^{\top} \mathbf{X}$ is smaller than $p$, and, consequently, it is singular. A square matrix that does not have an inverse is called *singular*. A matrix $\mathbf{A}$ is singular if and only if its determinant is zero: $\mbox{det}(\mathbf{A}) = 0$. - - - -Consider the matrix $\mathbf{A}$ given by: +Super-collinearity of an $(n \times p)$-dimensional design matrix $\mathbf{X}$ implies +that the inverse of the matrix $\hat{X}^T\hat{x}$ (the matrix we needto invert to solve the linear regression equations) is non-invertible. If we have a square matrix that does not have an inverse, we say this matrix singular. The example here demonstrates this !bt \begin{align*} -\mathbf{A} & = \left( +\hat{X} & = \left[ \begin{array}{rr} -1 & 2 +1 & -1 \\ -2 & 4 -\end{array} \right) +1 & -1 +\end{array} \right]. \end{align*} !et -Clearly, $\mbox{det}(\mathbf{A}) = a_{11} a_{22} - a_{12} a_{21} = 1 \times 4 - 2 \times 2 = 0$. Hence, $\mathbf{A}$ is singular and its inverse is undefined. +We see easily that $\mbox{det}(\hat{X}) = x_{11} x_{22} - x_{12} x_{21} = 1 \times (-1) - 1 \times (-1) = 0$. Hence, $\mathbf{X}$ is singular and its inverse is undefined. +This equivalent to saying that the matrix $\hat{X}$ has at least an eigenvalue which is zero. -As $\mbox{det}(\mathbf{A})$ is equal to the product of the eigenvalues -$\nu_j$ of $\mathbf{A}$, the matrix $\mathbf{A}$ is singular if one -(or more) of the eigenvalues of $\mathbf{A}$ is zero. To see this, -consider the spectral decomposition of $\mathbf{A}$: - -!bt -\begin{align*} -\mathbf{A} & = \sum_{j=1}^p \nu_j \, \mathbf{v}_j \, \mathbf{v}_j^{\top}, -\end{align*} -!et -where $\mathbf{v}_j$ is the eigenvector corresponding to $\nu_j$. The inverse of $\mathbf{A}$ is then: -!bt -\begin{align*} -\mathbf{A}^{-1} & = \sum_{j=1}^p \nu_j^{-1} \, \mathbf{v}_j \, \mathbf{v}_j^{\top}. -\end{align*} -!et - -The right-hand side is undefined if $\nu_j =0$ for any $j$. - - -Matrix $\mathbf{A}$ has eigenvalues $\nu_1 =5$ and $\nu_2=0$. According to the spectral decomposition, the inverse of $\mathbf{A}$ is: -!bt -\begin{align*} -\mathbf{A}^{-1} & = \frac{1}{5} \, \mathbf{v}_1 \, \mathbf{v}_1^{\top} + \frac{1}{0} \, \mathbf{v}_2 \, \mathbf{v}_2^{\top}. -\end{align*} -!et - -This expression is undefined as we divide by zero in the second summand on the right-hand side. - - - -In summary, the columns of a high-dimensional design matrix -$\mathbf{X}$ are linearly dependent and this super-collinearity causes -$\mathbf{X}^{\top} \mathbf{X}$ to be singular. Now recall the ML -estimator of the parameter of the linear regression model: +!split +===== Fixing the singularity ===== +If our design matrix $\hat{X}$ which enters the linear regression problem !bt \begin{align} -\hat{\beta} & = (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{Y}. +\hat{\beta} & = (\hat{X}^{T} \hat{X})^{-1} \hat{X}^{T} \hat{y}, \end{align} !et +has linearly dependent column vectors, we will not be able to compute the inverse +of $\hat{X}^T\hat{X}$ and we cannot find the parameters (estimators) $\beta_i$. +The estimators are only well-defined if $(\hat{X}^{T}\hat{X})^{-1}$ exits. +This is more likely to happen when the matrix $\hat{X}$ is high-dimensional. In this case it is likely to encounter a situation where +the regression parameters $\beta_i$ cannot be estimated. - -This estimator is only well-defined if $(\mathbf{X}^{\top} -\mathbf{X})^{-1}$ exits. Hence, when $\mathbf{X}$ is high-dimensional -the regression parameter $\beta$ cannot be estimated. - - -Above only the practical consequence of high-dimensionality is presented: the expression $( \mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{Y}$ cannot be evaluated numerically. But the problem arising from the high-dimensionality of the data is more fundamental. To appreciate this, consider the normal equations: - +The *ad hoc* which was introduced in the 70s was simply to add a diagonal component to the matrix to invert, that is we change !bt -\begin{align*} -\mathbf{X}^{\top} \mathbf{X} \beta & = \mathbf{X}^{\top} \mathbf{Y}. -\end{align*} -!et +\[ +\hat{X}^{T} \hat{X} \rightarrow \hat{X}^{T} \hat{X}+\lambda \hat{I}, +\] +!et +where $\hat{I}$ is the identity matrix. +!split +===== The matrix $\mathbf{X}^{\top} \mathbf{X}$ is of rank $n$, while $\beta$ is a vector of length $p$. Hence, while there are $p$ unknowns, the system of linear equations from which these are to be solved effectively comprises $n$ degrees of freedom. If $p > n$, the vector $\beta$ cannot uniquely be determined from this system of equations. To make this more specific let $U$ be the $n$-dimensional space spanned by the columns of $\mathbf{X}$ and the $p-n$-dimensional space $V$ be orthogonal complement of $U$, i.e. $V = U^{\perp}$. Then, $\mathbf{X} \mathbf{v} = \mathbf{0}_{p}$ for all $\mathbf{v} \in V$. So, $V$ is the non-trivial null space of $\mathbf{X}$. Consequently, as $\mathbf{X}^{\top} \mathbf{X} \mathbf{v} = \mathbf{X}^{\top} \mathbf{0}_{p} = \mathbf{0}_{n}$, the solution of the normal equations is: