• Understand the difference between a linear and a generalised linear model
• Be able to run and interpret a binary logistic regression model in \(R\)
• Transform values between different scales (probabilities, odds, log odds)
• Make predictions using model coefficients and \(R\) functions such as predict()

You will need to watch the lecture and read Chapter 10 in Andrews (2021) (pages 346 – 363) to fully grasp the concepts.

• regression models are models that model the probability distribution of one variable, known as the outcome variable, and varies as a function of other variables, known as predictor variables.
• the most basic type of regression model is the normal linear model.
• in normal linear models, we assume that the outcome variable is normally distributed around a mean \(\mu\); the mean varies linearly with changes in a set of predictor variables.

Where \(\mu\) is a linear function of the intercept \(\beta_0\) and its slope \(\beta_1\), relative to the value that the predictor \(x\) takes on.

A simple linear regression – with a single predictor \(x\) – has the form

\[ \mu_i = \beta_0 + \beta_1 \times x_i \]

For multiple linear regression predictors, with \(K\) predictors, the form is

\[ \mu_i = \beta_0 + \beta_1 \times x_{1i} + \beta_2 \times x_{2i} + \cdots + \beta_K \times x_{Ki} \]

or for short

\[ \mu_i = \beta_0 + \sum_{k=1}^K \beta_k \times x_{ki} \]

Essentially, \(\mu\) is the model’s predicted average of \(y\) for a given (set of) \(x\).

\[ y_i \sim \mathcal{N}(\mu_i, \sigma^2), \quad i \in \{1, \dots, N\} \] \(\in\) means “is an element of”

• As useful as this model is, it is clearly limited
• What (sort of data) is the normal linear model limited to?

• Data are continuous (rather than discrete, categorical).
• Data are assumed to come from a normal distribution.
• Cannot be applied to categorical outcomes or counts.
• Generalised linear models allow us to use other distributions with link functions to map different types of outcome variables to a continuous numeric space.
• For example, binary logistic regression – a specific type of logistic regressions – uses a Binomial / Bernoulli distribution with a logit link to model binary outcomes (yes / no, correct / incorrect).

The normal distribution isn’t appropriate because the data are not continuous.

\[ \begin{aligned} y_i &\sim \mathcal{N}(\mu, \sigma^2) \end{aligned} \] We can use another distribution that can handle binary outcomes.

\[ \begin{aligned} y_i \sim \text{Bernoulli}(\theta) \\ \end{aligned} \]

Normal linear model

# lm =  linear model
lm(outcome ~ predictor_1 + predictor_2, data = data)

Binary logistic regression

# glm = generalised linear model
glm(outcome ~ predictor_1 + predictor_2, data = data, family = binomial(link = "logit"))

binomial(link = "logit") models data as coming from a Binomial / Bernoulli distribution where the logit links maps probabilities to a continuous numeric space.

\(\theta\) controls the relative heights of each bar

# load packages
library(tidyverse)

# load smoking data (link can be found on NOW)
smoking_df <- read_csv("https://raw.githubusercontent.com/mark-andrews/ntupsychology-data/main/data-sets/smoking.csv") 

# select variables of interest
smoking_df <- select( smoking_df, cigs, age )

# look at the data
glimpse( smoking_df )

Rows: 807
Columns: 2
$ cigs <dbl> 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 10, 20, 30, 0, 0, 20, 0, 30, 0, 20,…
$ age  <dbl> 46, 40, 58, 30, 17, 86, 35, 48, 48, 31, 27, 24, 71, 44, 22, 29, 3…

# Make the cigs a binary variable
smoking_df <- mutate(smoking_df, is_smoker = cigs > 0)

# How many smokers and nonsmokers are in the data?
count( smoking_df, is_smoker )

# A tibble: 2 × 2
  is_smoker     n
  <lgl>     <int>
1 FALSE       497
2 TRUE        310

# Probability of smokers (same as `mean(smoking_df$is_smoker)`)
310 / ( 310 + 497 )

[1] 0.3841388

# Overall probability of smokers in odds
0.3841388 / ( 1 - 0.3841388 )

[1] 0.6237425

# Overall probability of smokers in log-odds (logits)
log( 0.6237425 )

[1] -0.4720177

\(\dots\) transforms probability to log-odds.

