Logistic Regression I

Learning Objectives

At the end of this lab, you will:

  1. Understand when to use a logistic model
  2. Understand how to fit and interpret a logistic model

What You Need

  1. Be up to date with lectures

Required R Packages

Remember to load all packages within a code chunk at the start of your RMarkdown file using library(). If you do not have a package and need to install, do so within the console using install.packages(" "). For further guidance on installing/updating packages, see Section C here.

For this lab, you will need to load the following package(s):

  • tidyverse
  • patchwork
  • kableExtra
  • psych
  • sjPlot

Presenting Results

All results should be presented following APA guidelines.If you need a reminder on how to hide code, format tables/plots, etc., make sure to review the rmd bootcamp.

The example write-up sections included as part of the solutions are not perfect - they instead should give you a good example of what information you should include and how to structure this. Note that you must not copy any of the write-ups included below for future reports - if you do, you will be committing plagiarism, and this type of academic misconduct is taken very seriously by the University. You can find out more here.

Lab Data

You can download the data required for this lab here or read it in via this link https://uoepsy.github.io/data/drunkdoor.csv

Study Overview

Research Question

Is susceptibility to change blindness influenced by age, level of alcohol intoxication, and perceptual load?

Drunk door codebook

Setup

Setup
  1. Create a new RMarkdown file
  2. Load the required package(s)
  3. Read the drunkdoor dataset into R, assigning it to an object named drunkdoor 

#load packages
library(tidyverse)
library(psych)
library(kableExtra)
library(patchwork)
library(sjPlot)

#read in data
drunkdoor <- read_csv("https://uoepsy.github.io/data/drunkdoor.csv")

Exercises

Study & Analysis Plan Overview

Question 1

Firstly, examine the dataset, and perform any necessary and appropriate data management steps.

Next, think about the scales that the variables are currently on, with a particular focus on BAC, age, and condition. Consider:

  • Do you want BAC on the current scale, or could you transform it somehow? Consider that instead of the coefficient representing the difference when moving from 0% to 1% BAC (1% blood alcohol is fatal!), we might want to have the difference associated with 0% to 0.01% BAC (i.e, a we want to talk about effects in terms of changing 1/100th of a percentage of BAC)
  • Do you want age to be centred at 0 years (as it currently is), or could you re-centre to make it more meaningful?
  • In your data management, you will hopefully make condition a factor, but have you considered the reference level? It would likely make most sense for this to be set as “Low”.

  • The str() function will return the overall structure of the dataset, this can be quite handy to look at.

  • When considering the scale of your BAC and age variables, review the data transformation flashcards.

Data Management - Convert categorical variables to factors, and if needed, provide better variable names - Label factors appropriately to aid with your model interpretations if required
- Check that the dataset is complete (i.e., are there any NA values?). We can check this using is.na()

Note that all of these steps can be done in combination - the mutate() and factor() functions will likely be useful here.

Let’s examine the data:

