TU Wien:Advanced Methods for Regression and Classification VU (Filzmoser)/Possible Exam Questions and Answers

Notes from lecture: Questions

List of questions given in last lecture:

1. Some formulas should be known: (1.1), (2.1), F-test, comparing RSS for model selection / variable selection
2. $R^{2}$ : fit measure, for fit of model to training data; adjusted $R^{2}$ $\rightarrow$ possibility for model selection
3. Information criteria: derivation not important, but result. AIC vs BIC: different penalty, $p$ vs. $log(n)*p$
4. Resampling: CV: 5-, 10-fold. Bootstrapping
5. Least Squares: Linear model, LS estimator: ${\boldsymbol {\hat {\beta }}}=({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}{\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {y}}$. When problematic? Near singularity of $({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})$
6. How to deal with this? $\rightarrow$ Variable selection: AIC, BIC; forward, backward, both; best subset regression
7. PCR principal comp. regression. What and why?
8. Partial Least Squares. Involves also the response variable, unlike PCR!
9. Shrinkage Methods: when is Sum of $\beta _{j}^{2}$ large? In case of near singularity of $({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})$, some $|\beta |$ may be very large, we don’t want that.
10. Formula for ${\boldsymbol {\hat {\beta }}}_{ridge}$ should be known! ${\boldsymbol {\hat {\beta }}}_{ridge}=({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}}+\lambda {\boldsymbol {I}})^{-1}{\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {y}}$, regularization of the solution!
11. Lasso: shrinkage of $|\beta |$. No explicit formula for solution. Why do we get zeros for some of the $\beta _{j}$ ($\rightarrow$ this aims at the figure with the norm balls for $l_{1}$ norm, where LS solution is more likely to hit corners, where some $\beta$ are zero)
12. Classification: Basic formulation using Bayesian theorem, what is posterior, prior prob., likelihood. Assumptions on distributions? $\rightarrow$ mostly gaussians
13. LDA vs QDA: assumptions on the covariance matrices; discriminant functions.
14. RDA: compare to LDA, QDA
15. Logistic Regression: formula: No distribution assumptions, we get a weighted LS problem (Newton-Raphson algorithm)
16. Nonlinear Methods: Splines, what is that? Definition on p.67 (bottom).
17. Problems/parameters: order of splines, number and placement of knots
18. Natural cubic splines, linear extrapolation outside the data
19. Smoothing splines; (7.1) Formula for RSS, penalization of curvature. How to select the degrees of freedom? What is the role of $\lambda$?
20. GAM, GLM: what is that? Formula for PRSS(...) = ....
21. Tree based methods: Regression. p 88: objective function. $Q_{m}(T)=...$, $c_{\alpha }(T)=...$.
22. How to tune between under-, overfitting.
23. Classification: Gini, deviance measure?
24. Random Forests: what are they doing?
25. Bootstrapped data, out of bag data, variable importance?
26. SVM: objective functions, constraints
27. What is the kernel-trick?

Answers

Some Formulae

Classical linear model $y=f(x_{1},x_{2},\ldots ,x_{p})+\varepsilon$ with the linear function: $f(x_{1},x_{2},\ldots ,x_{p})=\beta _{0}+\sum _{j=1}^{p}x_{j}\beta _{j}$ (1.1)

In matrix-vector notation: ${\boldsymbol {y}}={\boldsymbol {X}}{\boldsymbol {\beta }}+{\boldsymbol {\varepsilon }}$

[p. 14f]

The residuals $\varepsilon$ are the difference between the independent (target) variable $y$ and the model prediction ${\hat {y}}$, for every data sample $i=1,\ldots ,n$: $\varepsilon _{i}=y_{i}-{\hat {y_{i}}}=y_{i}-{\hat {\beta _{0}}}+\sum _{j=1}^{p}x_{ij}{\hat {\beta }}_{j}$

This is the cost-function that is to be minimized to get an optimal linear fit (in the least-squares sense) to the data points: $\mathrm {RSS} ({\boldsymbol {\beta }})=\sum _{i=1}^{n}(\underbrace {y_{i}-{\hat {y}}_{i}} _{\varepsilon _{i}})^{2}=(y_{i}-{\hat {\beta _{0}}}+\sum _{j=1}^{p}x_{ij}{\hat {\beta }}_{j})^{2}$

This is one method for model selection.

