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

Description

For an example of a change blindness experiment, watch the following video:

Or review this paper: Simons, D. J., & Levin, D. T. (1997). Change blindness. Trends in cognitive sciences, 1(7), 261-267.. You may well have already heard of these series of experiments, or have seen similar things on Netflix.

Method

Researchers conducted a study in which they approached 120 people, recruited from within the vicinity of a number of establishments with licenses to sell alcohol to be consumed on-premises. Initially, experimenter A approached participants and asked if they were interested in participating in a short study, and obtained their written consent. While experimenter A subsequently talked each participant through a set of questions on multiple pieces of paper (with the pretense of explaining what the participant was required to do), experimenters B and C carrying a door passed between the participant and experimenter A, with experimenter C replacing A (as can be viewed in the video).

The perceptual load of the experiment was manipulated via a) the presentation of the door and b) the papers held by the experimenters. For 60 of these participants, the door was painted with some detailed graffiti and had a variety of pieces of paper and notices attached to the side facing the participants. Additionally, for these participants, the experimenters handled a disorganised pile of 30 papers, with the top pages covered in drawings around the printed text. For the remaining 60, the door was a standard MDF construction painted a neutral grey, and the experimenters handled only 2 sheets of paper which had minimal printed text on them and nothing else.

Measures

After experimenters A and C had successfully swapped positions, the participant was asked (now by C) to complete small number of questions taking approximately 1 minute. Either after this set of questions, or if the participant made an indication that they had noticed the swap, the experimenters regrouped and the participant was explicitly asked whether they had noticed the swap.

Immediately after this, participants were breathalysed, and their blood alcohol content was recorded.

Data Dictionary

variable	description
id	Unique ID number
bac	Blood Alcohol Content (BAC) - A BAC of 0.0 is sober, while in the United States 0.08 is legally intoxicated, and above that is very impaired. BAC levels above 0.40 are potentially fatal.
age	Age (in years)
condition	Perceptual load created by distracting object (door) and details and amount of papers handled in front of participant. Levels = Low, High
notice	Whether or not the participant noticed the swap (No = 0 vs Yes = 1)

Preview

The first six rows of the data are:

id	bac	age	condition	notice
ID1	0.067437348	43	Low	0
ID2	0.003312794	64	Low	0
ID3	0.003233583	44	Low	1
ID4	0.079840000	67	High	0
ID5	0.066760207	62	High	0
ID6	0.015507668	45	Low	1

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

Solution

#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”.

Hint

• Convert categorical variables to factors
• Label appropriately factors to aid with your model interpretations if required
• If needed, provide better variable names

Solution

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.

Hint

• 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
• Specify your chosen significance (\(\alpha\)) level
• State your hypotheses

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

Solution

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.

Hint

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
3. Make sure to comment on your observations

Solution

Numeric

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

Visual

Notice

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.

Notice & Age

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.

Notice & BAC

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.

Notice & Task

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?

Solution

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”.

Solution

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.

Hint

Consider whether or not your models are nested. The flowchart in Semester 1 Week 4 Lab, Q10 may be helpful to revisit.

Solution

#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.

Hint

The opposite of the natural logarithm is the exponential, and in R these functions are log() and exp().

Recall that we can obtain our parameter estimates using various functions such as summary(),coef(), coefficients(), etc. Thus, we want to exponentiate the coefficients from our model in order to translate them back from log-odds.

Solution

exp(coefficients(changeblind_mdl))

  (Intercept)        age_mc        bac100 conditionHigh 
    0.5970388     0.7978495     1.4190910     0.2034874 

(Intercept)

\(\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.

age_mc

\(\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.

bac100

\(\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.

conditionHigh

\(\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.

Hint

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

Solution

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.

Hint

Use tab_model() from the sjPlot package.

Remember that you can rename your DV and IV labels by specifying dv.labels and pred.labels.

Solution

#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 your regression table.

Solution

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
---

What to do if you cannot knit to PDF

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

Hiding Code and/or Output

Hiding R Code

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

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

Hiding R Output

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

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

Hiding R Code and Output

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

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

Tinytex

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()

Solution

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.