#look at structure of data
str(drunkdoor)
spc_tbl_ [120 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ id       : chr [1:120] "ID1" "ID2" "ID3" "ID4" ...
 $ bac      : num [1:120] 0.06744 0.00331 0.00323 0.07984 0.06676 ...
 $ age      : num [1:120] 43 64 44 67 62 45 50 45 48 61 ...
 $ condition: chr [1:120] "Low" "Low" "Low" "High" ...
 $ notice   : num [1:120] 0 0 1 0 0 1 0 1 0 0 ...
 - attr(*, "spec")=
  .. cols(
  ..   id = col_character(),
  ..   bac = col_double(),
  ..   age = col_double(),
  ..   condition = col_character(),
  ..   notice = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 
#check for NAs - there are none - all FALSE
table(is.na(drunkdoor))

FALSE 
  600

Next we need to specify Condition as a factor, and set ‘Low’ as the reference group:

#re-assign categorical IVs as factors
drunkdoor$condition <- factor(drunkdoor$condition,
                              levels = c("Low", "High"))
#set 'low' as reference group
drunkdoor$condition <- fct_relevel(drunkdoor$condition, "Low")

It may be more useful to have blood alcohol (BAC) in terms of 100ths of percentages, rather than percentages, for the reasons noted above. It would also be more useful for interpretation to have age centered on the mean (note that this won’t change anything other than the intercept in our model). Let’s make both these changes:

drunkdoor <- drunkdoor %>% 
  mutate(
    bac100 = bac*100,
    age_mc = age - mean(age)
  )


Question 2

Provide a brief overview of the study design and data, before detailing your analysis plan to address the research question.

  • Give the reader some background on the context of the study
  • State what type of analysis you will conduct in order to address the research question
  • Specify the model to be fitted to address the research question (note that you will need to specify the reference level of your categorical variable)
  • Specify your chosen significance (\(\alpha\)) level
  • State your hypotheses

Much of the information required can be found in the Study Overview codebook.

The drunkdoor dataset contained information on 120 hypothetical participants who were approached within the vicinity of a number of establishments with licenses to sell alcohol to be consumed on-premises. Using a between-subjects design experimental design, experimenters manipulated perceptual load in two ways (door appearance and number of papers held by experimenters) to create low- and high-load conditions. In the low-load condition, the door was a standard MDF construction painted in a neutral grey, and the experimenters handled only 2 sheets of paper which had minimal printed text. In the high-load condition, the door was painted with detailed graffiti and had numerous pieces of paper and notices attached to the side facing the participants, and the experimenters held a pile of 30 papers containing drawings and printed text. Experimenters recorded whether or not participants noticed the mid-conversation experimenter swap (no = 0, yes = 1). The researchers also collected information on participants’ blood alcohol level (BAC; percentage of alcohol in an individual’s blood stream) and age (in years).

To aid with interpretation, age was mean centered. The BAC measure was also transformed to represent a one-unit change as 0.01%, as opposed to 1%. Low-load was specified as the reference group of the perceptual load condition.

Density plots and bar plots will be used to individually visualise the associations between whether the mid-conversation swap was noticed and BAC, age, and condition.

To investigate whether the probability of noticing the swap differed as a function of age, BAC, and perceptual load, a binary logistic regression model was used, where the following model was specified:

\[ \begin{aligned} \ln \left(\frac{p}{1 - p}\right) ~=~ \beta_0 + \beta_1 \cdot Age + \beta_2 \cdot BAC + \beta_3 \cdot Condition_{High} \end{aligned} \]

\[ \begin{aligned} \text{where}~{p}~ &=~ \text{probability of noticing the mid-conversation swap} \end{aligned} \]

Our hypotheses are:

\(H_0:\) All \(\beta_j = 0\) (for \(j = 1, 2, 3\))

Susceptibility to change blindness is not influenced by age, level of alcohol intoxication, or perceptual load.

\(H_1:\) At least one \(\beta_j \neq 0\) (for \(j = 1, 2, 3\))

Susceptibility to change blindness is influenced by age, level of alcohol intoxication, and/or perceptual load.

Effects will be considered statistically significant at \(\alpha = .05\).

Descriptive Statistics & Visualisations

Question 3

Provide a table of descriptive statistics and visualise your data.

Remember to interpret these in the context of the study.

Review the many ways to numerically and visually explore your data by reading over the data exploration flashcards.

For examples, see flashcards on descriptives statistics tables - categorical and numeric values examples and categorical x categorical example - visualise data.

Specifics to consider:
1. For your table of descriptive statistics, both the group_by() and summarise() functions will come in handy here 2. For your visualisations, you will need to specify as_factor() when plotting the notice variable since this is numeric, but we want it to be treated as a factor only for plotting purposes

Let’s first produce a descriptive statistics table, grouped by whether or not participants noticed the swap, and what perceptual load condition they were assigned to:

door_stats <- drunkdoor %>% 
    group_by(notice, condition) %>%
    summarise(
        M_Age = mean(age), 
        SD_Age = sd(age),
        SE_Age = sd(age) / sqrt(n()),
        Min_Age = min(age),
        Max_Age = max(age),        
        M_BAC = mean(bac), 
        SD_BAC = sd(bac),
        SE_BAC = sd(bac) / sqrt(n()),
        Min_BAC = min(bac),
        Max_BAC = max(bac)) %>%
    kable(caption = "Descriptive Statistics", digits = 2) %>%
    kable_styling()

door_stats
Table 1: Descriptive Statistics
notice condition M_Age SD_Age SE_Age Min_Age Max_Age M_BAC SD_BAC SE_BAC Min_BAC Max_BAC
0 Low 61.62 8.09 1.65 41 74 0.03 0.02 0 0.00 0.07
0 High 60.21 7.96 1.39 46 79 0.04 0.02 0 0.00 0.08
1 Low 52.42 6.78 1.13 41 65 0.04 0.03 0 0.00 0.08
1 High 47.93 7.05 1.36 34 59 0.05 0.02 0 0.01 0.08 
door_plt1 <- ggplot(data = drunkdoor, aes(x=as_factor(notice), fill=as_factor(notice))) +
  geom_bar() +
  labs(x = "Noticed Swap (0 = No, 1 = Yes)", fill = "Noticed Swap (0 = No, 1 = Yes)", y = "Frequency")
door_plt1

Figure 1: Counts of Notice

From Figure 1, we can see that there are slightly more individuals that noticed the mid-conversation swap than those who did not.

door_plt2 <- ggplot(data = drunkdoor, aes(x=age, fill=as_factor(notice))) +
  geom_density() +
  labs(x = "Age (in years)", fill = "Noticed Swap (0 = No, 1 = Yes)")
door_plt2

Figure 2: Association between Notice and Age

From Figure 2, we can see that there appear to be more older adults who did not notice the mid-conversation swap.

door_plt3 <- ggplot(data = drunkdoor, aes(x=bac, fill=as_factor(notice))) +
  geom_density() +
  labs(x = "BAC", fill = "Noticed Swap (0 = No, 1 = Yes)")
door_plt3

Figure 3: Association between Notice and BAC

From Figure 3, we can see that higher BAC appears to be positively associated with noticing a mid-conversation swap.

door_plt4 <- ggplot(data = drunkdoor, aes(x=condition, fill=as_factor(notice))) +
  geom_bar(position = "dodge") +
  labs(fill = "Noticed Swap (0 = No, 1 = Yes)", x = "Condition")
door_plt4

Figure 4: Association between Notice and Perceptual Load Condition

From Figure 4, in the Low perceptual load condition, more participants noticed the swap than did not, whilst the opposite pattern emerged in the High perceptual group.


Question 4

For the moment, lets just consider the association between notice and age. Visually following the line from the plot produced below, what do you think the predicted model value would be for someone who is aged 30? What does this value mean?

Hopefully you can see that the model predicted value for someone aged 30 is approximately 1.30.

But let’s pause to contemplate what this means, keeping in mind the way in which the ‘Notice’ variable was measured (i.e., no = 0, yes = 1).

Does this mean then that a 30 year old participant will… Notice 1.30 experimenter-swaps? Have a 130% probability of noticing the swap?

Neither of these interpretations are quite right - participants could only notice one swap, and probability is between 0 and 1.

Therefore, it seems that we cannot capture this as a linear relationship - we instead need a logistic regression model.

Model Fitting & Interpretation

Question 5

Fit your model using glm(), and assign it as an object with the name “changeblind_mdl”.

For an overview of how to fit a binary logistic regression model using the glm() function, see the binary logistic regression flashcard.

For an example, review the binary logistic regression - example flashcard.

changeblind_mdl <- glm(notice ~ age_mc + bac100 + condition, data = drunkdoor, family = "binomial")
summary(changeblind_mdl)

Call:
glm(formula = notice ~ age_mc + bac100 + condition, family = "binomial", 
    data = drunkdoor)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -0.51577    0.53536  -0.963  0.33534    
age_mc        -0.22584    0.04214  -5.359 8.35e-08 ***
bac100         0.35002    0.11128   3.145  0.00166 ** 
conditionHigh -1.59215    0.54008  -2.948  0.00320 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 166.06  on 119  degrees of freedom
Residual deviance: 100.99  on 116  degrees of freedom
AIC: 108.99

Number of Fisher Scoring iterations: 5


Question 6

Conduct a model comparison of your model above against the null model. Report the results of the model comparison in APA format.

Consider whether or not your models are nested. The flowchart may be helpful to revisit, as well as the model comparisons - logistic regression flashcards

#fit null model
changeblind_mdl_null <- glm(notice ~ 1, data = drunkdoor, family = "binomial")

#run model comparison
anova(changeblind_mdl_null, changeblind_mdl, test= "Chisq")
Analysis of Deviance Table

Model 1: notice ~ 1
Model 2: notice ~ age_mc + bac100 + condition
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       119     166.06                          
2       116     100.99  3   65.065 4.857e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

At the 5% significance level, the addition of information about the participants’ age, BAC, and perceptual load resulted in a significant decrease in model deviance \(\chi^2(3) = 65.07, p < .001\).

Hence, we have strong evidence that the participants’ age, BAC, and perceptual load condition are helpful predictors of whether or not they were likely to notice a mid-conversation swap.


Question 7

Interpret your coefficients in the context of the study. When doing so, it may be useful to translate the log-odds back into odds.

For an overview of how probability, odds, and log-odds are related, review the probability, odds, and log-odds flashcard.

If you need help interpreting coefficients, review the interpretation of coefficients flashcard.

For an example, review the binary logistic regression - example flashcard.

exp(coefficients(changeblind_mdl))
  (Intercept)        age_mc        bac100 conditionHigh 
    0.5970388     0.7978495     1.4190910     0.2034874

\(\beta_0\) = (Intercept) = 0.6

  • The odds of noticing a mid-conversation person-switch for a sober person of average age, with Low perceptual load was 0.6.

\(\beta_1\) = age_mc = 0.8

  • For every year older someone was, the odds of noticing a mid-conversation swap were multiplied by 0.8, after accounting for differences in blood alcohol levels and perceptual load.

\(\beta_2\) = bac100 = 1.42

  • For every 1/100th of a percentage increase in BAC, the odds of noticing a mid-conversation person switch were multiplied by 1.42, after accounting for differences in age and perceptual load.

\(\beta_3\) = conditionHigh = 0.2

  • The odds of noticing the swap were significantly different depending upon the perceptual load. For those in the High perceptual condition, the odds of noticing a change decreased by a factor of 0.2 compared to when the condition was Low, after accounting for differences in age and BAC.

Visualise Model

Question 8

Plot the predicted model estimates.

Here you will need to use plot_model() from the sjPlot package. To get your estimates, you will need to specify type = "est".

plot_model() with type = "est" gives a nice way of visualising the model odds ratios and confidence intervals:

plot_model(changeblind_mdl,
           type = "est")

Figure 5: Model Estimates

Writing Up & Presenting Results

Question 9

Provide key model results in a formatted table.

Use tab_model() from the sjPlot package. For a quick guide, review the tables flashcard.

#create table for results
tab_model(changeblind_mdl,
          dv.labels = "Notice",
          pred.labels = c("age_mc" = "Age (in years)",
                          "bac100" = "Blood Alcohol Level (BAC)",
                          "conditionHigh" = "Condition - High"),
          title = "Regression table for Notice Model")
Table 2: Regression table for Notice model
  Notice
Predictors Odds Ratios CI p
(Intercept) 0.60 0.20 – 1.69 0.335
Age (in years) 0.80 0.73 – 0.86 <0.001
Blood Alcohol Level (BAC) 1.42 1.15 – 1.79 0.002
Condition - High 0.20 0.07 – 0.56 0.003
Observations 120
R2 Tjur 0.475


Question 10

Interpret your results in the context of the research question and report your model in full.

Make reference to the regression table.

For an example of coefficient interpretation, review the interpretation of coefficients flashcard.

Whether or not participants noticed the swap mid-conversation (binary 0 vs 1) was modeled using logistic regression, with blood alcohol content (BAC; measured in 100th of percentages blood content), perceptual load condition (low load vs high load, with low as the reference level), and age (in years). See Table 2 for full model results, and Figure 5 for a visualisation of the model estimates and confidence intervals.

Age was found to be associated with susceptibility to change-blindness, where for every one year increase in age, the odds of a person noticing the mid-conversation swap decreased, with \(odds = 0.8\) (\(95\%\, CI\, [0.73, 0.86]\)), after accounting for differences in blood alcohol levels and perceptual load.

In contrary to what might have been expected, change-blindness appeared to decrease with higher alcohol intoxication, where for every 1/100th of a percentage increase in blood alcohol content, there were 1.42 times (\(95\%\, CI\,[1.15, 1.79]\)) increased odds of noticing the swap when holding age and perceptual load constant.

After accounting for age and blood alcohol levels, the odds of noticing the swap were significantly different depending upon the perceptual load. For those with a high perceptual load, the odds of noticing a swap decreased by a factor of 0.2 (\(95\%\, CI\,[0.07, 0.56]\)).

In summary, older age increased the susceptibility to change blindness, high perceptual load was found to decrease the odds of people being blind to change, and increased levels of alcohol intoxication appeared to be associated with greater odds of noticing change.

Compile Report

Compile Report

Knit your report to PDF, and check over your work. To do so, you should make sure:

  • Only the output you want your reader to see is visible (e.g., do you want to hide your code?)
  • Check that the tinytex package is installed
  • Ensure that the ‘yaml’ (bit at the very top of your document) looks something like this:
---
title: "this is my report title"
author: "B1234506"
date: "07/09/2024"
output: bookdown::pdf_document2
---

If you are having issues knitting directly to PDF, try the following:

  • Knit to HTML file
  • Open your HTML in a web-browser (e.g. Chrome, Firefox)
  • Print to PDF (Ctrl+P, then choose to save to PDF)
  • Open file to check formatting

To not show the code of an R code chunk, and only show the output, write:

```{r, echo=FALSE}
# code goes here
```

To show the code of an R code chunk, but hide the output, write:

```{r, results='hide'}
# code goes here
```

To hide both code and output of an R code chunk, write:

```{r, include=FALSE}
# code goes here
```

You must make sure you have tinytex installed in R so that you can “Knit” your Rmd document to a PDF file:

install.packages("tinytex")
tinytex::install_tinytex()

You should end up with a PDF file. If you have followed the above instructions and still have issues with knitting, speak with a Tutor.