More Logistic Regression

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
3. Understand how to evaluate model fit

What You Need

1. Be up to date with lectures
2. Have completed previous lab exercises from Week 6

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.

Lab Data

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

Study Overview

Research Question

Is attempting to quit tobacco products associated with an individuals intentions?

Smoking Codebook

Description

Kalkhoran and colleagues (2015)[^1] investigated predictors of dual use of cigarettes and smokeless tobacco products. Although their original sample size was large (\(n\) = 1324), they were interested in running separate logistic regression analyses within subgroups. Most of the sample used only cigarettes. One of the smaller subgroups contained subjects who used both cigarettes and smokeless tobacco products (\(n\)=61). For this problem, we focus on this smaller subgroup of dual cigarette and smokeless tobacco users.

Data Dictionary

variable	description
S	Subject
Q	Quit Attempt - Whether the subject made an attempt to quit using tobacco products (0 = no attempt to quit; 1 = attempted to quit).
I	Intention to Quit - Individuals intention to quit (1 = Never intend to quit; 2 = May intend to quit but not in the next 6 months; 3 = Intend to quit in the next 6 months; 4 = Intend to quit in the next 30 days)

Preview

The first six rows of the data are:

S	Q	I
1	0	1
2	1	3
3	1	2
4	0	1
5	1	2
6	1	2

Setup

Setup

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

Solution

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

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

Exercises

Study & Analysis Plan Overview

Question 1

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

Solution

#look at structure of data:
str(smoke)