Log-odds and odds are another way of talking about probability (e.g., kilometers and miles are both measures of distance, but one might make more sense to you than the other).

Importantly,

• Log-odds transformation ensures that \(\theta\) (probability of outcome) stays between 0 and 1 while log-odds are unbound.
• \(\theta\) is not a linear function of predictors but the log odds transformation of \(\theta\) is. \[ \begin{aligned} \log\left(\frac{\theta}{1 - \theta}\right) = \beta_0 + \beta_1 \times x_i. \end{aligned} \]

\[ \log\!\left(\frac{\theta}{1-\theta}\right) \]

# probability to logit
logit <- function( x ) log( x / ( 1 - x ) )

logit( .6 ) # prob to logit

[1] 0.4054651

# or R function
qlogis( .6 )

[1] 0.4054651

\[ \frac{1}{1+e^{-\theta}} \]

# logit to probability
ilogit <- function( x ) 1 / ( 1 + exp( -x ) )

ilogit( 0.4054651 )

[1] 0.6

# or R function
plogis( 0.4054651 )

[1] 0.6

If the probability of bad weather is 0.75, what are the odds that we’ll have bad weather?

\[ \text{odds} = \frac{\theta}{1 - \theta} \] Use \(R\) as a calculator or solve it by hand!

If the probability of bad weather is \(\theta = 0.75\), what are the odds that we’ll have bad weather?

p <- 0.75
p / ( 1 - p ) 

[1] 3

Probability of \(\theta = 0.75\) means that we are three times more likely to have bad weather than good weather.

If the odds are 4 (i.e. 4:1) to have bad weather, what is the probability that we will have bad weather?

\[ \theta = \frac{\text{odds}}{1 + \text{odds}} \] Use \(R\) as a calculator!

If the odds are 4 (i.e. 4:1) to have bad weather, what is the probability that we will have bad weather?

odds <- 4
odds / ( 1 + odds )

[1] 0.8

The log-odds on the previous slide were 2.1972246. Convert this value into

1. Odds (the opposite of log(x) is exp(x)).
2. Probabilities (using plogis(x)).

(replace x correctly)

The log-odds on the previous slide were 2.1972246. Convert this value into

Odds (the opposite of log(x) is exp(x)):

exp(2.197225)

[1] 9.000004

Probabilities (using plogis(x)):

plogis(2.197225)

[1] 0.9

What are the log-odds of being a smoker at any given age?

What are the odds of being a smoker at any given age?

What is the probability of being a smoker at any given age?

# Extract model coefficients
b <- coef(m_smoke)
b

(Intercept)         age 
 0.02396177 -0.01215431 

But what does log-odds of -0.279 actually mean?

What is the probability of being a smoker at 25 years of age?

The inverse logit function, allows us to convert from log-odds to probability.

\[ \theta = \frac{1}{1 + e^{-\text{logit}(\theta)}} \] We implemented ilogit(x) before, or use the base-R function plogis(x).

\[ \log\left(\frac{\text{is_smoker}}{1 - \text{is_smoker}}\right) = \beta_0 + \beta_1 \times 25 \]

which is

logit_smoke_age25 <- b[1] + b[2] * 25 # linear function of age = 25

plogis( logit_smoke_age25 ) # using plogis to convert log odds to prob

[1] 0.4304793

1 / ( 1 + exp( -logit_smoke_age25 ) ) # using the formula we implemted in ilogit(x)

[1] 0.4304793

What are the log-odds of being a smoker at 50 years of age?

What are the odds of being a smoker at 50 years of age?

What is the probability of being a smoker at 50 years of age?

# Step 1. Make a new data frame with range of values of age.
smoking_df2 <- tibble(age = seq(from = 15, to = 80, by = 10))

# Step 2. Use model from earlier to generate predictions.
predict(m_smoke, newdata = smoking_df2) # in log-odds

         1          2          3          4          5          6          7 
-0.1583529 -0.2798960 -0.4014392 -0.5229823 -0.6445254 -0.7660686 -0.8876117 

Specify type = 'response' to get predictions in probability.

# Base R function for vector of predictions
predict(m_smoke, newdata = smoking_df2, type = 'response')

# Alternative using `modelr` for output as data frame
library(modelr)
add_predictions(data = smoking_df2, model = m_smoke, type = 'response')