Why estimate f?

• prediction: in cases where X are readily available but the output Y cannot be easily obtained.
  • \(\hat{Y}=\hat{f}(X)\)
  • the accuracy of \(\hat{Y}\) as a prediction of Y depends on two quantities: reducible error and irreducible error

In general, \(\hat{f}\) is not a perfect estimate for \(f\), and this inaccuracy will introduce some error. This error is reducible since we can potentially improve the accuracy of \(f\) by using the most appropriate technique. However, since Y is a function of \(\epsilon\), which by definition, cannot be predicted using X, and this is known as the irreducible error.

Why is the irreducible error larger than zero? The quantity \(\epsilon\) may contain unmeasured variables that are useful in predicting Y. The focus of this book is on techniques for estimating \(f\) with the aim of minimizing the reducible error.

• inference: We now are interested in understanding the association between Y and \(X_1, X_2, \dots, X_p\). In this situation, we wish to estimate \(f\), but our goal is not to make predictions for Y. Therefore, \(\hat{f}\) cannot be treated as a black box because we need to know its exact form.

  • which predictors are associated with the response?
  • what is the relationship between the response and each predictor?
  • Can the relationship between Y and each predictor be adequately summarized using a linear equation, or is the relationship more complicated?

How do we estimate f?

• parametric methods: involve a two-step model-based approach

  • first, we make an assumption about the functional form, or shape of \(f\). for example, we assume \(f\) is linear:

  \[ f(X)=\beta_0+\beta_1X_1+\dots+\beta_pX_p \]

instead of estimating the black-box function \(f(X)\), we only need to estimate \(p+1\) coefficients \(\beta_0, \beta_1, \dots, \beta_p\). * second, we need a procedure to train our model using training data, that is we want to find values of estimates of parameters such that

\[ Y \approx \beta_0+\beta_1X_1+\dots+\beta_pX_p \]

the most common approach to fitting/training the linear model is referred to as (ordinary) least squares.

The model-based approach is referred to as parametric: it reduces the problem of estimating f down to one of estimating a set of parameters. Assuming a parametric form for f simplifies the problem of estimating f because it is generally much easier to estimate a set of parameters than it is to fit an entirely arbitrary function f. The downside of parametric approach is that the model we choose will usually not match the true unknown form of f. If the choice of f is too off, then the estimate will be poor. To overcome this, we can choose flexible models that can fit many different possible functional forms for f. But in general, fitting a flexible model requires estimating a greater number of parameters, which can lead to a phenomenon known as overfitting the data, meaning they follow the error o noise too closely.

• non-parametric methods: no assumption about the functional form of f

Instead, they seek an estimate of f that gets as close to the data points as possible. The advantage over parametric approach is by avoiding the assumption of a particular form for f, they have the potential to accurately fit a wider range of possible shapes for f. But non-parametric approaches do suffer from a major disadvantage: since they do not reduce the problem of estimating f to a small number of parameters,a very large number of observations is required in order to obtain an accurate estimate for f, such as splines.

The trade-off between prediction accuracy and model interpretability

Why choosing inflexible method instead of a very flexible approach? * when inference is the goal, linear model may be a good choice since it is quite easy to understand the relationship between Y and X; in contrast, flexible approaches, such as splines and boosting methods can lead to such complicated estimates of f that it is difficult to understand how any individual predictors is associated with the response.

trade-off.

supervised vs. unsupervised

• supervised
  • linear regression
  • logisitic regression
  • GAM
  • boosting
  • support vector machines
• unsupervised
  • clustering
• semi-supervised learning
  • for part of the observations, we don’t have response

regression vs. classification

There is no clear-cut distinction:

• we tend to refer to problems with a quantitative response as regression problems
• those involving qualitative response as classification problems

Least squares linear regression is used with a quantitative response, where as logistic regression is typically used with a qualitative response (binary). Thus, despite its name, logistic regression is a classification method. But since it estimates class probabilities, it can be thought of as a regression method as well.

We tend to select statistical methods based on the fact that if the response is quantitative or qualitative. However, whether the predictors are qualitative or quantitative is generally considered less important, provided that ant qualitative predictors are properly coded before analysis is performed.

measuring the quality of fit

To measure the performance of a statistical learning method on a given data set, we need some ways to measure how well its predictions actually match the observed data. That is, we need to quantify the extent to which the predicted response value for a given observation is close to the true response value for that observation.

in regression setting, the most commonly-used measure is the mean squared error (MSE): \[ MSE=\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{f}(x_1))^2 \] * training MSE * test MSE: we care about this!

When a given method yields a small training MSE but a large test MSE, we are said to be overfitting the data. That is, our statistical learning procedure is working too hard to find patterns in the training data and may be picking up some patterns that are just caused by random chance rather than by true properties of the unknown function of f.

bias-variance trade-off

The expected test MSE, for a given value \(x_0\), can always be decomposed into the sum of three fundamental quantities:
* the variance of \(\hat{f}(x_0)\) * the squared bias of \(\hat{f}(x_0)\) * the variance of the error

A good statistical learning method would simultaneously achieve low variance and low bias. * variance: the amount by which \(\hat{f}\) would change if we estimate it using a different training data set. Since the training data are used to fit the statistical learning method different training data sets will result in a different \(\hat{f}\). If a method has high variance then small changes in the training data can result in large changes in \(\hat{f}\). In general, more flexible statistical methods have higher variance.

• bias: the error introduced by approximating a real-life problem by a much simpler model. Generally, more flexible methods result in less bias.

As a general rule, as we use more flexible methods, the variance will increase and the bias will decrease. The relative rate of change of these two quantities determines whether the test MSE increases or decreases.

The relationship between bias, variance, and test MSE is referred to as the bias-variance trade-off. Good test set performance of a statistical leanring method required low variance as well as low squared biasa. The challenge lies in finding a method for which both the variance and the squared bias are low. in a real-life case in which \(f\) is unobserbed, it is generally not possible to explicitly compute the test MSE, bia, or variance for a statistical learning method. Nevertheless, we should alway keep in mind the bias-variacne trade-off. Cross-validation, which is a way to estimate the test MSE using the training data. Bias-variance trade-off transfer to the classification setting with only some modifications.

classification setting

• training error rate
• test error rate

Suppose that we seek to estimate f on the basis of training observations \({(x_1, y_1), \dots, (x_n,y_n)}\), where \(y_1, \dots, y_n\) are qualitative. The most common approach for quantifying the accuracy of our estimate f is the training error rate, the proportion of mistakes that are made if we apply our estimate \(\hat{f}\) to the training observations:

\[ \frac{1}{n}\sum_{i=1}^{n}I(y_i\ne\hat{y_i}) \]

I is an indicator variable that equals 1 if observed y does not euqal predicted y and 0 if it does. A good classifier is one for which the test error rate is smallest

\[ Ave(I(y_0\ne\hat{y_0})) \]

The Bayes classifier: we should simply assign a test observation with predictor vector \(x_0\) to the class j for which \[ Pr(Y=j|X=x_0) \] is largest. This is a conditional probability and it is the probability that Y=j given the observed predictor vector \(x_0\). This very simple classifier is called Bayes classifier. in a two-class problem where there are only two possible response values, the Bayes classifier corresponds to predicting class one if \(Pr(Y=1|X=x_0)>0.5\), and class two otherwise.

K-nearest Neighbors: * in theory, we would always like to predict qualitative responses using the Bayes classifier. But for real data, we do not know the conditional distribution of Y given X, and so computing the Bayes classifier is impossible. That means, the Bayes classifier serves as an unattainable gold standard against which to compare other methods. * K-nearest neighbors (KNN), like many other approaches, attempts to estiamte the conditional distribution of Y given X. Despite KNN is a very simple approach, it can often produce classifiers that are surprisingly close to the optimal Bayes classifier.

estimating the regression coefficients

We choose \(\beta_0, \beta_1,\dots,\beta_p\) to minimize the sum of squared residuals.

some important questions

when performing multiple linear regression, we usually are interested in answering a few important questions:

• is at least one of the predictors useful in predicting the response?
  • F-test
• do all the predictors help to explain Y, or is only a subset of the predictors useful?
  • forward selection
  • backward selection
  • mixed selection
• how well does the model fit the data?
  • RSE
  • \(R^2\)
• given a set of predictor values, what response value should we predict, and how accurate is our prediction?
  • confidence intervals
  • prediction intervals

qualitative predictors

So far, we have assumed that all variables in our linear regression model are quantitative. In practise, this is not necessarily the case.

• predictors with two levels: indicator or dummy variables using one-hot encoding
• qualitative predictors with more than two levels

extensions of the linear model

• additive: we assume the association between a predictor \(X_j\) and the response Y does not depend on the values of the other predictors
• linear: the change in response associated with a one-unit change in \(X_j\) is constant, regardless of the value of \(X_j\).

interaction term. non-linear relationship: polynomial functional form. * collinearity: correlation matrix for pairs; better way using variance inflation factor(VIF) * outliers: unusual \(y_i\) * high leverage: unusual \(x_i\)

The logistic model

How should we model the relationship between \(p(X)=Pr(Y=1|X)\) and X? we consider using a linear regression model to represent these probabilities: \[ p(X)=\beta_0+\beta_1X \]

The problem with this approach is that with for X close to 0, we predict a negative probability of Y; if we were to predict for very large X, we would get values bigger than 1. These predictions are not sensible, since the true probabilities must fall between 0 and 1.

To avoid this problem, we must model p(X) using a function that gives outputs between 0 and 1 for all values of X. There are many candidate functions that meet this criteria. In logistic regression, we use the logistic function: \[ p(X)=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}} \] To fit the model, we use a method called maximum likelihood. We rewrite the logistic function: \[ \frac{p(X)}{1-p(X)}=e^{\beta_1+\beta_1X} \]

The quantity \(\frac{p(X)}{1-p(X)}\) is called odds, taking the range from \([0, \infty]\). By taking the logarithm of both sides, we arrive at \[ log(\frac{p(X)}{1-p(X)})=\beta_1+\beta_1X \] The left-hand side is called the log odds or logit. We see that the logistic regression model has a logit that is linear in X. What we can observe from the relationship between p(X) and X are:

• the relationship is not linear

• the change in p(X) caused by one unit change in X depends on the actual value of X

• p(X) and X is positively associated if \(\beta_1\) is positive; and vice versa.

estimating the regression coefficients

The coefficients are unknown and must be estimated based on available training data. non-linear least squares can be used to fit the model but maximum likelihood is preferred because of its better statistical properties. We seek estimates for \(\beta_0, \beta_1\) such that the predicted probability \(\hat{p}(x_i)\) of \(x_i\) for each individual corresponds as closely as to the individual’s observed X status. In other words, we try to find such estimates that plugging these estimates into the model for p(X) yields a number close to 1 for all individuals who are classified as group 1, and a number close to 0 for all individuals who are classified as group 0. Formulating it in mathematical equation in a likelihood function: \[ \ell(\beta_0, \beta_1)=\prod_{i:y_i=1} p(x_i)\prod_{i':y_{i'}=1}(1-p(x_{i'})) \] The estimates \(\hat{\beta}_0\) and \(\hat{\beta_1}\) are chosen to maximize this likelihood function. maximum likelihood is a very general approach that is used to fit many of the non-linear models. In the linear regression setting, the least squares approach is a special case of maximum likelihood. One can use qualitative predictors with the logistic regression model using the dummy variable approach.

multiple logistic regression

The results obtained using one predictor may be different from those obtained using multiple predictors when there is correlation among the predictors. In general, the phenomenon is known as confounding.

multinomial logistic regression

when we have a response variable with more than two levels, we extend the two-level case to multinomial logistic regression. To do this, we first select a single level (the Kth level) to serve as the baseline, and then we rewrite the model : \[ Pr(Y=k|X=x)=\frac{e^{\beta_{k0}+\beta_{k1}x_1+\dots+\beta_{kp}x_p}}{1+\sum_{l=1}^{K-1}e^{\beta_{l0}+\beta_{l1}x_1+\dots+\beta_{lp}x_p}} \] for \(k=1,\dots,K-1\), and \[ Pr(Y=K|X=x)=\frac{1}{1+\sum_{l=1}^{K-1}e^{\beta_{l0}+\beta_{l1}x_1+\dots+\beta_{lp}x_p}} \] it can be shown that for \(k=1,\dots,K-1\): \[ log\frac{Pr(Y=k|X=x)}{Pr(Y=K|X=x)}=\beta_{k0}+\beta_{k1}x_1+\dots+\beta_{kp}x_p \] This equation again indicates that the log odds between any pair of classes is linear in the features/predictors.

classifiers with different estimates of the density function

• linear discriminant analysis for p=1
  • assumption: \(f_k(x)\) is normal or Gaussian
  • estimates omitted

LDA classifier plugs in the estimates and assigns an observation \(X=x\) to the class for which \[ \hat{\delta}(x)=x\frac{\hat{\mu_k}}{\hat{\sigma^2}}-\frac{\hat{\mu^2_k}}{\hat{2\sigma^2}}+log(\hat{\pi}_k) \] is largest. The word linear in the classifier’s name stems from the fact that the discriminant function \(\hat{\delta}(x)\) are linear function of x as compared to a more complex form of x.

• linear discriminant analysis for p>1
  • assumption: we assume that \(\vec{X}=(X_1, \dots, X_p)\) is drawn from a multivariate Gaussian distribution, with a class-specific mean vector and a common covariance matrix.

The multivariate Gaussian distribution assumes that each individual predictor follows a one-dimensional normal distribution, with some correlation between each pair of predictors. Put it in mathematical term, we write \(X \sim N(\mu, \mathbb{\Sigma})\). here \(E(X)=\mu\) is the mean of X (a vector of p components), and \(Cov(X)=\mathbb{\Sigma}\) is the \(p \times p\) covariance matrix of X.

in practise, a binary classifier can make two types of errors: * incorrectly assign an individual who belongs to class 1 to class 2 * incorrectly assign an individual of class 2 to class 1

A confusion matrix is a convenient way to display this information.

• sensitivity: the percentage of true class Yes that are predicted
• specificity: the percentage of class No that are correctly identified

The Bayes classifier works by assigning an observation to the class for which the posterior probability \(p_k(X)\) is greatest. In the two-class case this amounts to assigning an observation to class 1 if the posterior probability of class 1 is greater than 0.5. However, if we are concerned about the incorrectly predicting the class 1 status for individuals who belong to class 1, then we can consider lowering the threshold. Choice on the best threshold must be based on domain knowledge, such as detailed information about the costs associated with class 1.

ROC curve is a popular graphic for simultaneously displaying two types of errors for all possible thresholds. ROC is an acronym for receiver operating characteristics.The overall performance of a classifier, summarized over all possible thresholds, is given by the area under the ROC curve (AUC). An ideal ROC curve will hug the top left corner, so the larger the AUC the better the classifier. We expect a classifier that performs no better than chance to have an AUC of 0.5 on test sets. ROC curves are useful for comparing different classifiers since they take into account all possible thresholds.

When we apply a classifier, possible results are:
* true positive (TP) and true negative (TN) * false positive (FP) and false negative (FN)

• false positive rate = FP/N, type I error, 1-specificity
• true positive rate = TP/P, 1-type II error, power, sensitivity
• positive predicted value: TP/P*, precision, 1-false discovery proportion
• negative predicted value: TN/N*

quadratic discriminant analysis

LDA assumes that the observations within each class are drawn from a multivariate Gaussian distribution withi a class specific mean vector and a covariance matrix that is common to all K classes. Quadratic discriminant analysis (QDA) provides an alternative approach.

• assumption: observations are drawn from a Gaussian distribution
• unlike LDA, QDA assumes that each class has its own covariance matrix.

that is, QDA assumes that an observation from the kth class is of the form \(X \sim N(\mu_k, \Sigma_k)\), where \(\Sigma_k\) is a covariance matrix for the kth class. For QDA, the quantity x appears as a quadratic function in the discriminant function. this is where QDA gets its name.

Why assuming different covariance matrices for different classes of Y? or why prefer LDA over QDA? The answer lies in the bias-variance trade-off.

• when there are \(p\) predictors, then estimating a covariance matrix requires estimating \(p(p+1)/2\) parameters.
• QDA estimates a separate covariance matrix for each class, for a total of \(Kp(p+1)/2\) parameters.
• LDA estimates \(Kp\) linear coefficients

Consequently, LDA is a much less flexible classifier than QDA, and has substantially lower variance, which could potentially lead to improved prediction performance. However, if the K classes share a common covariance matrix is badly off, then LDA can suffer from high bias.

In general, LDA tends to be a better bet than QDA if there are relatively few training observations and so reducing variance is crucial. QDA is recommended if the training data set is very large where the variance of the classifier is not a major concern, or the assumption of a common covariance matrix for all K classes is clearly untenable.

naive Bayes

naive Bayes classifier takes a different tract for estimating \(f_1(x), f_2(x), \dots, f_K(x)\). Instead of assuming that these functions belong to a particular family of distributions, e.g., multivariate normal distribution, it makes a single assumption: within the kth class, the p predictors are independent.

Stated mathematically, this assumption means that for \(k=1,\dots, K\), \[ f_k(x)=f_{k1}(x_1) \times f_{k2}(x_2) \times \dots \times f_{kp}(x_p) \] where \(f_{kj}\) is the density function of the jth predictor among observations in the kth class. This is a powerful assumption. Essentially, estimating a p-dimensional density function is challenging because we must consider not only the marginal distribution of each predictor-that is the distribution of each predictor on its own-but also the joint distribution of the predictors-that is, the association among different predictors. in the case of multivariate normal distribution, the association among predictors is summarized by the off-diagonal elements of the covariance matrix. However, in general, this association can be hard to characterize, and exceedingly challenging to estimate. by assuming independence of predictors, we completely eliminate the worry about the association.

In most settings, we do not believe the naive Bayes assumption that all predictors are independent. Essentially, the naive Bayes assumption introduces some bias, but reduces variance, leading to a classifier that works quite well in practise as a result of bias-variance trade-off. When p is large or n is small, naive Bayes relative to LDA or QDA presents a greater pay-off in where reducing variance is very important.

K-nearest neighbors (KNN)

No assumptions are made in KNN so it’s a non-parametric approach: no assumptions are made about the shape of the decision boundary.

• because KNN is completely non-parametric, we can expect this approach to dominate LDA and logistic regression when the decision boundary is highly non-linear, provided that n is very large and p is small
• in order to provide accurate classification, KNN required a lot of observations relative to the number of predictors.
• in settings where the decision boundary is non-linear but n is only modest ,or p is not very small, QDA may be preferred to KNN because QDA can provide non-linear boundary while taking parametric forms, meaning it required a smaller sample size for accurate classification, relative to KNN.
• unlike logistic regression, KNN does not tell us which predictors are important: we don;t get a table of coefficients as from model summary of a logistic regression model.

Poisson regression

Suppose a random variable Y takes on non-negative integer values, and it follows the Poisson distribution, then \[ Pr(Y=k)=\frac{e^{-\lambda}\lambda^k}{k!}, k=0,1,2\dots \] important properties of Poisson distribution are: \(E(Y)=Var(Y)=\lambda\). IN stead of modeling the count directly, we model the mean to vary as a function of the covariates. IN particular, we consider the following model for the mean \(\lambda=E(Y)\), which we write as \(\lambda(X_1,\dots,X_p)\) to emphasize that it is a function of the covariates \(X_1,\dots,X_p\): \[ log\lambda(X_1,\dots,X_p)=\beta_0+\beta_1X_1+\dots+\beta_pX_p \] we take a log of \(\lambda(X_1,\dots,X_p)\) to be linear in predictors in order to ensure that \(\lambda(X_1,\dots,X_p)\) takes on non-negative values for all values of the covariates. To estimate the coefficients \(\beta_0,\beta_1,\dots,\beta_p\), we use the same maximum likelihood approach that we adopt for logistic regression. Specifically, given n observations from the Poisson regression model, the likelihood takes the form \[ \ell(\beta_0, \beta_1, \dots,\beta_p)=\prod_{i=1}^{n}\frac{e^{-\lambda(x_i)}\lambda(X_i)^{y_i}}{y_i!} \] where \(\lambda(x_i)=e^{\beta_0+\beta_1x_{i1}+\dots+\beta_px_{ip}}\). we estimate the coefficients that maximize the likelihood \(\ell(\beta_0, \beta_1, \dots,\beta_p)\).

• interpretation: an increase in \(X_j\) by one unit is associated with a change in \(E(Y)=\lambda\) by a factor of \(exp(\beta_j)\).
  
• mean-variance relationship: by modeling with Poisson regression, we explicitly assume that, in Bike example, mean bike usage in a given hour equals the variance of bike usage during that hour. By contrast, under a linear regression model, the variance of bike usage always takes on constant value. IN general, Poisson regression model is able to handle the mean-variance relationship in a way that the linear regression model is not.
  
• non-negative fitted values: there are no negative predictions using the Poisson regression model. By contrast, when we fit a linear regression model to the Bike data set, almost 10% of the prediction are negative.
  
• over-dispersion: in the Bike example, the variance appears to be much higher than the mean, a situation referred to as over-dispersion. this cause the Z-values to be inflated.

generalized linear models in greater generality

we have covered three types of regression models: linear, logistic, and Poisson. they share some characteristics:

• each approach uses predictors \(X_1, \dots, X_p\) to predict a response Y. we assume that conditional on \(X_1, \dots, X_p\), Y belongs to a certain family of distributions. for linear regression, we typically assume that Y follows a Gaussian or normal distribution; for logistic regression, we assume that Y follows a Bernoulli distribution; for Poisson regression, we assume that Y follows a Poisson distribution.
  
• each approach models the mean of Y as a function of the predictors. in linear regression, the mean of Y takes the form \[ E(Y|X_1, \dots, X_p)=\beta_0+\beta_1X_1+\dots+\beta_pX_p \] i.e., it is a linear function of the predictors. For logistic regression, the mean takes the form \[ E(Y|X_1, \dots, X_p)=Pr(Y=1|X_1, \dots, X_p)=\frac{e^{\beta_{0}+\beta_{1}x_1+\dots+\beta_{p}x_p}}{1+e^{\beta_{0}+\beta_{1}x_1+\dots+\beta_{p}x_p}} \]

while for Poisson regression it take the form \[ E(Y|X_1, \dots, X_p)=\lambda(X_1, \dots, X_p)=e^{\beta_{0}+\beta_{1}x_1+\dots+\beta_{p}x_p} \] Equations above can be expressed using a link function, \(\eta\), which applies a transformation to \(E(Y|X_1, \dots, X_p)\) so that the transformed mean is a linear function of the predictors. that is \[ \eta(E(Y|X_1, \dots, X_p))=\beta_{0}+\beta_{1}x_1+\dots+\beta_{p}x_p \]

The link function for linear, logistic, and Poisson regression are

• \(\eta(\mu)=\mu\)
• \(\eta(\mu)=log(\frac{\mu}{1-\mu})\)
• \(\eta(\mu)=log(\mu)\)

The Gaussian, Bernoulli, and Poisson distributions are all members of a wider class of distributions, known as the exponential family. Other well-known members of this family include exponential distribution, the Gamma distribution, and the negative binomial distribution. In general, we can perform a regression by modeling the response Y as coming from a particular member of the exponential family, and then transforming the mean of the response so that the transformed mean is a linear function of the predictors via the link function. Any regression approach that allows this very general recipe is known as a generalized linear model (GLM).

validation set approach

The validation set approach randomly divides the available set into two parts, a training set and a validation set or hold-out set. The statistical model is fit on the training set, and the fitted model is used to predict the responses for the observations in the validation set. The resulting validation set error rate-mean squared error in the case of quantitative response-provides an estimate of the test error rate.

• the validation estimate of the test error rate can be highly variable, depending on precisely which observations are included in the training set and which observations are included in the validation set.
  
• only observations included in the training set are used to fit the model. That is, only part of the original traning observations are used to fit the model. Statistical learning methods tend to perform worse when trained on fewer observations, which suggests that the validation set error rate may overestimate the test error rate for the model fit on the entire data set.

cross-validation

Cross-validation is a refinement of the validation set approach that addresses these two issues.

best subset selection

to perform best subset selection, we fit a separate least squares regression for each possible combination of the predictors. that is, we fit p models that contain one predictor, \(\binom{n}{k}=p(p-1)/2\) models that contain exactly 2 predictors, and so forth. we then look at all of the resulting models to identify the best predictors.

algorithms is described below:

• let \(M_o\) denote the null model, with no predictors. this model simply predicts the sample mean for each observation.
• for \(k \in 1,2,\dots,p\):
  • fit all \(\binom{p}{k}\) models that contain exactly k predictors
  • pick the best among the \(\binom{p}{k}\) models,a nd call it \(M_k\). best is defined as having the smallest residual sum of squares (RSS), or equivalently largest \(R^2\)-training error.
• select a single best model from among \(M_0, \dots, M_p\) using cross-validated test error, \(C_p(AIC), BIC,\) or adjusted \(R^2\).

in the case of logistic regression, instead of ordering models by RSS in step 2, we use deviance, a measure that plays the role of RSS for a broader class of models. The deviance is negative 2 times the maximized log-likelihood; the smaller the deviance, the better the fit.

disadvantage of best subset selection:
* computational limitation: number of possible models (\(2^p\)) grows rapidly as p increase. if \(p=10\), there are approximately 1000 possible models to be considered; if \(p=20\), then there are over one million possibilities.

solution to computational challenge: step-wise selection. for computational reasons, best subset selection cannot be applied with very large p.

choosing the optimal model

to use above methods to select the best model, we need a way to determine which of these models is best.

note from chapter 2, training error is a poor estimate of test error. therefore training RSS and \(R^2\) are not suitable for selecting the best model among a collection of models with different numbers of predictors since they are related to training error.

to select the best model with respect to test error, we need to estimate test error. there are two common approaches:

• indirectly estimate test error by adjusting to the training error to account for the bias due to over-fitting.
  • \(C_p\)
  • AIC
  • BIC
  • adjusted \(R^2\)
• directly estimate the test error, using either a validation set approach or a cross-validation approach.
  • problem with using validation set approach and cross-validation approach is that by changing the validation set or using another k-fold split, the precise best model will change.
  • solution: select the model using the one-standard-error rule. we first calculate the standard error of the estimated test MSE (RSS/n) for each model size, and then select the smallest model for which the estimated test error is within one standard error of the lowest point on the curve.

ridge regression

recall that the least squares fitting procedure estimates \(\beta_0, \beta_1, \dots, \beta_p\) to minimize \[ \rm RSS= \sum_{i=1}^{n}(y_i-\beta_0-\sum_{j=1}^{p}\beta_jx_{ij})^2 \]

ridge regression is similar to least squares, except that the coefficients are estimated to minimize a slightly different quantity \[ \sum_{i=1}^{n}(y_i-\beta_0-\sum_{j=1}^{p}\beta_jx_{ij})^2+\lambda\sum_{j=1}^{p}\beta^2_j=\rm RSS + \lambda\sum_{j=1}^{p}\beta^2_j \] where \(\lambda \ge 0\) is a tuning parameter, to be determined separately. this equation trades off two different criteria:
- as with least squares, ridge regression seeks to minimize RSS - the second term, called shrinkage penalty, is small when \(\beta_1, \dots, \beta_p\) are close to 0. the tuning parameter \(\lambda\) serves to control the relative impact of these two terms on the regression coefficient estimates.

as \(\lambda \rightarrow \infty\), the impact of the shrinkage penalty grows; when \(\lambda =0\), ridge regression produces least squares estimates. For each \(\lambda\), ridge regression produces a different set of coefficient estimates. so selecting a good value for \(\lambda\) is critical.

when ridge regression has coefficient estimates shrinkun to 0, this correspond to the null model.

term: \(||\hat{\beta}_\lambda^R||_2/||\hat{\beta}||_2 \in [0,1]\). here the notation \(\hat{\beta}\) denotes the vector of least squares coefficient estimates; \(||\hat{\beta}||_2\) denotes the \(\ell_2\) norm of a vector, and is defined as \(||\beta||_2=\sqrt{\sum_{j=1}^{p}\beta_j^2}\), which measures the distance of \(\beta\) from \((0,0,\dots, 0)_{p-1}\).

as \(\lambda\) increases, the \(\ell_2\) norm of \(\hat{\beta}_\lambda^R\) will always decrease, and so will \(||\hat{\beta}_\lambda^R||_2/||\hat{\beta}||_2\). when \(\lambda=0\) when ridge regression penalty term equals RSS, the quantity equals 1; when \(\lambda \rightarrow \infty\) and therefore the ridge regression coefficient estimate is a vector of 0s and \(\ell_2\) norm equal to 0, the quantity equal to 0. A very small value of this quantity indicates that the ridge regression coefficient estimates have been shrunken towards 0.

attention: ridge regression is sensitive to the scaling of predictors. In fact, the value of \(X_j\beta_{j,\lambda}^R\) will depend not only on the value of \(\lambda\), but also on the scaling of the \(j\)th predictor. therefore, it’s best to apply ridge regression after standardizing the predictors so that all standardized predictors all have a standard deviation of 1.

why does ridge regression improve over least squares?

• ridge regression’s advantage over least squares is rooted in the bias-variance trade-off.
  • as \(\lambda\) increases, the flexibility of the ridge regression fit decreases, leading to decreased variance but increased bias.
  • in general, when the relationship between the response and predictors is close to linear, the least squares estimates will have low bias but may have high variance, which means s small change in the training data can cause a large change in the least squares coefficient estimates. Particularly, when \(p \rightarrow n\), the least squares estimates will be extremely variable; if \(p > n\), then the least squares estimates do not even have a unique solution, whereas ridge regression can still perform well by trading off a small increase in bias for large decrease in variance.
  • therefore, ridge regression works best when the least squares estimates have high variance.
  • ridge regression had computational advantages over best subset selection, which requires search for \(2^p\) models. in contrast, for a given \(\lambda\), ridge regression fits only one model and it can be done very quickly.

disadvantages of ridge regression:

• unlike best subset, forward stepwise and backward stepwise selection which generally select models that contain a subset of predictors, ridge regression will include all predictors in the final model
  • the \(\lambda\sum_{j=1}^{p}\beta_j^2\) will shrink all coefficient towards 0, but it will not set any of them exactly to 0 (unless \(\lambda = \infty\)).
  • it may not be a problem for prediction accuracy, but it causes problems to interpretation especially when p is large.

lasso

lasso is a recent alternative to ridge regression that gets rid of noise variables: \[ \sum_{i=1}^{n}(y_i-\beta_0-\sum_{j=1}^{p}\beta_jx_{ij})^2+\lambda\sum_{j=1}^{p}|\beta_j|=\rm RSS + \lambda\sum_{j=1}^{p}|\beta_j| \] comparing lasso to ridge regression, lasso has a similar formulation. the only difference is that we replace \(\beta_j^2\) in the ridge regression penalty with \(|\beta_j|\) in the lasso. in statistical parlance, lasso uses \(\ell_1\) penalty instead of \(\ell_2\) penalty. the \(\ell_1\) norm of a coefficient vector \(\beta\) is defined by \(||\beta||=\sum_{j=1}^{p}|\beta_j|\).

like ridge regression, lasso shrinks the coefficient estimates towards 0. however, the \(\ell_1\) penalty will force some of the coefficient estimates to be exactly 0 when the tuning parameter \(\lambda\) is sufficiently large. Hence, much like best subset selection, the lasso performs variable selection. consequently, models generated from lasso are generally easier to interpret than those produced by ridge regression.

we say lasso yields sparse models-that is, models that contain only a subset of predictors.

• when \(\lambda=0\), lass0 simply gives the least squares estimates
• when \(\lambda \rightarrow \infty\), lasso give the null model where every coefficient estimates are 0
• in between two extremes, lasso is quite different from ridge regression, depending on the value of \(\lambda\), lasso can produce a model involving any number of variables.

re-formularize ridge regression and lasso \[ \rm \]

we can think of above equation as follows: when we perform lasso we are trying to find the set of coefficient estimates that lead to the smallest RSS, subject to the constraint that there is a budget of \(s\) for how large \(\sum_{j=1}^{p}|\beta_j|\) can be.

• when \(s\) is extremely large, then this budget is not very restrictive, and the coefficient estimates can be large. in fact, if s is large enough that the least squares falls within the budget, then it will simply yield the least squares solution
• if s is small, then the penalty must be small to avoid violating the budget.

principal components regression

The first principal component \(Z_1\) direction of the data is that along which the observations vary the most. that is, if we project the observations onto this line, then then resulting projected observations would have the largest possible variance.

the second principal component \(Z_2\) is a linear combination of the variable that is uncorrelated with \(Z_1\), and has largest variance subject to this constraint. it turns out that zero correlation condition of two directions is equivalent to the condition that the direction must be perpendicular, or orthogonal, to the first principal component direction.

advantages of PCR:

• PCR tends to do well in cases when the first few principal components are sufficient to capture most of the variation in the predictor as well as the relationship with the response.
  • when many principal components are required to adequately model the response, PCR performs worse.
  • overall, PCR offers a significant improvement over least squares
  • PCR and ridge regression slightly outperform lasso.
  • PCR is not performing variable selection, since each principal component used in PCR is a linear combination of all the predictors.
  • PCR and ridge regression are closely related.

disadvantage

• PCR suffers from a drawback: there is no guarantee that the directions that best explain the predictors will also be the best directions to use for predicting the response.

partial least squares

PCR involves identifying linear combinations, or directions, that best capture the variance. these directions are identified in an unsupervised way, since the response Y is not used to help determine the principal component directions.

partial least squares(PLS) is an supervised alternative to PCR. likewise, PLS is a dimension reduction method, whcih first identifies a new set of features \(Z_1, \dots, Z_M\) that are linear combinations of the original features, and then fits a linear model via least squares using these M new features. But unlike PCR, PLS identifies these new features in a supervised way with the response involved so tha the M new features are not only approximate the original features well, but also are related to the response. Overall, PLS attempts to find directions that help explain both the response and the predictors.

implementation of PLS:

• after standardizing all predictors, PLS computes the first direction \(Z_1\) by setting each \(\phi_{j1}\) equal to the coefficient fromthe simple linear regression of Y onto \(X_j\). Hence, in computing \(Z_1=\sum_{j=1}^{p}\phi_{j1}X_J\), PLS places the highest weight on the variables that are most strongly related to the response.

• to identify the second PLS direction,

  • we first adjust each of the variable for \(Z_1\), by regressing each variable on \(Z_1\) and taking residuals. these residuals can be interpreted as the remaining information that has not been explained by the first PLS direction.
  • we then compute \(Z_2\) using this orthogonalized data in exactly the same fashion as \(Z_1\) was computed based on the original data.
  • this iterative approach can be repeated M times to identify multiple PLS components
  • finally, we use least squares to fit a linear model to predict Y using all PLS components.

what goes wrong in high dimensions?

for classical statistical methods, least squares for example, that are not intended for high-dimensional setting, when the number of features is as large as, or larger than, the number of observations \(p>n, \rm or, p \approx n\), least squares cannot be performed. the reason is that:

• regardless of whether or not there truly is a relationship between the features and the response, least squares will yield a set of coefficient estimates that result in a perfect fit to the data, such that the residuals are 0.
  • this is problematic because this perfect fit will almost certainly lead to over-fitting.
  • that is although it’s possible to perfectly fit the training data in high-dimensional setting, the resulting linear model will perform extremely poor on an independent test set, and therefore does not constitute a useful model.
  • thus a simple least squares regression line is too flexible and hence overfits the data.
  • \(R^2, C_p\), AIC, BIC are all not applicable in high-dimensional setting

regression in high-dimensions

many of the methods introduced so far, such as forward stepwise selection, ridge regression, lasso, and principal components regression, are particularly useful for performing regression in the high-dimensional setting.

essentially, these methods avoid overfitting bu using a less flexible fitting approach than least squares.

• regularization or shrinkage plays a key role in high-dimensional problems.

• appropriate tuning parameter selectionis crucial for good predictive performance.

• test error tends to increase as the dimensionality of the problem increases, unless the additional features are truly associated with the response.

  curse of dimensionality: one would expect a better performance from a model fitted with more features; however, that is not the case.
  

• adding additional signal features that are truly associated with the response will improve the fitted model, with reduced test set error.

• adding noise that are not truly associated with the response will lead to a deterioration in the fitted model, and consequently an increase test set error. this is because noise features increase the dimensionality, exacerbating the risk of over-fitting.

interpreting results in high dimensions

in high-dimensional settings, the multi-collinearity problem is extreme: any variables in the model can be written as a linear combination of all of the other variables.

essentially, this means that we can never know exactly which variable (if any) truly are predictive of the outcome. at most, we can hope to assign large regression coefficients to variables that are correlated with the variables that truly are predictive of the outcome.

for example, suppose that we are trying to predict blood pressure on the basis of half a million SNPs, and that forward stepwise selection suggests that 17 of those SNPs lead to a good predictive model on the training data. it would be incorrect to conclude that these 17 SNPs predict pressure more effectively than the other SNPs not included in the model. there are likely to be many sets of 17 SNPs that would predict blood pressure just as well as the selected model. if we were to obtain another independent data set and perform the same statistical procedure, we would lilely to obtain a model containing a different, and perhaps even no n-overlapping, set of SNPs. However, this does not detract from the value of the first model obtained. the model may be very effective in predicting blood pressure on an independent set of patients, and might be clinically useful for physicians.

• we must be careful not to overstate the results obtained, and to make it clear that what we have identified is simply one of many possible models for predicting blood pressure, and that it must be further validated on independent data sets.
• it is also important to be particularly careful in reporting errors and measures of model fir in high-dimensional setting. we have seen that when \(p>n\), it is easy to obtain a useless model that has 0 residuals. therefore one should never use sum of squared errors, p-values, \(R^2\) statistics, or other traditional measures of model fit on the training data as evidence of a good model fit in the high-dimensional setting.
  • it is easy to obtain an \(R^2=1\) when \(p>n\).
  • it is important to instead report results on an independent test set, or cross-validation errors. for instance, the MSE or \(R^2\) on an independent test set is a valid measure of model fit, but MSE on the training set certainly is not.

piecewise polynomials

instead of fitting a high-order polynomial over the entire range of X, piecewise polynomial regression fits separate low-degree polynomials over different regions of X. for example, a piecewise cubic polynomial works by fitting a cubic regression model \[ y_i=\beta_0+\beta_1x_1+\beta_2x_i^2+\beta_3x_i^3+\epsilon_i \] where the coefficients differ in different parts of the range of X. the points where the coefficient change are called knots.

for instance, a piecewise cubic with no knots is just a standard cubic polynomial, with \(d=3\); a piecewise cubic polynomial with a single knot at point c takes the form: \[ y_i=\begin{cases} y_i=\beta_{01}+\beta_{11}x_1+\beta_{21}x_i^2+\beta_{31}x_i^3+\epsilon_i, x_i<c \\ y_i=\beta_0+\beta_1x_1+\beta_2x_i^2+\beta_3x_i^3+\epsilon_i, x_i \ge c \end{cases} \] put it in another way, we will fit two different polynomial functions to the data, one on the subset of the observations with \(x_i<c\), and one on the subset of \(x_i \ge c\). note that two functions have different regression coefficients, and can be obtained using least squares. using more knots leads to more flexible piecewise polynomial. in general, for K knots, we will have \(K+1\) different polynomial functions.

problem: when plotting out two polynomial functions, we notice that the function is discontinuous.

solution: add a constraint

constraints and splines

to remedy this problem, we can fit a piecewise polynomial function under the constraint that the fitted curve must be continuous. but still it looks unnatural. we could add another constraint: we not only want the polynomial function to be continuous at the knots, but also very smooth. that is, we want the first and second derivatives of the polynomials are continuous at the knots. this gives us a cubic spline. in general, a cubic spline with K knots uses a total of \(4+K\) degree of freedom.

the general definition of a degree-d spline is that it is a piecewise degree-d polynomial, with continuity in derivatives up to degree \(d-1\) at each know.

spline basis representation

how can we fit a piecewise degree-d polynomial under the constraint that it be continuous?

we can use the basis function model to represent a regression spline. A cubic spline with K knots can be modeled as \[ y_i=\beta_0+\beta_1b_1(x_1)+\beta_2b_2(x_i)+\dots+\beta_{K+3}b_{K+3}(x_i)+\epsilon_i \] for an appropriate choice of basis functions \(b_1, b_2, \dots, b_{K+3}\). the coefficients then can be fit using least squares.

to represent cubic splines, the most direct way is to start off with a basis for a cubic polynomial-namely, \(x, x^2, x^3\)-and then add one truncated power basis function per knot. a truncated power basis function is defined as \[ h(x,\xi)=(x-\xi)^3_+=\begin{cases}(x-\xi)^3, x>\xi\\ 0, \rm otherwise \end{cases} \] where \(\xi\) is the knot. adding a term of the form \(\beta_4h(x,\xi)\) to the cubic regression model will lead to a discontinuity in only the third derivative at \(\xi\); the function will remain continuous, with the first and second derivatives, at each of the knots.

in other words, to fit a cubic spline to a data set with K knots, we perform least squares regression with an intercept and \(3+K\) predictors, of the form \(X, X^2, X^3, h(X, \xi_1), h(X, \xi_2), \dots, h(X, \xi_K))\). this amounts to estimating a total of \(K+4\) regression coefficients; and fitting a cubic spline wit K knots uses \(K+4\) degrees of freedom.

choosing the number and locations of the knots

• when we fit a spline, where should we place the knots?

the splines is most flexible in regions that contain a lot of knots, because in those regions the polynomial coefficients can change rapidly. hence, we could place more knots in places where we feel the function might vary most rapidly, and place fewer knots where it seems more stable. but in practice, it is common to place knots in a uniform fashion. one way to do this is to specify the desired degrees of freedom, and then have the software automatically place the corresponding number of knots at uniform quantile of the data.

• how many knots should we use, or equivalently how many degree of freedom should our spline contain?

one option is to try out different numbers of knots and see which produces the best looking curve. another more objective approach is to use cross-validation. we hold out a proportion of the data, use a spline with a certain number of knots to the remaining data, and then use the spline to make predictions for the held-out data. we repeat this process multiple times until each observation has been left out once., and then compute the overall cross-validated residual sum of squares (RSS). this procedure can be repeated for different numbers of knots K, then the value of K giving the smallest RSS is chosen.

comparing spline to polynomial regression

• polynomial regression produces undesired results at the boundaries, while natural spline still provides a reasonable fit to the data.
• regression splines often give superior results compared with polynomial regression.
  • this is because instead of using high-order, e.g., \(X^{15}\) to produce flexible fits, splines introduces flexibility by increasing the number of knots but keeping the degree fixed.
• generally, spline produces more stable estimates.

in sum, to fit regression splines, we

• specify a set of knots
• produce a sequence of basis functions
• estimate spline coefficients using least squares

choosing the smoothing parameeter

a smoothing spline is simply a natural cubic spline with knots at every unique value of \(x_i\).

in fitting a smoothing spline, we aim at finding the value of \(\lambda\) that makes the cross-validation RSS as small as possible.

GAMs for regression problems

a natural way to extend the multiple linear regression model

\[ y_i=\beta_0+\beta_1x_{i1}+\dots+\beta_px_{ip}+\epsilon_i \]

to allow for non-linear relationshin between predictor and the response is to replace each linear coefficient \(\beta_j x_{ij}\) with a non-linear function \(f_j(x_{ij})\):\[ y_i=\beta_0+f_1(x_{i1})+\dots+f_p(x_{ip})+\epsilon_i \]

it is called an additive model because we calculate a separate \(f_j\) for each \(X_j\), and then add together all of their contributions

pros and cons of GAMs

• GAMs allow us to fit a non-linear \(f_j\) to each predictor, so that we can automatically model non-linear relationships, meaning we don’t have to manually try out transformations on variables.

• non-linear fit can potentially make more accurate predictions for Y.

• able to examine the effect of each predictor while holding all others fixed.

• smoothness of the function \(f_j\) can be summarized via \(df\).

• no interactions are allowed, that is, important interactions can be missed. solution is to manually add interaction terms.

regression trees

Example

• rules: Year < 4.5, Hits < 117.5.
• leaves: sample sub-spaces that are divided by the rules are called terminal nodes or leaves of the tree.

regression trees have an advantage over other types of regression models such as standard least squares regression, ridge regression, lasso: it’s easier to interpret with a graphical representation.

• prediction via stratification of the feature space

  • we divide the feature space into J distinct and non-overlapping regions, \(R_1, R_2, \dots, R_J\).
  • we assign the same prediction value (mean/mode value of the training response in region \(R_j\)) to every observation that falls in the same region \(R_j\).

• objective function: the goal of regression trees is to find a set of regions \(R_1, \dots, R_J\) that minimize \[ \rm RSS = \sum_{j=1}^{J}\sum_{i \in R_j}(y_i-\hat{y}_{R_j})^2 \] where \(\hat{y}_{R_j}\) is the mean response for the training observations in \(j\)th region.

• recursive binary splitting: a top-down, greedy approach

  • top-down: it begins at the top of the tree (with all observations), and then successively splits the predictor space
  • greedy: at each step of the tree-building process, the best split is determined locally at that particular step, rather than globally and looking ahead for future steps

• procedure

  • first select the predictor \(X_j\) and the cutpoint to split the predictor space into regions {\(X|X_j<s\)} and {\(X|X_j \ge s\)} that lead to the greatest reduction in RSS.
  • next repeat the process, and look for the best predictor and cutpoint to split the data further so as to minimize the RSS within each of the resulting regions.
  • after identifying regions, we predict the response for a given test observation by assigning the mean of the training observations within the region where the test observation belongs.

• tree pruning

procedure above may overfit the data and produces poor test set performance due to high complexity of the resulting tree. in contrast, a smaller tree with fewer splits might lead to lower variance and better interpretation at the price of slightly increased bias.

• one alternative is to build the tree only so long as the decrease in RSS due to each split exceeds some threshold.
• a better approach would be building a large tree, and then prune it back to obtain a sub-tree.

how do we decide the best way tp prune the tree? intuitively, we could select a subtree that leads to the lowest test error rate using cross-validation or validation set approach. But estimating cross-validation error for every possible sub-tree is too much since we can have an extremely large number of sub-trees.

• cost complexity pruning, aka weakest link pruning, provides a walk-around.

• algorithm

  • use recursive binary splitting to grow a large tree on training data.
  • apply cost complexity pruning to obtain a sequence of best sub-trees, as a function of \(\alpha\).
  • use k-fold cross-validation to choose \(\alpha\) that minimize the average error.
  • return the sub-tree from step 2 that correspond to the chosen value fo \(\alpha\).

classification trees

similar to regression trees except that it is intended for qualitative responses. prediction for an observation is achieved by assigning the most common class in the region to the observation that falls in the same region.

• objective function
  • classification error rate: not sufficiently sensitive for tree-growing.
  • Gini index: a measure of node purity as it takes on a small value when \(\hat{p}_{mk}\)s are 0 or 1.

\[ \rm G= \sum_{k=1}^{K}\hat{p}_{mk}(1-\hat{p}_{mk}) \] - Entropy: another measure of purity since it approximates 0 when nodes are pure.

\[ \rm D=-\sum_{k=1}^{K}\hat{p}_{mk}log\hat{p}_{mk} \]

Geni index and entropy are preferred when we prune the tree; but if prediction accuracy is the goal, then classification error rate is preferred.

trees vs. linear models

linear regression assumes a model of the form \[ f(X)=\beta_0+\sum_{j=1}^{p}X_j\beta_j \] whilst regression trees assume a model of the form \[ f(X)=\sum_{m=1}^{M}c_m*1_{X \in R_m} \] where \(R_1, \dots,R_M\) represent a partition of feature space. no method dominate the other.

• trees are very easy to explain than linear regression

• trees are believed to more closely mimic human decision-making than classical regression and classification methods

• trees has graphic representations

• no need to create dummy variables for trees

• trees generally have no matching predictive accuracy as non-tree based regression and classification approaches

• trees can be sensitive to small changes, or have high variance, meaning if we split the training set into two subsets at random and fit a decision tree to both subsets, the results could be quite different. A low-variance procedure, such as linear regression if the ratio of n to p os moderately large, will produce similar results when fitting on distinct data repeatedly.

k-means clustering

objective function: a good clustering is one for which the within-cluster variation is as small as possible.

\[ \begin{align*} \text{minimize} \quad & \{\sum_{k=1}^{K}W(C_k)\} \\ \end{align*} \]

to solve this problem, we have to define the within-cluster variation, \(W(C_k)\). the most common choice involves squared Euclidean distance:

\[ W(C_k)=\frac{1}{|C_k|}\sum_{i,i'\in C_k}\sum_{j=1}^{p}(x_{ij-x_{i'j}})^2 \]

where \(|C_k|\) denotes the number of observations in the kth cluster, that is, the within-cluster variation for the kth cluster is the sum of all of the pairwise squared Euclidean distance between the obervatiosn in the kth cluster, divided by the total number of observations in the kth cluster. the optimization problem that defined K-means clustering is given by

\[ \min_{X_1, \dots,X_K} \{ \sum_{k=1}^{K}\frac{1}{|C_k|}\sum_{i,i'\in C_k}\sum_{j=1}^{p}(x_{ij}-x_{i'j})^2 \} \]

hierarchical clustering

one of potential disadvantages of K-means clustering is that it requires pre-specification of the number of clusters K. alternatively, hierarchical clustering does not have such requirement.

• bottom-up/agglomerative clustering: starting from the leaves and combining clusters up to the trunk

• interpretation of hierarchical clustering

  • the height of the fusion on the y-axis indicates how different the two observations are. observations that fuse at the bottom of the tree are quite similar to each other, whereas observations that fuse at the top of the tree will be quite different.
  • misinterpretation: on the basis of the relative location on the y-axis, it is tempting but incorrect to conclude two observations are similar.

• assumption: the term hierarchical refers to the fact that clusters obtained by cutting the dendrogram at a given height are necessarily nested within the clusters obtained by cutting the dendrogram at any greater height.

This is not necessarily true!

• algorithm of hierarchical clustering
  • dissimilarity: measurement of dissimilarity between each pair of observations are defined, such as Euclidean distance.
  • starting at the bottom where each observation is treated as a cluster, and two clusters that are most similar are then fused resulting in n-1 clusters.
  • next two two clusters that are most similar to each other are fused again, resulting in n-1 clusters.
  • proceed until each observation belongs to one single cluster.
• linkage: we need to determine dissimilarity between two clusters/groups
  • complete
  • average
  • single
  • centroid
• choice of dissimilarity measure
  • Euclidean distance
  • correlation-based distance

practical issues in clustering

• small decisions with big consequences
  • should the features be standardized?
  • what distance measure should be used? what type of linkage should be used? where should we cut the dendrogram to obtain clusters?
  • how many clusters should we look for in the data for K-mean clustering?

type I and type II errors

• type I error: the probability of making a type I error given that \(H_0\) holds, i.e., the probability of incorrectly rejecting \(H_0\).
• type II error: we do not reject \(H_0\) when \(H_0\) is in fact falls
• power: 1-type II error, i.e., the probability of not making a type II error, given \(H_a\) holds or the probability of correctly rejecting \(H_0\).

control of family-wise error rate

• Bonferroni method let \(A_j\) denote the event that we make a Type I error for the jth null hypothesis, for \(j=1,\dots, m\), then \[ \begin{aligned} \rm FWER(\alpha)=Pr(\text{falsely reject at least one null hypothesis}) \\ & = \rm Pr(\bigcup_{j=1}^{m}A_j) \le \sum_{j=1}^{m}Pr(A_j) \end{aligned} \]

