Introduction to logistic regression

We have, by now, covered linear regression which is used when an outcome variable is continuous. For count outcomes, we use Poisson regression and its derivatives, where the outcome variable is a count or integer. For outcomes that are binary (“yes” or “no” for example), logistic regression is frequently the answer.

Getting started

We will be leaning on the following libraries to do these analyses:

• regplot
• PerformanceAnalytics
• jtools
• corrplot
• finalfit
• car
• arsenal

Background

Marijuana by Matthew Brodeur from unsplash.com

Imagine that we have been assigned a task to consider risk factors for marijuana use, and we have in mind some programmes to help reduce harmful use of the drug, but we need to know where to target this resource. How can we do that?

Well, we have a survey of marijuana use in a representative population, so we will analyse the data to see what the major risk behaviours are.

This is available in the NHANES US national health survey.

What sort of factors do you think are likely to increase the risk of taking recreational marijuana?

Summary of analysis

1. Load libraries and import data (here the data is imported from an installed library).

2. Check data integrity (we will skip through this)

3. Subset data. It doesn’t make sense to ask very young children about their marijuana use! In an epidemiological study, the denominator should consist of people at risk of the event of interest!

4. Plot associations between various demographic and behavioural risk factors and marijuana use. Conduct crude analyses, using odds ratios. We’ll be using histograms, boxplots (lattice package) and Euler plots (eulerr package). Check the bivariable associations and correlations between variables PerformanceAnalytics::chart.Correlation()). For this, we will need to convert character or factor variables ("Yes" or "No", or TRUE or FALSE) to numeric variables (1 or 0).

5. If writing a paper, you would normally make a ‘Table 1’ in which the various risk factors are divided by the outcome with column percentages and appropriate statistical tests. We will skip through this, but include it for completeness. This gives us similar information as in the previous point, but limited only to associations between exposures, confounders and outcome, rather than also between exposures and confounders as well. We will use the arsenal::tableby() function.

6. Logistic regression, use finalfit() and regplot()to finalise and visualise models. jtools is useful for making initial model output interpretable!

7. Check model fit (car::marginalModelPlots()) and adjust the model if necessary. The Hosmer-Lemeshow test is also relevant for checking the calibration between observed and predicted values.

8. Interpret the study results. Which variables are significantly associated Could there be other explanations other than a causal explanation? Could it be: reverse causation, measurement error, social desirability bias?

Loading libraries and data

Again, we use the pacman library to load other libraries.

if(!require(pacman)) install.packages("pacman")


pacman::p_load(magrittr, epiDisplay, tidyverse,
               rio, NHANES, visdat, skimr, finalfit, knitr, 
               kableExtra, readr, dplyr, lattice, jtools, 
               ggstance, broom.mixed, car, eulerr, regplot, 
               AF,# Model based attributable fractions 
               PerformanceAnalytics, # correlations    
              arsenal) # Making tables

## Respond with a "Y" if asked for a personal library.
update.packages(ask = FALSE)

Run the next line of code to increase the width of the output in the console.

## Setup work space
options(width = 600)

Load our data, which is the NHANES dataset.

## Load data----
data("NHANES")

Use ? to find out more about this dataset.

?NHANES

It doesn’t really make sense to ask children about marijuana use. So we’ll focus on those aged 20 or over. Run this line of code, which creates a new dataset called Na that includes only those aged ≥20 years.

## Restrict to adults
Na <- NHANES |> subset(Age >= 20) |> data.frame()

Univariable distributions and associations

In this section, we’ll look at the distribution of age and behaviour group, and see what groups are associated with marijuana use.
Use the epiDisplay::tab1() function to view the groups sizes for our outcome, Marijuana.

Na$Marijuana |> tab1()

Notice that there are a lot of missing responses (NAs)?
Remove these with the next line of code which creates a new data.frame called M:
- is.na() selects NA (missing) values and ! means “not”
- So, !is.na() findings selects non-missing values.

Restrict to those with complete responses to Marijuana variable

M <- Na |> subset(!is.na(Na$Marijuana))

Use the skim() function to get a summary of our dataset

M |> skim()

Check for any implausible values (check p0 and p100).

M$Marijuana |> tab1()

M$Depressed |> tab1()

M$Marijuana |> str()

cc(M$Marijuana, M$Depressed) # Outcome first, exposure second

# cc(dataset$outcome, dataset$exposure)

cc(M$Marijuana, M$Gender)

lattice::bwplot(Age ~ Marijuana, data = M)

# lattice::bwplot(variable1 ~ variable2, data = dataset)

lattice::bwplot(BMI ~ Marijuana, data = M)

## Convert SmokeNow from factor to character
M$SmokeNow <- M$SmokeNow |> as.character()

## Missing data for Smokenow indicate 'never smokers'
M$SmokeCurrent <- ifelse(is.na(M$SmokeNow), "Never", M$SmokeNow)

M$SmokeCurrent |> tab1()

tab1 <- arsenal::tableby(Marijuana ~ # column
                            SmokeCurrent,
                            cat.stats="countrowpct", #row percents are best for estimating risk ratios.
                            data = M)

summary(tab1, text = TRUE)

lattice::bwplot(AlcoholYear ~ Marijuana, data = M,
                ylab = "Average number of days in year alcohol consumed")

## These values are missing but should be rather zero
M[M$AlcoholYear == 0,]$AlcoholDay |> summary()

## Change them to zero rather than missing
M$AlcoholDay_corrected <- ifelse(M$AlcoholYear == 0, 0, M$AlcoholDay)

M$AlcoholDay |> summary()

lattice::bwplot(AlcoholDay_corrected ~ Marijuana, data = M,
                ylab = "Average number of alcoholic drinks consumed 
                in a drinking session")

car::powerTransform(I(AlcoholDay_corrected + 0.5) ~ 1, data = M,
                          family = "bcPower") |> summary()

car::symbox(~ I(AlcoholDay_corrected + 0.5), data = M)

summary(M$AlcoholDay)

lattice::bwplot(I(log(AlcoholDay_corrected + 0.5)) ~ Marijuana, data = M,
                ylab = "ln(Average number of alcoholic drinks consumed 
                in a drinking session)")

lattice::bwplot(HHIncomeMid ~ Marijuana, data = M,
                ylab = "Household income (US$ per year)")

cc(M$Marijuana, M$HardDrugs)

Euler plots

Euler plots are Venn diagrams that allow us to visualise associations between categorical variables, such as exposures and outcomes. These show the overlap between groups, which are represented as circles or ellipses.

When interpreting a Euler plot please note that:
- The area of the ellipse is proportional to the counts of individuals with that characteristic.
- More overlap between variables indicates stronger statistical association between them.

Create a Euler plot to visualise associations between Marijuana, HardDrugs and Alcohol12PlusYr.

M$All <- 1 # Assigns `1` to all rows. We’ll use this to create an ellipse for all participants

# Euler chooses the reference category as that to plot.
# We want to show the "Yes" category, so we'll change the reference
# from "No" to "Yes"

M$Marijuana_factor <- M$Marijuana |> factor() |> relevel(ref = "Yes")
M$HardDrugs_factor <- M$HardDrugs |> factor() |> relevel(ref = "Yes")
M$Alcohol_factor <- M$Alcohol12PlusYr |> factor() |> relevel(ref = "Yes")

# Euler plot
fit1 <- eulerr::euler(M[, c("Marijuana_factor",
                           "HardDrugs_factor",
                           "Alcohol_factor"  ,
                           "All")] |> na.omit(),# subset data and remove missing
                            shape = "ellipse")

plot(fit1, 
     main = "Marijuana, Alcohol and Hard drugs")

Note:
- 3 variables in the plot end in _Yes. This is because the relevel(ref = "Yes") in our code above makes _Yes the reference group and the group that a Euler plot shows is the reference group
- The area of the Alcohol_factor_Yes (dark pink) group is similar to the area of the All (white) group. This indicates that a high % of the sample were in the Alcohol_factor=Yes group.

M$Marijuana |> tab1()
M$HardDrugs |> tab1()
M$Alcohol12PlusYr |> tab1()

Let’s add a legend to make it easier to see what the colours represent. Also, to help quantify the degree of overlap between the groups, we’ll add count and % values.

plot(fit1, 
     legend = list(labels = c( "All", "Marijuana","HardDrugs", "Alcohol")),
     quantities = list(type = c("percent", "counts")))

Let’s now focus on the association between Alcohol_factor and Marijuana_factor by dropping HardDrugs_factor from our plot.

fit1 <- eulerr::euler(M[,c("Marijuana_factor",
                           "Alcohol_factor",
                           "All")] |>  na.omit(),
                          shape = "ellipse")

plot(fit1, main = "Marijuana and Alcohol")

Correlations between variables

In the previous section, we created binary variables. This allows us to look at correlations between them and other variables.

# Select variables for corelation analysis
my_data <- M[, c("Male",
                  "Marijuana_bin", "Alcohol_bin", "SmokeNow_bin",
                  "HardDrugs_bin", "Age")]

PerformanceAnalytics::chart.Correlation(my_data, 
                                        histogram=TRUE, pch=19)

Publication-quality by-group frequency tables (Table 1)

We can use the tableby() function in the arsenal library to produce a publication-quality table of count and % values for variables by Marijuana use. Usually, this is Table 1 in a publication.

tab1 <- arsenal::tableby(Marijuana ~ # column
                            SmokeCurrent + # rows
                            HardDrugs + 
                            Race1 + 
                            Age + 
                            Gender + 
                            Depressed + 
                            Alcohol12PlusYr, 
                            cat.stats="countrowpct", #row percents are best for estimating risk ratios.
                            data = M)

summary(tab1, text = TRUE)

In posit, this table can be copied to a clipboard (accessible from the File tab) and pasted into a document using the code below. (Below, we’ve named the file trash). This means you don’t have to create a table manually by typing in values – so it saves you time!
write2html(tab1, "./trash.html", title = "My table in html")
write2word(tab1, "./trash.doc", title = "My table in Word")

Logistic model

We’ll build models using logistic regression. Before we start, our categorical variables should be of the type, factor rather than character - check if they are using the str() function.

M$HardDrugs |> str()
M$SmokeNow |> str()
M$Race1 |> str()

Run this code to look at the association between HardDrugs and Marijuana.
- glm() stands for generalised linear model
- family = binomial(link = "logit") indicates that we’re using logistic regression (logit link and binomial errors).

The basic logistic model is produced with the following code:

model <- glm(Marijuana ~ HardDrugs, data = M,
             family = binomial(link = "logit"))
summary(model)

jtools::summ(model, exp = TRUE)

model <- glm(Marijuana ~ SmokeCurrent + HardDrugs + Race1 + Age + 
               Gender + Depressed + Alcohol12PlusYr, data = M,
             family = binomial(link = "logit"))

jtools::summ(model, exp = TRUE) # Exponentiate the beta coefficients to get odds ratios

Regression nomogram

This is a visual summary of the contributions of predictors to an outcome. We covered regression nomograms in the previous tutorial.

The following code draws a regression nomogram for our model from the logistic model section.

regplot(model)

The output is shown below

Nomogram

The nomogram depicts the structure of the model as a weighted average.

Each variable contributes a score β(X-m) terms based on the horizontal deviation. These are summed for each individual to give a Total score.
Reading vertically down from the Total score, then gives a probability of Marijuana use for each individual in the model (ranging between 0 and 1).

The horizontal distance between boxes or degree of spread of density plots for variables represents the relative influence of the variable on the model outcome.

The size of each box shows the distribution of categorical variables and the density plot shows the distribution of continuous variables.

Model fit

How well a model fits the data can be assessed with a marginal model plot, which shows the relationship between model predicted and observed data.

# Note, all predictor variables must be factors for marginal model plot to work....
model <- glm(Marijuana ~ factor(SmokeCurrent) + factor(HardDrugs) + factor(Race1) + Age + 
               factor(Gender) + factor(Depressed) + factor(Alcohol12PlusYr), data = M,
             family = binomial(link = "logit"))

assign("M", M, envir = .GlobalEnv)
assign("model", model, envir = .GlobalEnv)

## Hosmer-Lemeshow test of model fit

glmtoolbox::hltest(model)

## marginal model plot
car::marginalModelPlots(model)

Publication-quality regression tables (Table 2)

The code below allows us to produce publication-quality regression tables. Usually, a regression table is Table 2 in a publication. Note that the dependent variable must be a factor variable with two levels for this code to produce a logistic regression model, rather than a linear or other regression.

dependent <- "Marijuana"
explanatory <- c(
  "SmokeCurrent", 
  "HardDrugs","Race1", "Age", "Depressed", "Alcohol12PlusYr")

M |>
 finalfit(dependent, explanatory, metrics = FALSE) -> regression_table
  
  
regress_tab <- knitr::kable(regression_table, row.names = FALSE,
                            
                            align = c("l", "l", "r", "r", "r", "r"))

regress_tab |> kableExtra::kable_styling()


M |>
 finalfit(dependent, explanatory, metrics = TRUE) -> metrics_table

## show metrics only for display
metrics_table[[2]]

Summary

Today, we have extended regression to help us estimate the strength of association between risk factors and binary outcomes which are often used in epidemiology. Logistic regression is used to predict the probability of a binary outcome such as disease status.

• Logistic regression is often used when the outcome is a binary variable and there are confounders to adjust for.
• The principles are similar to usual regression, however the beta coefficients are exponentiated and interpreted as odds ratios (either crude or adjusted).
• The outcome variable must have only two levels (1 and 0, Yes or No, or TRUE or FALSE).
• We often present results in a table with both crude and adjusted odds ratios to see the effect of adjusting for confounders on an association.

Homework

Repeat the above analyses looking at risk factors for alcohol use. For the outcome, use the response to the question: consumption of at least 12 drinks of any type of alcoholic beverage in any one year (Alcohol12PlusYr == "Yes").

Please complete this in posit.