Part I

In the first part of this class activity, we will explore quantile residuals as a diagnostic tool for logistic regression. Quantile residuals can be calculated for a fitted GLM in R using the qresid() function in the statmod package (you may need to install this package).

In class, we used the following code to simulate data for which the logistic regression model assumptions hold:

# simulate a single explanatory variable from a Normal distribution
x <- rnorm(1000)

# create P(Y = 1 | X) for each entry in x
# Here log odds = -1 + 2x
p <- exp(-1 + 2*x)/(1 + exp(-1 + 2*x))

# Finally, simulate a binary response at each x
y <- rbinom(1000, 1, p)

# fit the model and plot the quantile residuals against x
# add a smooth fit to see if there is a relationship
m1 <- glm(y ~ x, family = binomial)
data.frame(x = x, residuals = qresid(m1)) %>%
  ggplot(aes(x = x, y = residuals)) +
  geom_point() +
  geom_smooth() +
  theme_bw()

Now, let’s see what happens when we break a logistic regression assumption!

Part II

In the second part of this class activity, we will use diagnostic tools to help fit a logistic regression model on real data. Here we work with data from a website called ScienceForums.Net (SFN), which has been open since 2002 and hosts conversations on a range of topics from biological and physical science to religion and philosophy. Each row in the data represents one ‘thread’, which is comprised of a series of posts stemming from an initial post. For each thread, we have some information that SFN collects such as the number of views and the number of authors. The threads present in the data are a random sample of threads from 2002-2014, with the data collected in 2014. SFN moderators are interested in using this data to determine which threads warrant the most attention.

You can load the SFN data into R by

sfn <- read.csv("https://sta712-f22.github.io/class_activities/sfn.csv")

The sfn dataset contains the following columns:

• Age: the age of the thread (in days) when the data was collected in 2014, measured from the first post in the thread
• State: sometimes moderators close threads if they are inappropriate. closed indicates the thread has been closed, otherwise State is open
• Posts: the number of posts in the thread
• Views: the total number of views of the thread
• Duration: the number of days between the first and last posts in the thread
• Authors: the number of distinct authors posting in the thread
• AuthorExperience: the number of days the author of the first post in the thread had been registered on SFN when the thread began (0 indicates they registered that day)
• DeletedPosts: the number of posts in the thread that have been deleted by a moderator
• Forum: the forum in which the thread was posted (e.g., Science)
• AuthorBanned: whether the original author of the thread is currently banned from posting on SFN (at the time of data collection, not when the thread was first posted)

logodds_plot <- function(data, num_bins, bin_method,
                         x_name, y_name, grouping = NULL, reg_formula = y ~ x){
  
  if(is.null(grouping)){
    dat <- data.frame(x = data %>% pull(x_name), 
                      y = data %>% pull(y_name),
                      group = 1)
  } else {
    dat <- data.frame(x = data %>% pull(x_name), 
                      y = data %>% pull(y_name),
                      group = data %>% pull(grouping))
  }
  
  if(bin_method == "equal_size"){
    logodds_table <- dat %>%
      drop_na() %>%
      arrange(group, x) %>%
      group_by(group) %>%
      mutate(obs = y,
             bin = rep(1:num_bins,
                       each=ceiling(n()/num_bins))[1:n()]) %>%
      group_by(bin, group) %>%
      summarize(mean_x = mean(x),
                prop = mean(c(obs, 0.5)),
                num_obs = n()) %>%
      ungroup() %>%
      mutate(logodds = log(prop/(1 - prop)))
  } else {
    logodds_table <- dat %>%
      drop_na() %>%
      group_by(group) %>%
      mutate(obs = y,
             bin = cut(x, 
                       breaks = num_bins,
                       labels = F)) %>%
      group_by(bin, group) %>%
      summarize(mean_x = mean(x),
                prop = mean(c(obs, 0.5)),
                num_obs = n()) %>%
      ungroup() %>%
      mutate(logodds = log(prop/(1 - prop)))
  }
  
  if(is.null(grouping)){
    logodds_table %>%
      ggplot(aes(x = mean_x,
                 y = logodds)) +
      geom_point(size=2) +
      geom_smooth(se=F, method="lm", formula = reg_formula) +
      theme_bw() +
      labs(x = x_name,
           y = "Empirical log odds") +
      theme(text = element_text(size=25))
  } else {
    logodds_table %>%
      ggplot(aes(x = mean_x,
                 y = logodds,
                 color = group,
                 shape = group)) +
      geom_point(size=2) +
      geom_smooth(se=F, method="lm", formula = reg_formula) +
      theme_bw() +
      labs(x = x_name,
           y = "Empirical log odds",
           color = grouping,
           shape = grouping) +
      theme(text = element_text(size=25))
  }
  
}