the inequality is a result from the fact that for any two events A and B, \(\rm Pr(A \bigcup B) \le Pr(A)+Pr(B)\), regardless of whether A and B are independent.

Bonferroni method sets the threshold for rejecting each hypothesis test to \(\alpha/m\), so that for each null hypothesis \(\rm Pr(A_j)\le \alpha/m\). together, \[ \rm FWER(\alpha) \le m \times \alpha/m=\alpha \] so this procedure controls the FWER at level \(\alpha\). For example, to control the family-wise error rate at level 0.1 while testing \(m=100\), null hypotheses, the Bonferroni method requires us to control the Type I error for each null hypothesis at level \(0.1/100=0.001\), and rejecting all null hypotheses for which p-value is below 0.001. - advantages: easy to understand, easy to implement, and successfully controls for Type I error. - cost: being conservative so we reject few null hypotheses, and thus will make quite a few Type II errors.

• Holm’s step-down procedure/Holm-Bonferroni method

Holm’s method is less conservative than Bonferroni’s method while controlling for Type I error.

• special cases: achieving higher power while controlling for Type I error
  • Tukey’s method: when performing \(m=G(G-1)/2\) pairwise comparisons of G means, it allows us to control the FWER at level \(\alpha\) while rejecting all null hypotheses for which the p-value falls below \(\alpha_T\), for some \(\alpha_T>\alpha/m\).
  • Scheffe’s method: \(H_0: \frac{1}{3}(\mu_1+\mu_2+\mu_3)=\frac{1}{2}(\mu_4+\mu_5)\).

trade-off between FWER and Power

power is defined as the number of false null hypotheses that we reject divided by the total number of false null hypotheses \(S/(m-m_0)\). controlling the FWER at level \(\alpha\) guarantees that we are very unlikely to reject any true null hypotheses, or to have any false positives.

in practise, when m is large, we may be willing to tolerate a few false postiveis, in the interest of making more discoveries, i.e., more rejections of the null hypothesis.