Testing several of the model coefficients $\beta _{i}$ to be zero, where we have tow different models $M_{0}$ and $M_{1}$, where $p_{0}<p_{1}$. We have the $H_{0}$: “the smaller model $M_{0}$ is true”, i.e. the additional coefficients introduced in the model $M_{1}$ are not significantly different from 0 (CHECK). This yields the following test statistic: $F={\frac {\frac {\mathrm {RSS} _{0}-\mathrm {RSS} _{1}}{p_{1}-p_{0}}}{\frac {\mathrm {RSS} _{1}}{n-p_{1}-1}}}$ (2.1) Under the assumed $H_{0}$, Gaussian residuals ($\varepsilon \sim {\mathcal {N}}(0,\sigma ^{2}$), the $F$ statistic above is distributed according to the $F$-distribution: $F\sim F_{p_{1}-p_{0},n-p_{1}-1},$ with degrees of freedom $p_{1}-p_{0}$ (numerator of ) and $n-p_{1}-1$ (denominator), respectively. If we cannot reject the $H_{0}$, $M_{1}$ is not explaining the data significantly better, if we can reject the $H_{0}$, $M_{1}$ does explain the data significantly better.

$R^{2}$ as fit measure

This is the multiple $R^{2}$, the coefficient of determination, it describes the amount of variance explained by the model. The formula is given as:

${\begin{aligned}R^{2}&=1-{\frac {\mathrm {RSS} }{\mathrm {TSS} }}={\frac {\mathrm {RegSS} }{\mathrm {TSS} }}=\\&=1-{\frac {\sum (y_{i}-{\hat {y}}_{i})^{2}}{\sum (y_{i}-{\bar {y_{i}}})^{2}}}=\\&=\mathrm {Cor} ^{2}(y,{\hat {y}}).\end{aligned}}$

The last term is the square of the correlation between the data $y$ and ${\hat {y}}$, that is the fitted values (the predicted values on the training set variables/observations), note that $\mathrm {Cor} ^{2}(y,{\hat {y}})\in [0,1]$.

The multiple $R^{2}$ grows with the size $p$ of the model. Therefore, we use the adjusted $R^{2}$, that prevents this effect and is thus more suitable for comparing the fit of models of different sizes $p$ (number of variables). This is achieved by including the degrees of freedom in the calculation. ${\tilde {R}}^{2}=1-{\frac {\frac {\mathrm {RSS} }{n-p-1}}{\frac {\mathrm {TSS} }{n-1}}}$ Remember: $n\ldots$ number of observations, $p\ldots$ number of variables (not including the intercept).

$R^{2}$ is a measure mostly for evaluating fit to training data (PF’s comment in lecture recap) .

Information criteria AIC, BIC for model selection

Used for model selection; strictly only valid for nested models (what exactly does that mean?).

Value itself is meaningless, only in comparison between different models it is useful.

Tries to minimize the log-likelihood function of the parameter of a distribution (??), given data. In the lecture notes, there is a derivation for Gaussian likelihood function.

If the residual variance $\sigma ^{2}$ is unknown, it can be estimated using ML method: ${\hat {\sigma }}^{2}=\mathrm {RSS} /n$.

The log-likelihood function can then be expanded to this: ${\begin{aligned}\mathrm {AIC} &=n\log(\mathrm {RSS} /n)+2p\\&={\frac {\mathrm {RSS} }{\sigma ^{2}}}+2p\quad {\text{if }}\sigma ^{2}{\text{ is known}}.\end{aligned}}$

${\begin{aligned}\mathrm {BIC} &=-2\max \log L+\log(n)p\\&=n\log(\mathrm {RSS} /n)+\log(n)p\\&={\frac {\mathrm {RSS} }{\sigma ^{2}}}+\log(n)p\quad {\text{if }}\sigma ^{2}{\text{ is known}}.\end{aligned}}$

AIC and BIC thus look similar, but have a slightly different penalization of the number of variables $p$, and BIC also includes the log of the number of data points, $\log(n)$.

For $n\geq 8$, BIC penalizes more strongly, thus leads to smaller models.

Resampling: k-fold Cross Validation, Bootstrapping

Use loss function $L$, typically MSE.

Used instead of test set, to measure performance of model on unseen data (data it was not trained on).

The given $n$ data samples are split into $q$ parts. One part $k$ is chosen to be the test set, the other $q-1$ parts are used as training data. For the $k$-th part, we calculate the prediction error, e.g. MSE. So, we get basically $q$ different models (estimates ${\hat {f}}$ of the modelled function $f$), and thus $q$ different test set prediction errors.

Note that we get for each $n$ exactly one predicted value! From these we calculate a CV error, as average over the error for each example:

${\widehat {\mathrm {Err} }}_{\mathrm {CV} }={\frac {1}{n}}\sum _{i=1}^{n}L(y_{i},{\hat {y}}_{i})$

Types:

• $q=n$, leave-one-out CV, LOOCV. Has high computational effort, high variance due to the $n$ training sets being very similar
• $k$-fold CV with $k=5$ or $k=10$.

Randomly draw with replacement data sets from the original data. The bootstrap samples have the same size as the original data set. This is done $q$ times, i,e, $q$ bootstrap samples. To each of the bootstrap samples, a model is fit, and evaluated with the loss function. I.e., in contrast to CV, there is no (simulated) test set here? We simply evaluate the fit onto the bootstrap (training) data?

${\widehat {\mathrm {Err} }}_{\mathrm {Boot} }={\frac {1}{q}}{\frac {1}{n}}\sum _{k=1}^{q}\sum _{i=1}^{n}L(y_{i},{\hat {f}}_{k}({\boldsymbol {x}}_{i})).$

Bootstrap error is usually too optimistic.

There is also leave-one-out bootstrap, see lecture notes section 2.4.2, "Bootstrap", page 11 (lecture notes version WS 2018, version from January 2019).

Least Squares estimation, Linear model

In matrix-vector notation: ${\boldsymbol {y}}={\boldsymbol {X}}{\boldsymbol {\beta }}+{\boldsymbol {\varepsilon }}$ LS-estimator for the coefficient vector: ${\hat {\boldsymbol {\beta }}}=({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}{\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {y}}$ Fitted values ${\hat {\boldsymbol {y}}}$: ${\begin{aligned}{\hat {\boldsymbol {y}}}&={\boldsymbol {X}}{\hat {\boldsymbol {\beta }}}\\&={\boldsymbol {X}}\underbrace {({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}{\boldsymbol {X}}{^{\mathsf {T}}}} _{\hat {\boldsymbol {\beta }}}{\boldsymbol {y}}\\&=\underbrace {{\boldsymbol {X}}({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}{\boldsymbol {X}}{^{\mathsf {T}}}} _{{\text{hat matrix}}\,{\boldsymbol {H}}}{\boldsymbol {y}}\end{aligned}}$

How to deal with problems in ordinary LS, near singularity of ${\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}}$?

[p. 20]

For (highly) correlated variables, $({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})$ becomes almost singular, the inverse $({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}$ will become very large (large entries), and the statistical tests for the coefficients to be zero become unreliable. We can thus not make a reliable inference, and the coefficients are not uniquely defined, there is more than one solution.

Note that the inverse of the matrices shows up as the covariance matrix of the coefficient vectors’ distribution; thus, if this gets very large near singularity, the covariance will blow up: ${\hat {\boldsymbol {\beta }}}\sim {\mathcal {N}}_{p+1}({\boldsymbol {\beta }},\sigma ^{2}({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}).$ This is valid for our usual assumptions about the residuals, being zero mean iid uncorrelated Gaussian variables.

Variable selection can get around that problem, by selecting variables that are not correlated. We can use methods that take derived values of the input variables instead of the variables themselves, see PCR, PLS (and continuum regression); furthermore, regularization methods can be used (see Lasso, Ridge).

Principal Component Regression (PCR)

Basic idea: decompose the data matrix ${\boldsymbol {X}}$ into its principal components, and use only the subset explaining a certain amount of variance. The regression is then done on (a subset of) the principal components instead of the original data. This will eliminate from the estimation procedure the contributions of variables that are (highly) correlated. The degree of lost “information”, or rather explained variance of the data, can be tuned by selecting more or less principal components. For predicting on the original data space, the model needs to be transformed back. As PCA is not scale invariant, the data ${\boldsymbol {X}}$ should be centered and scaled beforehand (the centers and standard deviations should be kept for re-scaling and adding the means back after transforming back to the original data space, if desired).

In PCR, we only look at the data ${\boldsymbol {X}}$, but not on whether this makes sense for getting better predictions.

Eigenvalue, eigenvector decomposition of the covariance matrix of the data: ${\boldsymbol {\mathit {\Sigma }}}={\boldsymbol {V}}{\boldsymbol {A}}{\boldsymbol {V}}{^{\mathsf {T}}}$ where ${\boldsymbol {V}}$ is the orthogonal ($p\times p)$ matrix of the eigenvectors and ${\boldsymbol {A}}$ is the diagonal matrix of the eigenvalues of the covariance matrix. The eigenvalues are also the variances of the components $z_{j}$ ($j=1,\ldots ,p$) of ${\boldsymbol {z}}={\boldsymbol {V}}{^{\mathsf {T}}}{\boldsymbol {x}}$.

We can select only $q<p$ principal components for the regression, to use e.g. only a the components adding up to a certain total explained variance.

The model looks like:

${\begin{aligned}{\boldsymbol {y}}&={\boldsymbol {X}}{\boldsymbol {\beta }}+{\boldsymbol {\varepsilon }}\\&={\boldsymbol {X}}{\boldsymbol {V}}{\boldsymbol {V{^{\mathsf {T}}}}}{\boldsymbol {\beta }}+{\boldsymbol {\varepsilon }}\\&=\\&={\boldsymbol {Z}}{\boldsymbol {\theta }}+{\boldsymbol {\varepsilon }}\end{aligned}}$

The model using only $q<p$ variables: ${\boldsymbol {y}}={\boldsymbol {Z}}_{1:q}{\boldsymbol {\theta }}_{1:q}+{\tilde {\boldsymbol {\varepsilon }}}.$

LS estimate:

${\hat {\boldsymbol {\theta }}}={\boldsymbol {Z}}_{1:q}{^{\mathsf {T}}}{\boldsymbol {y}}.$

Transforming back

${\tilde {\hat {\boldsymbol {\beta }}}}={\boldsymbol {V}}_{1:q}{^{\mathsf {T}}}{\hat {\boldsymbol {\theta }}}_{1:q}.$

TO DO: CHECK and reformulate if necessary.

Partial Least Squares (PLS)

Data ${\boldsymbol {X}}$ and targets ${\boldsymbol {y}}$ should be centered. (WHY?)

In contrast to PCA, we try to find projections that maximize (besides the variance of the components?) the covariance between the projections and the predicted variable, but keeps different projections uncorrelated from each other.

For $k=1,2,\ldots ,q$: ${\boldsymbol {a}}_{k}={\underset {\boldsymbol {a}}{\arg \max }}\,{\mathrm {Cov} }({\boldsymbol {y}},{\boldsymbol {X}}{\boldsymbol {a}})$ under constraints $\vert {\boldsymbol {a}}_{k}\vert =1$ and ${\mathrm {Cov} }({\boldsymbol {X}}{\boldsymbol {a}}_{k},{\boldsymbol {X}}{\boldsymbol {a}}_{j})=0$ for $1\leq k<j$ (i.e. $i\neq j$).

Each ${\boldsymbol {a}}_{k}$ makes a weighted combination of the columns of ${\boldsymbol {X}}$ into a single column that is highly correlated with ${\boldsymbol {y}}$. They are then stacked into a matrix ${\boldsymbol {A}}=({\boldsymbol {a}}_{1},\ldots ,{\boldsymbol {a}}_{q})$.

The model looks like:

${\begin{aligned}{\boldsymbol {y}}&={\boldsymbol {T}}{\boldsymbol {\gamma }}+{\boldsymbol {\varepsilon }}\\&=({\boldsymbol {X}}{\boldsymbol {A}}){\boldsymbol {\gamma }}+{\boldsymbol {\varepsilon }}\\&={\boldsymbol {X}}({\boldsymbol {A}}{\boldsymbol {\gamma }})+{\boldsymbol {\varepsilon }}\\&={\boldsymbol {X}}{\boldsymbol {\beta }}+{\boldsymbol {\varepsilon }}\end{aligned}}$

Shrinkage Methods: Lasso, Ridge

Ridge and Lasso are methods for regularization, i.e., stabilizing LS methods for estimating regression coefficients in cases of (highly) correlated variables and thus cases where $({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}$ is close to singularity (see also above).

We penalize the coefficients ${\boldsymbol {\beta }}$. Why do we do that, why are we concerned about large values of / within ${\boldsymbol {\beta }}$? Because we get that if we are near singularity of $({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}})^{-1}$.

Ridge Regression, cost function

[p. 23f]

The cost function looks similar to classical RSS, but has as important difference an additive regularization term: ${\hat {\boldsymbol {\beta }}}_{\mathrm {Ridge} }={\underset {\boldsymbol {\beta }}{\mathrm {argmin} }}\left\{\underbrace {\sum _{i=1}^{n}\left(y_{i}-\beta _{0}-\sum _{j=1}^{p}x_{ij}\beta _{j}\right)^{2}} _{\text{``classical'' RSS}}+\lambda \sum _{j=1}^{p}\beta _{j}^{2}\right\}$ (3.5) $\lambda$ is called complexity parameter and controls the amount of shrinkage. Larger $\lambda$ means more shrinkage. Coefficients are shrunk towards zero and thus towards each other. The regularization term penalizes the sum of the squared coefficients, i.e. the $L_{2}$-norm of the coefficient vector: $\vert {\boldsymbol {\beta }}\vert _{2}^{2}=\sum _{j=1}^{p}\beta _{j}^{2}$ (Note: the intercept $\beta _{0}$ is not included in the penalization term! So this is a slight abuse of notation to write it as the coefficient vector here).

Can be re-written as classical RSS cost function to be solved under the constraint $\sum _{j=1}^{p}\beta _{j}^{2}\leq s,$ with some finite positive $s$.

More important is probably the solution to the cost function written in matrix-vector notation: ${\hat {\boldsymbol {\beta }}}_{\mathrm {Ridge} }=({\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}}+\lambda {\boldsymbol {I}})^{-1}{\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {y}}.$ This makes evident the adding of values $\lambda$ to the diagonal of ${\boldsymbol {X}}{^{\mathsf {T}}}{\boldsymbol {X}}$ before inverting, to ensure non-singularity.

• Ridge solution is not scale invariant, so standardize input.
• Intercept $\beta _{0}$ not included in penalization term.
• Both points above mean that the Ridge regression estimator is not affine equivariant (not sure this was discussed in this class, see e.g. Multivariate Statistik).

There is a relation between Ridge regression and PCR.

TO DO: how to choose $\lambda$, using CV or something? Check notes.

Lasso Regression

The cost-function is again classical RSS with a side constrain: ${\hat {\boldsymbol {\beta }}}_{\mathrm {Lasso} }={\underset {\boldsymbol {\beta }}{\mathrm {argmin} }}\sum _{i=1}^{n}\left(y_{i}-\beta _{0}-\sum _{j=1}^{p}x_{ij}\beta _{j}\right)^{2}$ (x.x) subject to: $\sum _{j=1}^{p}\vert \beta _{j}\vert \leq s.$ Note that again, the intercept is not included in the penalization term. Here, the sum of the absolute values is constrained, i.e. the $L_{1}$-norm of the coefficient vector (without intercept). This is the difference to the Ridge regression case (in terms of cost function).

We get a sparse solution, i.e. we get zeros in the coefficient vector ${\boldsymbol {\beta }}$. Zeros in the coefficient vector basically means variable selection, getting a smaller model.

But, in contrast to Lasso, we get no closed form solution, but have to resort to iterative solvers, e.g. coordinate descent aka. shooting algorithm (?)

TO DO: $s_{0}$, how to select good $s$, why do get coefficients to 0 more likely/faster (norm balls).

Classification, Bayes Theorem

Q: Name the involved parts, i.e. posterior, prior probability, likelihood function. Assumptions on the involved (family of) distributions?

A quick refresher of (one specific incarnation of) the Bayesian theorem. Here, $f(\cdot )$ is a probability density function (pdf), $P(\cdot )$ is a probability, $A$ is some event and ${\boldsymbol {x}}$ is a $p$-dimensional observation: $P(A|{\boldsymbol {x}})={\frac {f({\boldsymbol {x}}|A)P(A)}{\sum \limits _{i}f({\boldsymbol {x}}|A_{i})P(A_{i})}}.$ $P(A)$ is the prior probability of the event $A$, $f({\boldsymbol {x}}|A)$ is a conditional pdf of ${\boldsymbol {x}}$, conditioned on the event $A$; it is termed likelihood function in the context of Bayesian methods. $P(A|{\boldsymbol {x}})$ is the posterior (or a-posteriori) probability of the event $A$ given the observation ${\boldsymbol {x}}$. The denominator follows from the law of total probability and normalizes the numerator to the range of a probability, i.e. $[0,1]$ (the $A_{i}$ have to be mutually disjoint and their union must cover the whole probability space ${\mathit {\Omega }}$. The event $A$ from the numerator is of course included in the $A_{i}$). We will see that the denominator is not important here, don’t worry too much about it.

Let’s specialize this formula to the notation given in the lecture notes [p. 48]: $P(G=k|{\boldsymbol {x}})={\frac {h_{k}({\boldsymbol {x}})\pi _{k}}{\sum \limits _{l=1}^{K}h_{l}({\boldsymbol {x}})\pi _{l}}}.$ Our event $A$ here takes the form of the statement “the class (or group) $G$ is $k$”, where $k=1,2,\ldots ,K$ simply labels one of the $K$ possible groups that the observations ${\boldsymbol {x}}$ belong to. Now, $P(G=k|{\boldsymbol {x}})$ is the (posterior) probability that, given a specific ${\boldsymbol {x}}$, the group label is some number $k$. $h_{k}({\boldsymbol {x}})$ is the likelihood function, and simply uses a different notation, think of this as simplifying the notation like so: $f({\boldsymbol {x}}|G=k)\rightarrow h_{G=k}({\boldsymbol {x}})\rightarrow h_{k}({\boldsymbol {x}})$. This notation no longer makes it clear that this is a conditional pdf, but it is. For the prior probability, the notation is changed to: $P(G=k)\rightarrow \pi _{k}$.

Here, we assume the involved probability distributions to be multivariate Gaussians (note that a likelihood function, being a conditional pdf, has all the properties of a pdf). We denote the multivariate Gaussian pdf, fully parameterized by a mean vector ${\boldsymbol {\mu }}$ and covariance matrix ${\boldsymbol {\Sigma }}$, as $\varphi ({\boldsymbol {x}})$. Note that, as this is a conditional pdf conditioned on a certain group $k$, we index the pdf and the parameters accordingly ($\varphi _{k}$, ${\boldsymbol {\mu }}_{k}$, ${\boldsymbol {\Sigma }}_{k}$): $\varphi _{k}({\boldsymbol {x}})={\frac {1}{\sqrt {(2\pi )^{p}|{\boldsymbol {\Sigma }}_{k}|}}}\exp {\left\{-{\frac {1}{2}}({\boldsymbol {x}}-{\boldsymbol {\mu }}_{k}){^{\mathsf {T}}}{\boldsymbol {\Sigma }}_{k}^{-1}({\boldsymbol {x}}-{\boldsymbol {\mu }}_{k})\right\}},$ where $|{\boldsymbol {\Sigma }}_{k}|$ is the determinant of ${\boldsymbol {\Sigma }}_{k}$.

This is the basis we need for formulating Linear and Quadratic Discriminant Analysis (LDA, QDA). For the derivations, different assumptions on the covariance matrices ${\boldsymbol {\Sigma }}_{k}$ are made to achieve this, see below. In general, the necessary parameters $\pi _{k}$, ${\boldsymbol {\mu }}_{k}$ and ${\boldsymbol {\Sigma }}_{k}$ are estimated from a set of training examples ${\boldsymbol {x_{i}}}$ ($n_{k}$ many per group $k$) that are labelled with their corresponding group $k$ (we need these parameters for each group $k$).

As a teaser, we will see:

• LDA: all covariance matrices are assumed equal, i.e. ${\boldsymbol {\Sigma }}_{k}={\boldsymbol {\Sigma }}$.
• QDA: covariances matrices unequal for the groups: ${\boldsymbol {\Sigma }}_{k}$, for $k=1,2,\ldots K$.

LDA vs. QDA

The main assumption behind LDA (besides Gaussian distributions) is that the covariance matrix ${\boldsymbol {\Sigma }}$ is identical for all groups $k=1,2,\ldots K$. This means that we assume all groups (clusters) in the data to have the same shape (remember how the covariance matrix defines the shape and direction of axes of the Gaussian pdf). This may be a limitation of the capacity to model the data, but we need to estimate far less parameters in comparison to QDA (more on that later). Thus, there is a specific way to estimate the overall pooled covariance matrix ${\boldsymbol {\Sigma }}$ from the training data (it is not group specific and thus has no index $k$): ${\hat {\boldsymbol {\Sigma }}}={\frac {1}{n-K}}\sum _{k=1}^{K}\sum _{g_{i}=k}({\boldsymbol {x}}_{i}-{\hat {\boldsymbol {\mu }}}_{k})({\boldsymbol {x}}_{i}-{\hat {\boldsymbol {\mu }}}_{k}){^{\mathsf {T}}},$ where $n\ldots$ number of training examples, $K\ldots$ number of groups, $g_{i}\ldots$ indices $i$ of the training examples belonging to group $g_{i}=k$. Here, the inner sum estimates a covariance matrix for the respective group $k$ over the training examples belonging to that group (up to some scaling), and the outer sum averages all these per-group-covariance matrices to one pooled matrix.

The estimated covariance matrix can also be computed by first centering the columns of the data matrix ${\boldsymbol {X}}$ (i.e., subtracting the estimated mean of each column), and then arranging the data points belonging to each group $k$ into a distinct, centered data matrix ${\boldsymbol {X}}'_{k}$ (there will be $K$ of those). Then, the common pooled covariance matrix can be computed by averaging these matrices: ${\hat {\boldsymbol {\Sigma }}}={\frac {1}{n-K}}\sum _{k=1}^{K}{\boldsymbol {X}}'_{k}{^{\mathsf {T}}}{\boldsymbol {X}}'_{k}.$ Comment: I took this from the lecture note for Multivariate Statistik, if I remember correctly. It is presented there as a method for estimating the covariance matrix, and was adapted here, using averaging, to calculate the pooled common covariance matrix.

For a derivation of the linear discriminant function (it is linear in the observations ${\boldsymbol {x}}$) see the lecture notes. It has the form: $\delta _{k}({\boldsymbol {x}})={\boldsymbol {x}}{^{\mathsf {T}}}{\boldsymbol {\Sigma }}^{-1}{\boldsymbol {\mu }}_{k}-{\frac {1}{2}}{\boldsymbol {\mu }}_{k}{^{\mathsf {T}}}{\boldsymbol {\Sigma }}^{-1}{\boldsymbol {\mu }}_{k}+\log \pi _{k},$ and involves the group mean, the pooled covariance matrix, the prior probability, and the respective observation ${\boldsymbol {x}}$ that we now want to classify, i.e. assign the most probable class label $k$ to. We do this by plugging our specific ${\boldsymbol {x}}$ into the $\delta _{k}$ for all $K$ different $k$ and assign that $k$ that yields the maximum value of all the $\delta _{k}$. This can also be described formally as: $G({\boldsymbol {x}})={\underset {k}{\mathrm {argmax} }}\,\delta _{k}({\boldsymbol {x}}),$ where $G({\boldsymbol {x}})$ is then the $k$ yielding largest $\delta$. (This is a maximum a-posterior (MAP) decision rule.)

For QDA, we do not assume the covariance matrices of the different groups (clusters) to be identical, thus we have one specific covariance matrix ${\boldsymbol {\Sigma }}_{k}$ per group $k$. Clearly, this means we need to estimate more parameters than in the LDA case, but the model is more flexible as we can have different shapes of the groups.

The discriminant function is more complicated than in the LDA case, and is first and foremost no longer linear in the observations ${\boldsymbol {x}}$ (note the quadratic form involving ${\boldsymbol {x}}$, that is basically the squared Mahalanobis distance of ${\boldsymbol {x}}$ and ${\boldsymbol {\mu }}_{k}$): $\delta _{k}({\boldsymbol {x}})=-{\frac {1}{2}}\log |{\boldsymbol {\Sigma }}_{k}|-{\frac {1}{2}}({\boldsymbol {x}}-{\boldsymbol {\mu }}_{k}){^{\mathsf {T}}}{\boldsymbol {\Sigma }}_{k}^{-1}({\boldsymbol {x}}-{\boldsymbol {\mu }}_{k})+\log \pi _{k}.$

Regularized Discriminant Analysis (RDA) vs. LDA, QDA

Compromise between LDA and QDA. It uses a different version of a pooled common covariance matrix, created as a weighted average from the groups’ distinct covariance matrices ${\boldsymbol {\Sigma }}_{k}$ (for $k=1,2,\ldots ,K$ groups (note that these are the matrices estimated from the data): ${\hat {\boldsymbol {\Sigma }}}={\frac {1}{\sum _{k=1}^{K}n_{k}}}\left(\sum _{k=1}^{K}n_{k}{\hat {\boldsymbol {\Sigma }}}_{k}\right),$ where $n_{k}$ is the number of examples in the respective group.

Then, for each group $k$, a regularized (estimated) covariance matrix can be calculated: ${\hat {\boldsymbol {\Sigma }}}_{k}(\alpha )=\alpha {\hat {\boldsymbol {\Sigma }}}_{k}+(1-\alpha ){\hat {\boldsymbol {\Sigma }}},$ where $\alpha \in [0,1]$ balances between LDA ($\alpha =0$, only pooled common covariance matrix; note that this matrix, calculated as stated above, is slightly different as it uses weighting in the averaging) and QDA ($\alpha =1$, only individual covariance matrices). A value for the hyperparameter $\alpha$ can be estimated using cross-validation.

Logistic Regression

As we formulate this as a likelihood problem, we don’t make distributional assumptions (maybe check this).

Supervised method for classification, i.e. we have a labeled (with groundtruth) dataset available. We want to find a linear function that separates the groups (classes).

The posterior probability for a group (class) label, given a data point ${\boldsymbol {x}}$, is $P(G=k|{\boldsymbol {x}})$, for a group label $k=1,2,\ldots ,K$. Further down, the notation of the posterior will change to $P(G=g_{i}|{\boldsymbol {x}}_{i})=p_{g_{i}}({\boldsymbol {x}}_{i};{\boldsymbol {\beta }})$, more on that later. The subscript $g_{i}\in \{1,2\}$ is the group label for groups 1, 2, for each data point ${\boldsymbol {x}}_{i}$, where $i=1,2,\ldots ,n$.

The log-ratios of posterior probabilities for different groups (for each specific data point) are modeled using affine functions $\beta _{0}+{\boldsymbol {\beta }}{^{\mathsf {T}}}{\boldsymbol {x}}$ (see lecture notes). For the case $K=2$ groups (classes), we get: ${\begin{aligned}P(G=1|{\boldsymbol {x}})&={\frac {\exp {(\beta _{0}+{\boldsymbol {\beta }}_{1}{^{\mathsf {T}}}{\boldsymbol {x}})}}{1+\exp {(\beta _{0}+{\boldsymbol {\beta }}_{1}{^{\mathsf {T}}}{\boldsymbol {x}})}}}\\P(G=2|{\boldsymbol {x}})&={\frac {1}{1+\exp {(\beta _{0}+{\boldsymbol {\beta }}_{1}{^{\mathsf {T}}}{\boldsymbol {x}})}}}.\end{aligned}}$ The posteriors must sum to 1.

For estimating the parameters $\beta _{0}$, ${\boldsymbol {\beta }}$, we can us the maximum likelihood (ML) method. As stated above, the notation in the lecture notes for the posterior probabilities changes at this point like so: $P(G=g_{i}|{\boldsymbol {x}}_{i})=p_{g_{i}}({\boldsymbol {x}}_{i};{\boldsymbol {\beta }})$. Furthermore, the vector ${\boldsymbol {\beta }}$ is now assumed to include the intercept coefficient $\beta _{0}$ and is thus of dimension $(p+1)\times 1$. The data points ${\boldsymbol {x}}_{i}$ are then of course also supposed to contain the intercept (have a 1 inserted in the first entry), and have the same dimension as the coefficient vector.

The log-likelihood function $l({\boldsymbol {}})$ is now written as: $l({\boldsymbol {\beta }})=\sum _{i=1}^{n}\log p_{g_{i}}({\boldsymbol {x}}_{i};{\boldsymbol {\beta }}).$ This can be brought to an easier form, see lecture notes.

We sought to maximize the log-likelihood, i.e., finding the coefficients ${\boldsymbol {\beta }}$ that achieves this. This involves taking the derivative of the log-likelihood and setting it to zero.

An iterative solution using the Newton-Raphson algorithm can be found, this involves also the second derivative, the so called Hessian matrix.

The update equation for the coefficients can also be written using matrices and vectors directly. The algorithm used for solving is also called iteratively reweighted Least Squares (IRLS). Update equation (CHECK): ${\hat {\boldsymbol {\beta }}}_{\mathrm {new} }\leftarrow {\hat {\boldsymbol {\beta }}}_{\mathrm {old} }({\boldsymbol {z}}-{\boldsymbol {X}}{\hat {\boldsymbol {\beta }}}_{\mathrm {old} }){^{\mathsf {T}}}{\boldsymbol {W}}({\boldsymbol {z}}-{\boldsymbol {X}}{\hat {\boldsymbol {\beta }}}_{\mathrm {old} }),$ where ${\boldsymbol {z}}=\exp {(\beta _{0}+{\boldsymbol {\beta }}_{1}{^{\mathsf {T}}}{\boldsymbol {x}})}$. The matrix ${\boldsymbol {W}}$ involved is for $K\leq 2$ a diagonal matrix.

I.e., we get a weighted LS solution to the problem.

Note that the target values ${\boldsymbol {y}}$ $(n\times 1)$ that appear in the LS equation are simply indicators for the class membership, i.e., there is for each row ${\boldsymbol {x_{.i}}}$ in the data matrix ${\boldsymbol {X}}$ a corresponding indicator $y_{i}$. For the 2 class case, the indicators are 1 or 0, accordingly.

The lecture notes say that logistic regression is used for modeling and inference. “The goal is to understand the role of the input variables in explaining the outcome”.

TO DO: Check how the inference works, what kind of information is given by the method about the variables in ${\boldsymbol {X}}$?

It turns out that LDA also has a linear function modeling the log-ratio of the posterior probabilities. LDA achieves this by assumptions of equal covariance matrices for all groups, and the assumption of normally distributed data in the groups. However, the ways for estimating the coefficients are different. Logistic regression makes no assumptions about the prior probabilities of the data points, whereas LDA assumes a (multivariate) Gaussian mixture model.

Nonlinear Methods: Splines

Splines are piecewise polynomial functions. For order $M$, it has continuous derivatives up to order $M-2$. So, for up to 2nd derviative being continous, we need $M=4$, so called cubic splines.

There are different (formulations of?) bases available for splines, the B-splines are numerically more stable.

Parameters for Spline Regression

Parameters of splines:

• the order $M$ of the splines, typically $M\in \{1,2,4\}$
• number of knots $k$
• placement of the $k$ knots. Must (?) be manually chosen, e.g. at appropriate percentiles of $x$

Natural cubic splines

Have the additional constraint that they are linear beyond the boundary knots (outermost knots). Thereby, they avoid the possibly erratic behaviour continuation of the splines beyond the outermost knots (see figures in lecture notes).

Smoothing splines

From all functions $f(x)$ with two continuous derivatives, we seek the one that minimizes the penalized RSS: $\mathrm {RSS} (f,\lambda )=\sum _{i=1}^{n}\left(y_{i}-f(x_{i})\right)^{2}+\lambda \int f''(t)^{2}{\mathrm {d} }t.$ Here, $f''$ is the (continuous) second derivative of the function $f(x)$. $\lambda$ is the smoothing parameter that balances between the influence of the “classical” RSS (that measures the fit of the function to the data) and the penalization term, that penalizes the function for its curvature. Two special cases:

• $\lambda =0$: $f$ is any function that fits the data (so, it could be very curvy).
• $\lambda =\infty$: LS fit, no second derivative of the function $f$ is tolerable (i.e., only an affine function is possible, linear regression?).

The unique minimizing solution of the RSS cost function above is a natural cubic spline with knots at the $x$-values of the data. The data here are collections of pairs $(x_{i},y_{i})$, $i=1,2,\ldots ,n$?

We use the following setup for our function (natural spline): $f(x)=\sum _{j=1}^{n}N_{j}(x)\theta _{j},\quad x\in {\mathbb {R} }(?),$ i.e., we have $n$ many spline functions, each one with its own coefficients.

The cost function can be rewritten in matrix-vector notation: $\mathrm {RSS} (f,\lambda )=({\boldsymbol {y}}-{\boldsymbol {N}}{\boldsymbol {\theta }}){^{\mathsf {T}}}({\boldsymbol {y}}-{\boldsymbol {N}}{\boldsymbol {\theta }})+\lambda {\boldsymbol {\theta }}{^{\mathsf {T}}}{\boldsymbol {\mathit {\Omega }}}_{N}{\boldsymbol {\theta }}.$ Here: $\{{\boldsymbol {N}}\}_{ij}=N_{j}(x_{i})$, i.e., each of the $n$ basis functions is evaluated at one of the $n$ data points (=knots) $x_{i}$. Thus, ${\boldsymbol {N}}$ has dimension $n\times n$. Furthermore, we have: $\{{\boldsymbol {\mathit {\Omega }}}_{N}\}_{jk}=\int N_{j}''(t)N_{k}''(t){\mathrm {d} }t$, note that the integral yields a scalar. Note that ${\boldsymbol {\mathit {\Omega }}}_{n}$ has dimension $n\times n$ as well. Finally, $\theta$ has dimension $n\times 1$.

The solution is then: ${\hat {\boldsymbol {\theta }}}=({\boldsymbol {N}}{^{\mathsf {T}}}{\boldsymbol {N}}+\lambda {\boldsymbol {\mathit {\Omega }}}_{N})^{-1}{\boldsymbol {N}}{^{\mathsf {T}}}{\boldsymbol {y}}.$ this should be compared to the solution of e.g. ridge regression (see lecture notes). Note the matrix ${\boldsymbol {N}}$, the encoding of the data values into the spline basis functions, takes the place of the data matrix ${\boldsymbol {X}}$, and the regularization term involving ${\boldsymbol {\mathit {\Omega }}}_{N}$, weighted by $\lambda$, that is added before inverting.

Choosing the value for $\lambda$ (degrees of freedom) can be done using cross-validation.

CHECK: $\theta$ plays the role of the $\beta$ here? Check how the other equaitons here look like, etc.

Generalized Additive and Linear Models (GAM, GLM)

GLM $\ldots$ Generalized Linear Models
GAM $\ldots$ Generalized Additive Models
For GLMs, the right hand side is a linear combination of the input variables $x_{1},\ldots ,x_{p}$ (of a $p$-dimensional vector ${\boldsymbol {x}}$). For GAMs, the right hand side is an additive combination of (possibly nonlinear) functions of the input variables.

We have the conditional expectation of the target value $y$, given the corresponding data ${\boldsymbol {x}}$: ${\mathrm {E} }(y|x_{1},x_{2},\ldots ,x_{p})={\mathrm {E} }(y|{\boldsymbol {x}})=\mu _{y}({\boldsymbol {x}}).$ The GAM model looks like this, for a so called link function $g(\cdot )$: $g\left(\mu _{y}({\boldsymbol {x}})\right)=\alpha +f_{1}(x_{1})+f_{2}(x_{2})+\ldots +f_{p}(x_{p}).$ There are different link functions available for different purposes, see lecture notes.

For regression (with Gaussian residuals?), the link function is the identity, i.e. $g(\mu _{y}({\boldsymbol {x}}))=\mu _{y}({\boldsymbol {x}})$.

Cost function: penalized residual sum or squares (PRSS), for pairs $({\boldsymbol {x}}_{i},y_{i})$: $\mathrm {PRSS} (\alpha ,f_{1},f_{2},\ldots ,f_{p})=\sum _{i=1}^{n}\left(y_{i}-\alpha -\sum _{j=1}^{p}f_{j}(x_{ij})\right)^{2}+\sum _{j=1}^{p}\lambda _{j}\int f_{j}''(t_{j})^{2}{\mathrm {d} }t_{j}.$ Larger $\lambda _{j}$ lead to smoother functions.

Independent of the choice of the $\lambda _{j}$, an additive model with cubic splines minimizes the PRSS above.

A unique solution needs the following constraints:

• non-singularity of ${\boldsymbol {X}}$
• $\sum _{i=1}^{n}f_{j}(x_{ij})=0,\quad \forall j\quad \Rightarrow \quad {\hat {\alpha }}={\frac {1}{n}}\sum _{i=1}^{n}y_{i}=\mathrm {avg} (y_{i})$

The solution can be found using the backfitting algorithm. (TO DO: add the general idea of the backfitting algorithm).

Tree-based methods

Trees for regression and for classification.

We have data $({\boldsymbol {x}}_{i},y_{i})$, $i=1,2,\ldots ,n$ for $p$-dimensional vectors ${\boldsymbol {x}}$.

Regression Trees

Partitioning the space of the ${\boldsymbol {x}}$-variables into $M$ regions $R_{1},R_{2},\ldots ,R_{M}$. The response variable $y$ as a constant $c_{m}$ in each Region: $f({\boldsymbol {x}})=\sum _{m=1}^{M}c_{m}I({\boldsymbol {x}}\in R_{m}),$ where $I(\cdot )$ is the index function.

A region $R_{m}$ with $n_{m}$ observations inside is also referred to as the node $m$.

A tree can be grown to a certain size $T_{0}$, where the splitting is stopped when a minimum (terminal ?) node size $n_{m}$ is reached. By pruning, i.e. reducing inner (internal) nodes (an internal node represents a split not at a terminal position in the tree?), one can get sub-trees $T$. The number of terminal nodes (regions) in $T$ is denoted as $\vert T\vert$.

Estimating the constant approximating the response in each region $R_{m}$: ${\hat {c}}_{m}={\frac {1}{n_{m}}}\sum _{{\boldsymbol {x}}_{i}\in R_{m}}y_{i}.$ Node impurity at node $m$: $Q_{m}(T)={\frac {1}{n_{m}}}\sum _{{\boldsymbol {x}}_{i}\in R_{m}}(y_{i}-{\hat {c}}_{m})^{2}.$ Cost complexity measure, to be minimized (for all possible $T$ for a given $T_{0}$?): $c_{\alpha }(T)=\sum _{m=1}^{\vert T\vert }n_{m}Q_{m}(T)+\alpha \vert T\vert =\sum _{m=1}^{\vert T\vert }\sum _{{\boldsymbol {x}}_{i}\in R_{m}}(y_{i}-{\hat {c}}_{m})^{2}+\alpha \vert T\vert .$ For every $\alpha$ there exists a unique smallest sub-tree $T_{\alpha }\subseteq T_{0}$ that minimizes $c_{\alpha }$. The optimal $\alpha$ can be found with cross-validation. The actual tree for a chosen $\alpha$ is then found by a search? (CHECK this)

Classification Trees

Partition the (space of the) ${\boldsymbol {x}}$ variables into $1,2,\ldots K$ classes, according to their respective output variable $y$ of $K$ levels.

Proportion of examples belonging to class $k$ in node $m$: ${\hat {p}}_{m,k}={\frac {1}{n_{m}}}\sum _{{\boldsymbol {x}}_{i}\in R_{m}}I(y_{i}=k)$

The class $k$ for the node $m$, denoted as $k(m)$ is then assigned as follows: $k(m)={\underset {k}{\arg \max }}\,{\hat {p}}_{m,k}.$

Measures of node impurity (also denoted $Q_{m}(T)$?):

• Gini index
• Deviance or cross-entropy $\sum _{k=1}^{K}{\hat {p}}_{m,k}\log {\hat {p}}_{m,k}$
• Misclassification error

Cross-entropy and gini index are differentiable

Random Forests

Consist of many trees (hence forest). The randomness is due to using only bootstrap samples instead of the whole data set, and at each knot only use a random selection of the variables, so that every tree is different.

A bootstrap sample has about 2/3 of distinct observations, the rest is the so called Out of Bag (OOB) data.

For the OOB data, compute the impurity of the whole tree, permute for each variable the data values and compute another (different) tree impurity, use difference as measure for variable importance.

For classification, use majority voting for reaching decision. For regression, use average of forest for prediction.

SVM (Support Vector Machines)

The perceptron algorithm gives a separating hyperplane in case the data are linearly separable, but only one of (in this case) infinitely many possible solutions. Support Vector Machines (SVM) can give an optimal solution for such data.

Distance $M$ as minimum distance of any observation to the separating hyperplane: $M={\underset {i=1,2,\ldots ,n}{\min }}g_{i}\,f({\boldsymbol {x}}_{i}),$ with class label $g_{i}$ of observation ${\boldsymbol {x}}_{i}$.

The margin $2M$ is to be maximized, under side-constraints:

${\begin{aligned}{\underset {{\boldsymbol {\beta }},\beta _{0} \atop \vert {\boldsymbol {\beta }}\vert =1}{\max }}&M,\\{\text{s.t.}}\quad g_{i}\,({\boldsymbol {x}}_{i}{^{\mathsf {T}}}{\boldsymbol {\beta }}+\beta _{0})&\geq M,\quad i=1,2,\ldots ,n.\end{aligned}}$

This means, we seek the parameters ${\boldsymbol {\beta }}$, $\beta _{0}$, where $\vert {\boldsymbol {\beta }}\vert =1$, that maximize $M$ such that all (training) examples ${\boldsymbol {x}}_{i}$ are correctly classified (this is what the last line says). SVMs are sometimes called maximum-margin classifiers, as their cost function seeks to maximize the margin.

For data that are not linearly separable, there is the idea of soft-margin classifiers, that introduces slack variables in the optimization scheme.

TO DO: re-formulating the optimization problem.

The Kernel Trick

Transforming the data to a higher dimensional space, where they can be separated linearly.

Basis expansion of the data ${\boldsymbol {x}}$, e.g. $f({\boldsymbol {x}}=h({\boldsymbol {x}}){^{\mathsf {T}}}{\boldsymbol {\beta }}+\beta _{0}$ (this is the function defining the hyperplane).

${\begin{aligned}f({\boldsymbol {x}}&=h({\boldsymbol {x}}){^{\mathsf {T}}}{\boldsymbol {\beta }}+\beta _{0}\\&=\\&=\sum _{i=1}^{n}\alpha _{i}g_{i}\langle h({\boldsymbol {x}}),h({\boldsymbol {x}}_{i})\rangle +\beta _{0}.\end{aligned}}$

Inner product of functions of the data is repalced by a kernel function $K(\cdot ,\cdot )$: $\langle h({\boldsymbol {x}}_{i}),h({\boldsymbol {x}}_{j})\rangle =:K({\boldsymbol {x}}_{i},{\boldsymbol {x}}_{j}).$

Examples: Linear kernel: $K({\boldsymbol {x}}_{i},{\boldsymbol {x}}_{j})={\boldsymbol {x}}_{i}{^{\mathsf {T}}}{\boldsymbol {x}}_{j}.$ Radial Basis Function (RBF): $K({\boldsymbol {x}}_{i},{\boldsymbol {x}}_{j})=\exp {(-\gamma \vert {\boldsymbol {x}}_{i}-{\boldsymbol {x}}_{j}\vert ^{2})}.$