Update on regression slides

This commit is contained in:
mhjensen
2018-09-05 22:28:34 +02:00
parent 1b59ce5a0c
commit eb9bc8a6ea
+74 -123
View File
@@ -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: