Chapter 11 Multiple Logistic Regression

11.1 Packages Needed for Multiple Logistic Regression

This code will check that required packages for this chapter are installed, install them if needed, and load them into your session.

req <- substitute(require(x, character.only = TRUE))
libs<-c("sjPlot", "ggplot2", "jtools", "car", "blorr", "DescTools", "MASS")
sapply(libs, function(x) eval(req) || {install.packages(x); eval(req)})

11.2 Data Prep for Multiple Logistic Regression

Declare factor variables as such.

class(mydata$var) # will show you how R sees the specified variable (double, factor, etc.)
mydata$catvar <- factor(mydata$catvar) # will reclassify specified variable as a factor variable

Change the reference level on a factor variable (specify the desired level; not determined by assigned number)

mydata$catvar <- relevel(mydata$catvar, ref = #)

11.3 Basic Multiple Logistic Regression Commands

object <- glm(dv ~ iv1 + iv2 + iv3, data = mydata, family = "binomial")
summary(object) # results in logit coefficients
exp(cbind(OR = coef(object), confint(object))) # results in ORs and their CIs

11.4 Model Fit Statistics for Multiple Logistic Regression

blorr::blr_model_fit_stats(object) # Gives various fit statistics
blorr::blr_test_hosmer_lemeshow(object) # Hosmer Lemeshow gof test
blorr::blr_roc_curve(blr_gains_table(object)) # ROC curve
DescTools::Cstat(object) # C-Statistic (concordance statistic)

11.5 Diagnostics for Multiple Logistic Regression

Logistic regression assumes: 1) The outcome is dichotomous; 2) There is a linear relationship between the logit of the outcome and each continuous predictor variable; 3) There are no influential cases/outliers; 4) There is no multicollinearity among the predictors. For now, I just have two commands that will provide VIFs (multicollinearity detection).

car::vif(object)
DescTools::VIF(object) # an alternative

11.6 Producing Formatted Tables of Multiple Logistic Regression Results

sjPlot’s tab_model function works really well, especially if you only have one to three models.

sjPlot::tab_model(object)
sjPlot::tab_model(object1, object2, pred.labels = c("Intercept", "iv1 label", "iv2 label"),
dv.labels = c("Model 1", "Model 2"), string.pred = "Estimates", string.ci = "95% CIs", string.p = "p-values")

11.7 Graphing Coefficients and CIs for Multiple Logistic Regression

(using sjPlot’s plot_model function, dotwhisker’s dwplot function, or coefplot’s function.

sjPlot::plot_model(object) # Default will display odds ratios (the other two won’t)
dotwhisker::dwplot(object) # No intercept is displayed; in logit coefficients
coefplot::coefplot(object) # Intercept is included; in logit coefficients

11.8 Graphing Margins (Predicted Probabilities) for Multiple Logistic Regression

jtools::effect_plot(object, pred = iv1, interval = TRUE, y.label = "label", x.label = "label")

11.9 Interactions (modeling and graphing) for Multiple Logistic Regression

See sjPlot or interactions pages for more information and argument options.

When interacting a continuous variable with a categorical variable:

object <- glm(dv ~ intvar * catvar + var, data = mydata, family = "binomial")
summary(object)
interactions::interact_plot(object, pred = intvar, modx = catvar)

# Or

object <- glm(dv ~ intvar * catvar + var, data = mydata, family = "binomial")
sjPlot::plot_model(object, type = "pred", terms = c("intvar", "catvar"))

# When interacting two categorical variables

object <- glm(dv ~ catvar1 * catvar2 + var, data = mydata, family = "binomial")
sjPlot::plot_model(object, type = "pred", terms = c("catvar1", "catvar2")) # can switch the order

# Or

interactions::cat_plot(object, pred = catvar1, modx = catvar2) # catvar1 will be the x-axis var

11.10 Ordered Logistic Regression

object <- MASS::polr(dv ~ iv1 + iv2 + iv3, data = mydata, Hess = TRUE)

11.11 Consolidated Code for Multiple Logistic Regression

Below is the consolidated code from this chapter. One could transfer this code into an empty RScript, which also offers the option of find/replace terms. You can also download this generic Multiple Logistic Regression RScript file here

Placeholders that need replacing:

mydata – name of your dataset
var – general variable(s)
depvar, indvar1, indvar2, etc – general variables
catvar – name of your categorical variable
intvar – name of your interval or continuous variable
object(s) – whatever you want to call your object(s))
number – Specified number
labels/title – any titles, axis labels, category labels

# Multiple Logistic Regression -- Generic RScript


# Packages Needed

  req <- substitute(require(x, character.only = TRUE))
  libs<-c("sjPlot", "ggplot2", "jtools", "car", "blorr", "DescTools", "MASS")
  sapply(libs, function(x) eval(req) || {install.packages(x); eval(req)})


# Prep Work

## Declare factor variables as such.
  class(mydata$var)
  mydata$catvar <- factor(mydata$catvar)

## Change the reference level on a factor variable (specify the desired level; not determined by assigned number)
  mydata$catvar <- relevel(mydata$catvar, ref = number)


# Basic Logistic Regression Commands

    object <- glm(depvar ~ indvar1 + indvar2 + indvar3, data = mydata, family = "binomial")
   summary(object)          # results in logit coefficients
    exp(cbind(OR = coef(object), confint(object)))  # results in ORs and their CIs


# Model Fit Statistics

  blorr::blr_model_fit_stats(object)        # Gives various fit statistics
  blorr::blr_test_hosmer_lemeshow(object)   # Hosmer Lemeshow gof test
  blorr::blr_roc_curve(blr_gains_table(object)) # ROC curve
  DescTools::Cstat(object)              # C-Statistic (concordance statistic)


# Diagnostics

  car::vif(object)
  DescTools::VIF(object)    # an alternative


# Producing formatted tables of regression results

  sjPlot::tab_model(object)
  sjPlot::tab_model(object1, object2, pred.labels = c("Intercept", "indvar1 label", "indvar2 label"),
    dv.labels = c("Model 1", "Model 2"), string.pred = "Estimates", string.ci = "95% CIs", string.p = "p-values")


# Graphing Coefficients and CIs (options)

  sjPlot::plot_model(object)    # Default will display odds ratios (the other two won’t)
  dotwhisker::dwplot(object)    # No intercept is displayed; in logit coefficients
  coefplot::coefplot(object)    # Intercept is included; in logit coefficients


# Graphing margins/predicted values for specified "X"

  jtools::effect_plot(object, pred = indvar1, interval = TRUE, y.label = "label", x.label = "label")


# Interactions (modeling and graphing)

## When interacting a continuous variable with a categorical variable

  object <- glm(depvar ~ intvar * catvar + var, data = mydata, family = "binomial")
  summary(object)
  interactions::interact_plot(object, pred = intvar, modx = catvar)

## Or

  object <- glm(depvar ~ intvar * catvar + var, data = mydata, family = "binomial")
  sjPlot::plot_model(object, type = "pred", terms = c("intvar", "catvar"))


## When interacting two categorical variables

  object <- glm(depvar ~ catvar1 * catvar2 + var, data = mydata, family = "binomial")
  sjPlot::plot_model(object, type = "pred", terms = c("catvar1", "catvar2"))  # can switch the order

## Or

  interactions::cat_plot(object, pred = catvar1, modx = catvar2)  # catvar1 will be the x-axis var


# Ordered Logistic Regression

  object <- MASS::polr(depvar ~ indvar1 + indvar2 + indvar3, data = mydata, Hess = TRUE)