spc_tbl_ [62 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ S: num [1:62] 1 2 3 4 5 6 7 8 9 10 ...
 $ Q: num [1:62] 0 1 1 0 1 1 1 0 1 0 ...
 $ I: num [1:62] 1 3 2 1 2 2 3 3 3 2 ...
 - attr(*, "spec")=
  .. cols(
  ..   S = col_double(),
  ..   Q = col_double(),
  ..   I = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 

#check for NAs - there are none - all FALSE:
table(is.na(smoke))

FALSE 
  186 

#re-name variables to improve clarity:
smoke <- smoke %>%
    rename(PID = S,
           Quit_Attempt = Q,
           Intention = I)

#re-assign categorical IVs as factors, and give more appropriate labels to each level:
smoke$Intention <- factor(smoke$Intention,
                              levels = c(1,2,3,4),
                              labels = c("Never", "Maybe (not in next 6 months)", "Yes (within next 6 months)", "Yes (within next 30 days)"))

Descriptive Statistics & Visualisations

Question 2

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, since we have two categorical variables, the select() and table() functions will come in handy here.
2. Recall that when visualising a continuous outcome across several groups, geom_boxplot() may be most appropriate to use.
3. For your visualisations, you will need to specify as_factor() when plotting the Quit_Attempt variable since this is numeric, but we want it to be treated as a factor only for plotting purposes
4. Make sure to comment on any observed differences among the sample means of the four treatment conditions.

Solution

Numeric

Let’s first produce a descriptive statistics table:

smoke_stats <- smoke %>%
    select(Intention, Quit_Attempt) %>%
    table() %>%
    kable(caption = "Descriptive Statistics") %>%
    kable_styling()

smoke_stats

Table 1: Descriptive Statistics

	0	1
Never	7	1
Maybe (not in next 6 months)	9	17
Yes (within next 6 months)	5	18
Yes (within next 30 days)	1	4

Or alternatively we could view as proportions:

smoke_prop <- smoke %>%
  group_by(Intention, Quit_Attempt) %>%
  summarise(
    n = n()) %>%
  mutate(
    Prop = round(n/sum(n),2)
  ) %>%
    kable(caption = "Proportions") %>%
    kable_styling()

smoke_prop 

Table 2: Proportions

Intention	Quit_Attempt	n	Prop
Never	0	7	0.88
Never	1	1	0.12
Maybe (not in next 6 months)	0	9	0.35
Maybe (not in next 6 months)	1	17	0.65
Yes (within next 6 months)	0	5	0.22
Yes (within next 6 months)	1	18	0.78
Yes (within next 30 days)	0	1	0.20
Yes (within next 30 days)	1	4	0.80

Based on the proportions presented in Table 2, we can see for those who had some intention to quit (i.e., maybe - not in next 6 months, yes - within next 6 months, or yes - within 30 days), there was a larger proportion of individuals who attempted to quit using tobacco products.

Visual

We can visually explore the association between the two categorical variables as follows:

smoke_plt1 <- ggplot(data = smoke, aes(x=as_factor(Quit_Attempt), fill=Intention)) +
  geom_bar(position = "dodge") +
  labs(fill = 'Intention to Quit', x = "Attempted Quitting (0 = No, 1 = Yes)")
smoke_plt1

Figure 1: Association between Intention and Quitting Attempt

From Figure 1, we can see that for those who made an attempt to quit using tobacco products, the majority intended to quit - either in the immediate or distant future. For those who made no attempt to quit using tobacco products, the majority either had no intention of quitting or did not plan to quit in the near future.

Model Fitting & Interpretation

Question 3

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

Solution

smoke_mdl1 <- glm(Quit_Attempt ~ Intention, data = smoke, family = "binomial")
summary(smoke_mdl1)

Call:
glm(formula = Quit_Attempt ~ Intention, family = "binomial", 
    data = smoke)

Coefficients:
                                      Estimate Std. Error z value Pr(>|z|)   
(Intercept)                             -1.946      1.069  -1.820  0.06872 . 
IntentionMaybe (not in next 6 months)    2.582      1.146   2.253  0.02423 * 
IntentionYes (within next 6 months)      3.227      1.183   2.729  0.00636 **
IntentionYes (within next 30 days)       3.332      1.547   2.154  0.03123 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 80.648  on 61  degrees of freedom
Residual deviance: 68.659  on 58  degrees of freedom
AIC: 76.659

Number of Fisher Scoring iterations: 4

Question 4

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 (see here for more details if you are interested), 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(smoke_mdl1))

                          (Intercept) IntentionMaybe (not in next 6 months) 
                            0.1428571                            13.2222222 
  IntentionYes (within next 6 months)    IntentionYes (within next 30 days) 
                           25.2000000                            28.0000000 

• \(\beta_0\) = (Intercept) = 0.14
  • For a participant who reported that they “never intend to quit”, the probability that they did attempt quitting was 0.14 times the probability that they did not.
• \(\beta_1\) = Intention Maybe (not in next 6 months) = 13.22
  • Relative to those who “never intend to quit”, reporting an “intention to quit but not in the next 6 months” meant that an individuals odds of attempting to quit increased by a factor of 13.2.
• \(\beta_2\) = Intention Yes (within next 6 months) = 25.2
  • Relative to those who “never intend to quit”, reporting an “intention to quit in the next 6 months” increased the odds of attempting to quit by a factor of 25.2.
• \(\beta_3\) Intention Yes (within next 30 days) = 28
  • Relative to those who “never intend to quit”, reporting an “intention to quit in the next 30 days” was associated with an increase in the odds of attempting to quit by a factor of 28.

Model Fit

Question 5

Examine the below plot to determine if the deviance residuals raise concerns about outliers:

plot(rstandard(smoke_mdl1, type = 'deviance'), ylab = 'Standardised Deviance Residuals')

Based on this plot, are there any residuals of concern? Are there any additional plots you could check to determine if there are influential observations?

Hint

Deviance Residuals

Because logistic regression models don’t have the same expected error distribution (we don’t expect residuals to be normally distributed around the mean, with constant variance), checking the assumptions of logistic regression is a little different.

Typically, we look at the “deviance residuals”. But we don’t examine plots for patterns, we simply examine them for potentially outlying observations. If we use a standardised residual, it makes it easier to explore extreme values as we expect most residuals to be within -2, 2 or -3, 3 (depending on how strict we feel).

Cook’s Distance

To check for influential observations, we can use cooks.distance(), and plot this using plot(). Alternatively, you can specify which = 4 when plotting your fitted model. In logistic regression, we can use the arbitrary cut-offs of 0.5 (moderately influential) or 1 (highly influential) to describe influential points.

Solution

In the plot above, there appears to be 1 residual with a value >|2|, but none with a value >|3|. We will keep this in mind and check if they are also an influential point by plotting Cook’s Distance:

#cooks d - option 1 (only need to use one of these two options)
plot(cooks.distance(smoke_mdl1), ylab = "Cook's Distance")

#cooks d - option 2 (only need to use one of these two options)
plot(smoke_mdl1, which = 4)

There doesn’t appear to be any influential observations based on our Cook’s distance plot, since all values < .50.

Question 6

Perform a Deviance goodness-of-fit test to compare your fitted model to the null.

\[ \begin{aligned} M_0 &: \qquad \log \left( \frac{p}{1 - p}\right) = \beta_0 \\ M_1 &: \qquad \log \left( \frac{p}{1 - p}\right) = \beta_0 + \beta_1 Intention2 + \beta_2 Intention3 + \beta_3 Intention4 \end{aligned} \]

Report which model you think best fits the data.

Hint

Consider whether or not your models are nested. The flowchart in Question 10 of the Semester 2 Week 1 lab may be helpful to revisit.

Solution

First, let’s fit the null model:

smoke_mdl0 <- glm(Quit_Attempt ~ 1, data = smoke, family = "binomial")
summary(smoke_mdl0)

Call:
glm(formula = Quit_Attempt ~ 1, family = "binomial", data = smoke)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.5978     0.2654   2.252   0.0243 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 80.648  on 61  degrees of freedom
Residual deviance: 80.648  on 61  degrees of freedom
AIC: 82.648

Number of Fisher Scoring iterations: 4

Since our models are nested, we can compare using the likelihood ratio test:

anova(smoke_mdl0, smoke_mdl1, test = 'Chisq')

Analysis of Deviance Table

Model 1: Quit_Attempt ~ 1
Model 2: Quit_Attempt ~ Intention
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
1        61     80.648                        
2        58     68.659  3   11.989  0.00742 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

At the 5% significance level, the addition of information about the participants’ intention to quit resulted in a significant decrease in model deviance \(\chi^2(3) = 11.99, p = .007\).

Hence, we have strong evidence that the model the subjects’ intention to quit is a helpful predictor of whether or not they will attempt quitting in the future.

Question 7

Check the AIC and BIC values for smoke_mdl0 and smoke_mdl1 - which model should we prefer?

Solution

#AIC
AIC(smoke_mdl0, smoke_mdl1)

           df      AIC
smoke_mdl0  1 82.64844
smoke_mdl1  4 76.65904

#BIC
BIC(smoke_mdl0, smoke_mdl1)

           df      BIC
smoke_mdl0  1 84.77557
smoke_mdl1  4 85.16758

We used AIC and BIC model selection to distinguish between two possible models describing the association between attempting to quit and intentions. Our model with with intentions included as a predictor (AIC = 76.66) was better fitting than the null model (AIC = 82.65). However, the BIC values suggested that the model including intentions (BIC = 85.17) was a poorer fit than the null (BIC = 84.77). Based on the weight of evidence from both the Deviance goodness-of-fit test alongside the AIC and BIC values, we would conclude that the model with intentions was better fitting than the null.

Visualise Model

Question 8

Plot the following:

• predicted model estimates
• predicted probability

Hint

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

Solution

Predicted model estimates

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

plot_model(smoke_mdl1,
           type = "est")

Figure 2: Model Estimates

Predicted probability

plot_model() with type = "eff" allows us to visualise the predicted probability based on our model estimates:

plot_model(smoke_mdl1,
           type = "eff")

$Intention

Figure 3: Predicted Probability

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(smoke_mdl1,
          dv.labels = "Attempt to quit using tobacco products",
          pred.labels = c("Maybe (not in next 6 months)" = "Intend to quit but not in the next 6 months",
                          "Yes (within next 6 months)" = "Intend to quit in the next 6 months",
                          "Yes (within next 30 days)" = "Intend to quit in the next 30 days"),
          title = "Regression table for Smoking Model")

Table 3: Regression table for Smoking Model

	Attempt to quit using tobacco products
Predictors	Odds Ratios	CI	p
(Intercept)	0.14	0.01 – 0.80	0.069
IntentionMaybe (not in next 6 months)	13.22	1.93 – 268.17	0.024
IntentionYes (within next 6 months)	25.20	3.42 – 534.12	0.006
IntentionYes (within next 30 days)	28.00	1.95 – 1181.69	0.031
Observations	62
R2 Tjur	0.192

Question 10

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

Solution

Whether or not participants made an attempt to quit using tobacco products (binary 0 vs 1; 0 no attempt to quit, 1 attempted to quit) was modeled using logistic regression, with intention to quit (‘never intend to quit’, ‘may intend to quit but not in the next 6 months’, ‘intend to quit in the next 6 months’, ‘intend to quit in the next 30 days’, with ‘never intend to quit’ as the reference level) as the only predictor. See Table 3 for full model results, and Figure 2 for a visualisation of the model estimates and confidence intervals.

Having no intention to quit was found to be associated with attempting to quit (\(OR = 0.14,\,\, 95\%\, CI\, [0.01, 0.8]\)), where the probability that they did attempt to quit was 0.14 times the probability that they did not.

Relative to those who never intend to quit, individuals who may intend to quit but not in the next six months had increased odds of attempting to quit (\(OR = 13.22,\,\, 95\%\, CI\, [1.93, 268.17]\)).

In comparison to those with no intention to quit, those who intended to quit in the next six months had their the odds of quitting increased by a factor of 25.2 (\(95\%\, CI\, [3.42, 534.12]\))

Finally, intending to quit in the next 30 days was associated with an increase in the odds of attempting to quit by a factor of 28 (\(95\%\, CI\, [1.95, 1181.69]\)) in comparison to those who never intended to quit.

In summary, the closer in the future someone intended to quit smoking, the higher their odds of attempting to quit using tobacco products (see Figure 3).

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.