Data 311: Machine Learning

Lecture 6: Logistic Regression

Dr. Irene Vrbik

University of British Columbia Okanagan

Introduction

• Lately we’ve been focused on the regression problem–that is, the supervised machine learning problem when our response variable is a continuous quantity.

• Machine learning algorithms that aim at predicting a categorical (AKA qualitative) response are referred to as classification techniques.

• In this section, rather than trying to predict a continuous responses, we will be trying to classify observations in to categories (AKA classes)

we often first predict the probabilities

These the posts are the only ones that made sense to me:

• https://online.stat.psu.edu/stat504/lesson/6/6.1

• https://statisticaloddsandends.wordpress.com/2019/10/31/understanding-the-components-of-a-generalized-linear-model-glm/

Why not Linear Regression?

Consider coding categories as numeric, say we had recorded eye colour:

\[ Y = \begin{cases} 1 & \text{Green}\\ 2 & \text{Blue} \\ 3 & \text{Brown} \end{cases} \]

Suppose we fed this into a linear model.

Question: What would a prediction of 1.5 signify?

• Mid-way beween green and blue is brown?

• but avg between blue and green would be brown?

• any number smaller than 1 is green

• any number >3 is brown?

• 1.5-2.5 is blue. But the answer would change if we coded them a different number

• can’t treat numeric codings as numeric values

• think carfeully about how we set it up

• diffent colour 50-50 chance between green and brown.

• what would a 50/50 change been green and blue look like?

Class levels

• Numeric codings for a response fed into a continuous model assumes a natural order, or hierarchy.

• For multi-class categorical data, this is simply inappropriate.

• But what about a binary classification problem?

\[ Y = \begin{cases} 1 & \text{Green}\\ 0 & \text{Not green} \\ \end{cases} \]

Can we use this in a linear model such that a predicted value of 0.5 suggests a 50-50 chance that a person’s eye colour is green?

Promising? treat predicted ys as probablities

Example: body

• Consider the body (from the gclus package) data set from previous lectures

• Let’s attempt to model the recorded Gender as a linear function of Height.

# install.packages(gclus) # DONT ever include in Rmd/Rscript
library(gclus)            # cannot data() before library() 
data(body); attach(body)  # cannot attach() before data()
fit1 <- lm(Gender~Height) # cannot call Gender/Height before attach()
plot(Gender~Height)       # note the formula option Y~X

A scatterplot of Gender on the y-axis and Height on the X axis. Points fall along two horizontal lines: where Gender = 1 (indicating males) and Gender =0 (indicating females).

Example: body

Multiple Linear Regression

Recall the linear regression model: \[\begin{equation} Y = \beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p + \epsilon \end{equation}\]

• \(Y\) is the random, quantitative response variable
• \(X_j\) is the \(j^{th}\) random predictor variable
• \(\beta_0\) is the intercept
• \(\beta_j\) from \(j = 1, \dots p\) are the regression coefficients
• \(\epsilon\) is the error term

Two basic assumptions of this model are the additivity in the parameters, and the normally distributed errors https://bookdown.org/pingapang9/linear_models_bookdown/logistic.html

Fitted line

plot(Gender~Height, ylab = "Probability of Male")
abline(fit1, lwd=2, col=2)

Question: What’s wrong with this?

(see next slide for hint)

Same Plot (zoomed out)

plot(Gender~Height, ylim=c(-0.5, 1.5), ylab="Probability of Male")
abline(fit1, lwd=2, col=2)

• get predictions >1 and <0 (negative probabilities)

• cap it at 1?

Generalized Linear Model

Rather than modeling continuous \(Y \in (-\infty, \infty)\) using: \[\begin{equation} Y = \underbrace{\beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p}_{\eta} \end{equation}\]

A GLM consists of three components:

1. Systematic component
2. Random component
3. Link function

Notice that there isn’t an error term \(\epsilon\) (it is absorbed into the random component)

1. The systematic component (aka linear predictor \(eta\)) relates the expected value of the response variable to a linear combination of the predictor variables. It’s represented as:
2. The systematic component specifies the probability distribution of the response variable (dependent variable).
3. Link function denoted by \(g\) relates linear predictor (η) to the expected value of the response variable.

Systematic Component

• The systematic component represents the linear combination of the independent variables.

• It is typically expressed as \[\begin{equation} \eta = \beta_0 + \beta_1 x_1 + \dots \beta_p x_p \end{equation}\] where

• \(\eta\) is the linear predictor,

• \(\beta_0, \beta_1, \dots \beta_p\) are the coefficients, and

• \(x_1, x_2, \dots, x_p\) are the predictors.

• A GLM is a flexible and powerful framework for modeling the relationship between a dependent variable and one or more independent variables.
• It generalizes the classical linear regression model by allowing for different types of response variables, including non-continuous data, and by using a link function to connect the linear predictor to the response variable. G so we need a function to map \(- \infty < y < \infty\) to value from 0 to 1

Random Component

• We assume that \(y_1, \dots, y_n\) are samples of independent random variables \(Y_1, \dots, Y_n\) respectively.

• The random component specifies the probability distribution \(f(y_i; \theta_i)\) of the response variable.

• For GLM’s the probability distributions are assumed to arise from the Exponential family

exponential family is a parametric set of probability distributions of a certain form, specified below.

Exponential family

• The exponential family includes a large class of probability distributions e.g. normal, binomial, Poisson and gamma distributions, among others.

• The mean, of the distribution depends on the independent variables, \(X\), through:

  \[\begin{align*} \operatorname {E}[Y_i] &= g^{-1}(\eta)\\ \operatorname \eta &= g({E}[Y_i]) \end{align*}\]

  where \(g(\cdot)\) is the link function

Special cases of GLMs

Classical regression assumes:

• \(Y_i \sim \text{Normal}(\mu_i, \sigma^2)\)
• \(E[Y_i] = \mu_i = g^{-1}(\eta)\)

Logistic regression assumes

• \(Y_i \sim \text{Bernoulli}(\pi_i)\)
• \(E[Y_i] = \pi_i = g^{-1}(\eta)\)

We will studying binary (where \(Y\) is a binary random variable with two possible classes) but this can be extended to \(>2\) classes.

• as its name suggest, logistic regression this is closely related to classical linear regression, only now…

• y-axis Rather than modeling this response \(Y\) directly, logistic regression models the probability that \(Y\) belongs to a particular category.

• DATA - 410: Regression and Generalized Linear Models

• Let \(\mu_i = E[Y_i]\). Often we will not estimate the \(\theta\)s directly, but rather some function of \(\mu_i\) (as we will see soon).

This is the only random component in the model; there is not a separate error term.

where \(\pi\) is the probability that \(Y_i = 1\). Recall that Bernoulli(\(\pi\)) is a special case to the Binomial(\(n\), \(\pi\)) distribution; i.e. when \(n = 1\)

Link function

• The link function is a transformation that relates the mean of the response variable to the linear predictor.

\[\begin{equation} \eta = g({E}[Y_i]) \end{equation}\]

• If \(\eta = g(\mu_i) = \beta_0 + \beta_1 x_1 + \dots \beta_p x_p\) the link function is called the identity link and we get classical regression.

• If \(\eta = g(\pi_i) = \ln \left( \frac{\pi_i}{1-\pi_i} \right)\) we get the logit link function and get logistic regression.

• Odds are traditionally used instead of probabilities in horse-racing, since they relate more naturally to the correct betting strategy

• odds take on values from 0 (very small probs) to positive infinity (when prob = 1)

• log-odds will go from -infty to infty

• log-odds may feel weird but they have a nice interpretation

Logit Function

More generally if \(p\) is a probability, then

• \(\dfrac{p}{1 − p}\) is the corresponding odds and

• logit \(p\) = \(\ln \left( \dfrac{p}{1-p}\right)\) is the logit of the probability \(p\)

• this value is sometimes referred to as the log-odds or “logits”

Since the logit is the natural logarithm of odds it is sometimes referred to as the log-odds.

• odds take on values from 0 (very small probs) to positive infinity (when prob = 1)

• log-odds will go from -infty to infty

• log-odds may feel weird but they have a nice interpretation

Plotted log-odds

Probabilities converted to the logit (log-odds) scale. Notice how the logit functions allows us to map values from 0 to 1 to values from negative infinity to infinty.

If \(p = 0.5 \implies\) the odds \(\left(\frac{p}{1-p}\right)\) are even and the log-odds is zero.

Negative (resp. positive) logits represent probabilities below (resp. above) one half.

• “A curve representing log-odds on the y-axis and probability on the x-axis. The curve approaches negative infinity for very small probabilities (close to 0) and infinity for very large probabilities (close to 1). The curve is a smooth s-shaped rotated by 90 degrees and mirrored.

Activation Function

• The link function \(g(\cdot)\) is an invertible function that transforms \(\operatorname {E}[Y_i]\) to \(\eta\).

• However, in the machine learning community, we are often first introduced to the inverse of the link function \(g^{-1}()\).

• This is sometimes called the activation function (or inverse link function.

• For logistic regression, the activation function \(g^{-1}\) is the standard¹ logisitic function (aka sigmoid function).

1. with the supremum of the values of the function = 1, the logistic growth rate = 1, and function’s midpoint = 0.

Relationship between Link and Activation Function

Link function \[\begin{align} g(\pi_i) &= \eta\\ \text{logit } \pi_i &= \eta\\ \ln \left( \dfrac{\pi_i}{1 - \pi_i} \right) &= \eta\\ \end{align}\]

Activation function \[\begin{align} \pi_i &= g^{-1}(\eta)\\ &= \dfrac{1}{1 + e^{-\eta}}\\ &= \dfrac{1}{1+ e^{-(\beta_0 + \beta_1 X_1 + \cdots \beta_p X_p)}} \end{align}\]

\(g\) is the logit function and \(g^{-1}\) is the logistic function

Both do the same thing, linking the linear predictor \(\eta\) with the \(E[Y_i]\) (in this case the probability of \(Y = 1\))

Logistic Function

\[\begin{align} g(\pi_i) &= \eta \\ \ln \left( \frac{\pi_i}{1-\pi_i} \right) &= \beta_0 + \beta_1 X_1 + \cdots \beta_p X_p \\ \frac{\pi_i}{1-\pi_i} &= \exp\{\beta_0 + \beta_1 X_1 + \cdots \beta_p X_p\} \\ \implies \pi_i &= \dfrac{\exp\{\beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p\}}{1+ \exp\{\beta_0 + \beta_1 X_1 + \cdots \beta_p X_p\}}\\ &= \dfrac{1}{1+ e^{-(\beta_0 + \beta_1 X_1 + \cdots \beta_p X_p)}} = g^{-1}(\eta) = \dfrac{1}{1 + e^{-\eta}}\\ \end{align}\]

This allows us to map the linear predictor \(\eta\) to probabilities (i.e. \(\pi_i\) values from 0 to 1)

Rather than mapping probabilities to the real line, let’s map the linear combination of predictors \(\eta\) to the values ranging from 0 to 1

Notice now that this is not a linear function in terms of probability but in terms of “log-odds”

1. Logistic (or sigmoid) function that maps the linear combination of predictors to values ranging from 0 to 1
2. write this in terms of odds
3. take the log of both sides to represent as log-odds

Plotted Logistic Function

It is a S-shaped curve (sigmoid curve) that allows us to go back from logits to probabilities.

curve representing probabilities on the y-axis and log-odds on the x-axis. The curve approaches 0 large as lod-odds becomes more negative and approaches 1 as the log-odds becomes more positive (when log-odds = 0, probablity = 0.5). The curve is a smooth s-shape.

For any value of log-odds i can convert to a number between 0 and 1. The smaller the logodds, the smaller the probabilitily

From Line to S-curve

Instead of trying to fit a line to this data, let’s try to fit something more like this S-curve so that \(0 \leq \pi_i \leq 1\)

• Fitting a straight line doesn’t make sense so instead we fit an S-curve goes from 0 to 1

• so our predictions will fall on this blue sigmoid curve and correspond to the \(P(Y=1 \mid X)\)

  • so gives the probabilty of male given height

• We need a function that will transform any point on the red line to a point on this blue S-curve

Logistic Regression

Logistic Regression:

\[ \text{log odds} = \beta_0 + \beta_1 X_1 + \dots + \beta_p X_p \]

Assumptions

\[ \begin{align} Y_i &\sim \text{Bernoulli}(\pi_i)\\ g({E} [Y_i]) &= \beta_0 + \sum_{j=1}^p\beta_j X_j = \eta \\ g^{-1}( \eta) &= \dfrac{e^{\sum_{j=1}^p\beta_j X_j}}{1+ e^{\sum_{i=1}^p\beta_i X_i}} = \dfrac{1}{1+e^{-\sum_{j=1}^p\beta_j X_j}}\\ \end{align} \]

where \(g(\cdot)\) is the logit function.

Interpretation of Coefficients

Since probability \(\pi_i\) is non-linear as \(X\) varies¹, we cannot interpret the coefficients as we did in regular regression.

The best we can do is talk about the direction:

• Holding all other variables constant, if \(\beta_j\) is positive, then an increase in \(X_j\) is associated with an increase in \(\pi_i\)

• Holding all other variables constant, if \(\beta_j\) is negative, then an increase in \(X_j\) is associated with a decrease in \(\pi_i\).

• IOW a logistic regression model, increasing X by one unit changes the log odds by β1 (4.4).

• Equivalently, it multiplies the odds by \(e^{\beta_1}\)

• it’ll matter where we are in terms of X (see example)

1. in logistic regression, increasing \(X\) by one unit changes the log odds by \(\beta_1\)

Estimation

Let \(G_1\) be the set of observations where \(Y=1\), let \(G_0\) be the set of observations where \(Y=0\), and let \(\boldsymbol{\beta}\) be the set of coefficients. The likelihood is given by:

\[ l(\boldsymbol{\beta}) = \prod_{i \in G_1 }P(Y=1\mid x_i)\prod_{h \in G_0 }(1-P(Y=1\mid x_h)) \]

To fit the parameters of these models, we use maximum likelihood ¹ estimation.

Not through least squares anymore (we don’t have a concept of a residual anymore)

• beta = vector of beta

• likelihood of betas = product of probabilities given the predictors

• one of the most common ways we fit models

• the MLE estimates are the most likely values of the betas given the data we saw

• try to find the best fitting squiggle s

1. Maximum likelihood is a core concept of statistical inference, but is outside the scope of this course (see STAT 401, and some others, instead)

From lm() to glm()

You can fit a logistic regression using the glm() function.

glm(formula, family = gaussian, data, ...)

Arguments	Description
formula	a symbolic description of the model to be fitted, eg Y~ X1 + X2
family	a description of the error distribution and link function to be used in the model ?family
data	data frame (usually your training set)

Click things to go to the help file!

glm families

• To perform logistic regression with a binary outcome we need to set family = "binomial"¹.

• Other options (not be covered in this course) include:

  • family = "gaussian" same as lm()

  • family = "poisson" for predicting counts

  • family = multinomial for \(>2\) classes

  • family = binomial('probit') probit regression (a common alternative to logistic regression)

• poisson for modelling count data
• multinomial for modeling >2 levels in a categorical variable
• logit is not the only function we can use to map the linear predictor to a probability (common alternative is probit - inverse CDF of a normal distribution)
• In practice, both probit and logistic regression can yield similar results for many datasets, and the choice between them often depends on the specific problem, theoretical considerations, and the preference of the researcher.

Descriptions of families: https://www.youtube.com/watch?v=DDP62EUMRFs&ab_channel=KasperWelbers

glm are very flexible and powerful

1. Think binomial distribution (can be either “success” or “failure”)

Fitted Logistics Regression Plot

We are using the Logistic Function with the fitted \(\hat \beta\) values to plot the fitted S-curve

# Fit the logistic regression model
simlog <- glm(factor(Gender) ~ Height , family="binomial")

# store the beta hats for easy referencing
betas <- coef(simlog)  

# plot the data (renaming the y-axis)
plot(Gender~Height, ylab="Probility of being Male")

# plot the fitted curve p(x) = e^(b0+b2x)/(1 + e^(b0+b2x))
curve(
  (exp(betas[1] + betas[2]*x))/(1+exp(betas[1] + betas[2]*x)),
  add=TRUE, # superimposes onto scatterplot (otherwise new plot)
  lwd=2, # line width (make the line a bit thicker)
  col=2) # 2 = red  

This will be the “best” S curve in terms of the maximum likelihood function.

All probabilities lie between 0 and 1 when log-reg is used

to convert the linear combo of betas to probabilities we have to use the inverse logit, aka the logistic function

• i’m using curve() but you can plot the predict()ed values for a lengthy sequence of plausible x values

Fitted Logistics Regression Plot

Summary Output

summary(simlog)

...
Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -46.76328    4.00571  -11.67   <2e-16 ***
Height        0.27292    0.02339   11.67   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 702.52  on 506  degrees of freedom
Residual deviance: 389.59  on 505  degrees of freedom
AIC: 393.59

Number of Fisher Scoring iterations: 5
...

We can use this output to assess what predictors are useful in classifying samples - the p value corresponds to the Wald test (but similar interpretation to the t-test from regression) to test if it’s statistically significant

interpretation of coefs:

• \(\beta_0\) when height = 0 the log-odds of being a male is -46.76 (nonsensical)

• \(\beta_1\) For every one cm increase in height, you would expect the log-odd of male to increase by 0.2792

• notice how we don’t have an R-2 cuz we don’t have residuals We’ll to compare log reg models in a different way than we did for normal regression

Multiple Logistic Regression

As with linear regression we could just as easily add predictors.

# Fit the logistic regression model
mlr <- glm(factor(Gender) ~ Height + WaistG + BicepG, family="binomial")
summary(mlr)

...
Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -55.32670    5.51960 -10.024  < 2e-16 ***
Height        0.21534    0.02864   7.520 5.50e-14 ***
WaistG        0.05167    0.02891   1.787   0.0739 .  
BicepG        0.46541    0.08265   5.631 1.79e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 702.52  on 506  degrees of freedom
Residual deviance: 232.90  on 503  degrees of freedom
AIC: 240.9

Number of Fisher Scoring iterations: 6
...

we could also include interaction terms

Residual Deviance

• You will notice that we no longer have a \(R^2\) or adjusted R-square value as we did in the linear regression case.

• the concept of “variance explained” doesn’t not apply.
• furthermore, since the observed y values are 0 and 1, transforming those values to the linear form using log(p/1-p) will produce -\(\infty\) and \(\infy\) so there are no ways to calculate residuals

• Deviance is a lack of fit measure (the smaller the better) that plays the role of RSS for a broader class of models.

• Instead of sum of squares, logistic regression uses deviance
• other metrics exist (example psuedo R^2 value)

• The residual deviance is a measure of the deviance that remains unexplained after fitting the logistic regression model.

• The residual deviance tells us how well the response variable can be predicted by a model with p predictor variables.
• The lower the value, the better the model is able to predict the value of the response variable.
• model comparison When comparing two or more logistic regression models, you can use the likelihood ratio test (e.g. anova(mod1, mod2, test = Chisq)), which involves comparing the residual deviance of the models.
• we will see deviance come up again later in the course

The deviance is a function of the log-likelihood

Evaluation

• Evaluating a logistic regression model involves assessing its performance and determining how well it fits the data and makes predictions.

• In Lecture 3 we saw some metrics for assessing regression models

• Logistic regression, however, deals with the classification problem

• In the realm of classification there are different set of evaluation techniques and metrics are used.

Metrics for Classification

In addition to deviance there are other popular metrics that we may consider for the classification problem:

• error/misclassification rate
• accuracy
• precision
• recall
• specificity
• F1-score

These can all be computed from a so-called classification table (aka confusion matrix)

• Also the ROC curve Receiver Operating Characteristic (ROC) curves graphically represent the trade-off between true positive rate (recall) and false positive rate (1-specificity) for different classification thresholds.

• The Area Under the ROC Curve (AUC-ROC) measures the model’s ability to discriminate between classes. A higher AUC-ROC indicates better discrimination.

we will look at these metrics in more detail next lecture

Validation of Predicted Values

With a fitted logistic regression model the predict() function (see ?predict.glm) we can either predict the categorical response or output a probability.

predict(mod_fit, newdata=testing)
predict(mod_fit, newdata=testing, type="prob")

or more generally for glm models

From Probability to Classification

For predicting classification we consider

\[ \max \{ P(Y=0 \mid x), P(Y=1 \mid x) \} \]

We can obtain these probabilities in R using the predict() function on the fitted model (see lab)

Since we are considering only binary, \(P(Y=y \mid x)=0.5\) defines a the boundary.

Question: Where is this boundary for the body example?

• Although, logistic regression, tells the probability that a mouse is obese or not it’s usually used for classification

• For example if the probability of male is greater than 50% then we’ll classify it as obese

Like KNN classification we say, whoever has the highest prob gets assigned that obs (Bayes classifier)

Notation alert: we refer to ML algorithms that automatically order or categorize data into one or more “classes” as a classifier;

Plot of Probability Boundary

Highlight the observations that would be misclassified by this classifier

Height boundary

\[ \begin{align} 0.5 &=\dfrac{e^{\beta_0 + \beta_1 X}}{1+ e^{\beta_0 + \beta_1 X}}\\ \implies \log \left(\frac{0.5}{0.5-0}\right) = \log \left(1\right) = 0 &= {\beta_0 + \beta_1 X} \\ \implies X &= -\frac{\beta_0}{\beta_1} \end{align} \] So for Height \(> -\dfrac{ -46.7632791}{0.2729192} = 171.34\) we would predict Male.

Classification Table

For classification, we often summarize performance in a classification table (aka confusion matrix).

	predicted male	predicted female
1 - male	216	44
0 - female	45	202

In this case, the off-diagonals are errors (and diagonal are correctly classified obs.).

Error Rate

In words, the error rate (aka misclassification rate) calculates the proportion of misclassifications that are made by \(\hat f\):

\[\begin{align} \frac{1}{n} \sum_{i=1}^n I(y_i \neq \hat y_i) = \frac{\text{total misclassified}}{\text{total no. observations}} \end{align}\]

and \(\hat y_i\) is the predicted class label for the \(i\)th observation \(x_i\) using \(\hat f\) and \[\begin{equation} I(y_i \neq \hat y_i)= \begin{cases} 1, y_i \neq \hat y_i\\ 0, \text{otherwise} \end{cases} \end{equation}\]

The error rate for this classification table is \[\begin{align} \dfrac{44+45}{216+44+45+202} = \dfrac{89}{507} = 0.1755424 \approx 17.6% \end{align}\]

no r squared values! (theres about 10 ways people have thought about doing this)

There are other way to access predicted classifications later

Problems

• the regular diagnostics plots don’t really apply for this model

• multiple categories? (back to the reference label)

• we can extend this model (but there are better options )

• next lecture will be more mathy

Testing vs Training error rate

• Akin to our earlier discussions with MSE for regression, we are generally interested with the error rate from the testing set rather than the training set.

• That is, for a set of new observations \((x_0, y_0)\), a good classifier achieves the smallest test error on average: \[\begin{equation} E[I(y_0 \neq \hat y_0))], \end{equation}\] where \(\hat y_0\) is the predicted class label for \(x_0\) that uses \(\hat f\).

Comments

• Is it reasonable to assume a linear relationship between the log-odds and the predictors?

Your guess is as good as mine…(no more residuals to check)

• assuming a linear relationship between log odds and predictors (log odds actually get mapped to -infty and infty so we don’t have the concept of a residual so we can’t do the diagnostics.

• predict Y = 1 Probabilty 1 goes to infinicty, Y = 0 goes to negative infinity

• What if we have multiple categories for our response?

No problem, bit of a change in notation and mathematics, but R can handle it just fine. family = "multinomial"

✏️ Next class we’ll move on to another (more natural) model for classification (Bring some paper and a pen!)

• unstable parameters when the separation between two classes is large (next class of models we consider do not suffer from this problem)

• DA is a more naturally extensions into the multivariate case