Bernoulli + Binomial random variables (1/2)

Logistic regression is used to analyze data with two types of responses:

• Binary: These responses take on two values success \((Y = 1)\) or failure \((Y = 0)\), yes \((Y = 1)\) or no \((Y = 0)\), etc.

Bernoulli + Binomial random variables (2/2)

Logistic regression is used to analyze data with two types of responses:

• Binomial: Number of successes in a Bernoulli process, \(n\) independent trials with a constant probability of success \(p\).

Binary vs. Binomial data

For each example, identify if the response is a Bernoulli or Binomial response

1. Use median age and unemployment rate in a county to predict the percent of Obama votes in the county in the 2008 presidential election.

2. Use GPA and MCAT scores to estimate the probability a student is accepted into medical school.

3. Use sex, age, and smoking history to estimate the probability an individual has lung cancer.

4. Use offensive and defensive statistics from the 2017-2018 NBA season to predict a team’s winning percentage.

Logistic regression model

\(\log\Big(\frac{p}{1-p}\Big) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \beta_px_p\)

• The response variable, \(\log\Big(\frac{p}{1-p}\Big)\), is the log(odds) of success, i.e. the logit

• Use the model to calculate the probability of success \(\hat{p} = \frac{e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \beta_px_p}}{1 + e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \beta_px_p}}\)

• When the response is a Bernoulli random variable, the probabilities can be used to classify each observation as a success or failure

Logistic vs linear regression model

set.seed(0)
dat <- tibble(x=runif(200, -5, 10),
                  p=exp(-2+1*x)/(1+exp(-2+1*x)),
                  y=rbinom(200, 1, p),
                  y2=.3408+.0901*x,
                  logit=log(p/(1-p)))
dat2 <- tibble(x = c(dat$x, dat$x),
               y = c(dat$y2, dat$p),
               `Regression model` = c(rep("linear", 200),
                                      rep("logistic", 200)))

ggplot() + 
  geom_point(data = dat, aes(x, y)) +
  geom_line(data = dat2, aes(x, y, linetype = `Regression model`, 
                             color = `Regression model`)) + 
  labs(title = "Linear vs. logistic regression models for binary response data") + 
  scale_colour_manual(name = 'Regression model',
                      values = c('blue', 'red'), 
                      labels = c('linear', 'logistic'), guide ='legend')

Logit link

Bernoulli and Binomial random variables can be written in one-parameter exponential family form, \(f(y;\theta) = e^{[a(y)b(\theta) + c(\theta) + d(y)]}\)

Assumptions for logistic regression

The following assumptions need to be satisfied to use logistic regression to make inferences

COVID-19 infection prevention practices at food establishments

Researchers at Wollo Univeristy in Ethiopia conducted a study in July and August 2020 to understand factors associated with good COVID-19 infection prevention practices at food establishments. Their study is published in Andualem et al. (2022)

They were particularly interested in the understanding implementation of prevention practices at food establishments, given the workers’ increased risk due to daily contact with customers.

The data

“An institution-based cross-sectional study was conducted among 422 food handlers in Dessie City and Kombolcha Town food and drink establishments in July and August 2020. The study participants were selected using a simple random sampling technique. Data were collected by trained data collectors using a pretested structured questionnaire and an on-the-spot observational checklist.”

Response variable

“The outcome variable of this study was the good or poor practices of COVID-19 infection prevention among food handlers. Nine yes/no questions, one observational checklist and five multiple choice infection prevention practices questions were asked with a minimum score of 1 and maximum score of 25. Good infection prevention practice (the variable of interest) was determined for food handlers who scored 75% or above, whereas poor infection prevention practices refers to those food handlers who scored below 75% on the practice questions.”

Interpreting the results

• Is the response a Bernoulli or Binomial?

• What is the strongest predictor of having good COVID-19 infection prevention practices?

  • It’s often unreliable to look answer this question just based on the model output. Why are we able to answer this question based on the model output in this case?

• Describe the effect (coefficient interpretation and inference) of having COVID-19 infection prevention policies available at the food establishment.

• The intercept describes what group of food handlers?

Access to personal protective equipment

We will use the data from Andualem et al. (2022) to explore the association between age, sex, years of service, and whether someone works at a food establishment with access to personal protective equipment (PPE) as of August 2020. We will use access to PPE as a proxy for wearing PPE.

covid_df <- read_csv("data/covid-prevention-study.csv") |>
  rename(age = "Age of food handlers", 
         years = "Years of service", 
         ppe_access = "Availability of PPEs") |>
  mutate(sex = factor(if_else(Sex == 2, "Female", "Male"))) |>
  select(age, sex, years, ppe_access) 

covid_df |> slice(1:5) |> kable()

EDA for binary response (2/2)

library(viridis)
ggplot(data = covid_df, aes(x = sex, fill = factor(ppe_access))) + 
  geom_bar(position = "fill")  +
  labs(x = "Sex", 
       fill = "PPE Access", 
       title = "PPE Access by Sex") + 
  scale_fill_viridis_d()

Visualizing coefficient estimates

model_coef <- tidy(ppe_model, exponentiate = TRUE, conf.int = TRUE)

ggplot(data = model_coef, aes(x = term, y = estimate)) +
  geom_point() +
  geom_hline(yintercept = 1, lty = 2) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high))+
  labs(title = "Exponentiated model coefficients") + 
  coord_flip()

Data: Supporting railroads in the 1870s

The data set RR_Data_Hale.csv contains information on support for referendums related to railroad subsidies for 11 communities in Alabama in the 1870s. The data were originally analyzed as part of a thesis project by a student at St. Olaf College. The variables in the data are

• pctBlack: percentage of Black residents in the county
• distance: distance the proposed railroad is from the community (in miles)
• YesVotes: number of “yes” votes in favor of the proposed railroad line
• NumVotes: number of votes cast in the election

EDA (1/3)

rr <- rr |>
  mutate(pctYes = YesVotes/NumVotes, 
         emp_logit = log(pctYes / (1 - pctYes)))

rr |> head(5)

# A tibble: 5 × 10
  County   popBlack popWhite popTotal pctBlack distance YesVotes NumVotes pctYes emp_logit
  <chr>       <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>  <dbl>     <dbl>
1 Carthage      841      599     1440     58.4       17       61      110 0.555      0.219
2 Cedervi…     1774      146     1920     92.4        7        0       15 0       -Inf    
3 Five Mi…      140      626      766     18.3       15        4       42 0.0952    -2.25 
4 Greensb…     1425      975     2400     59.4        0     1790     1804 0.992      4.85 
5 Harrison      443      355      798     55.5        7        0       15 0       -Inf    

EDA (2/3)

p1 <- ggplot(data = rr, aes(x = distance, y = emp_logit)) + 
  geom_point() + 
  geom_smooth(method  = "lm", se = FALSE) + 
  labs(x = "Distance to proposed railroad", 
       y = " ")
  
p2 <- ggplot(data = rr, aes(x = pctBlack, y = emp_logit)) + 
  geom_point() + 
  geom_smooth(method  = "lm", se = FALSE) + 
  labs(x = "% Black residents", 
       y = "")
p1 + p2 + plot_annotation(title = "Log(odds yes vote) vs. predictor variables")

EDA (3/3)

rr <- rr |>
  mutate(inFavor = if_else(pctYes > 0.5, "Yes", "No"))

ggplot(data = rr, aes(x = distance, y = pctBlack, color = inFavor)) + 
  geom_point() + 
  geom_smooth(method  = "lm", se = FALSE, aes(lty = inFavor)) + 
  labs(x = "Distance to proposed railroad", 
       y = "% Black residents",
       title = "% Black residents vs. distance", 
       subtitle = "Based on vote outcome") + 
  scale_color_viridis_d(end = 0.85)

Model

Let \(p\) be the percent of yes votes in a county. We’ll start by fitting the following model:

\[\log\Big(\frac{p}{1-p}\Big) = \beta_0 + \beta_1 ~ dist + \beta_2 ~ pctBlack\]

Model in R

rr_model <- glm(cbind(YesVotes, NumVotes - YesVotes) ~ distance + pctBlack, 
                data = rr, family = binomial)
tidy(rr_model, conf.int = TRUE) |>
  kable(digits = 3)

Residuals

Similar to Poisson regression, there are two types of residuals: Pearson and deviance residuals

Plot of deviance residuals

rr_int_model <- glm(cbind(YesVotes, NumVotes - YesVotes) ~ distance + pctBlack +
                      distance*pctBlack, 
                data = rr, family = binomial)

rr_int_aug <- augment(rr_int_model, type.predict = "response", 
                        type.residuals = "deviance")

rr_int_aug |> slice(1:5) |> kable()

ggplot(data = rr_int_aug, aes(x = .fitted, y = .resid)) + 
  geom_point() + 
  geom_hline(yintercept = 0, color = "red") + 
  labs(x = "Fitted values", 
       y = "Deviance residuals", 
       title = "Deviance residuals vs. fitted", 
       subtitle = "for model with interaction term")

Goodness of fit

Similar to Poisson regression, the sum of the squared deviance residuals is used to assess goodness of fit.

\[\begin{aligned} &H_0: \text{ Model is a good fit} \\ &H_a: \text{ Model is not a good fit}\end{aligned}\]

• When \(m_i\) is large and the model is a good fit \((H_0 \text{ true})\) the residual deviance follows a \(\chi^2\) distribution with \(n - p\) degrees of freedom.
  • Recall \(n - p\) is the residual degrees of freedom.
• If the model fits, we expect the residual deviance to be approximately what value?

Adjusting for overdispersion (1/2)

• Overdispersion occurs when there is extra-binomial variation, i.e. the variance is greater than what we would expect, \(np(1-p)\).

• Similar to Poisson regression, we can adjust for overdispersion in the binomial regression model by using a dispersion parameter \[\hat{\phi} = \sum \frac{(\text{Pearson residuals})^2}{n-p}\]

  • By multiplying by \(\hat{\phi}\), we are accounting for the reduction in information we would expect from independent observations.

Adjusting for overdispersion (2/2)

• We adjust for overdispersion using a quasibinomial model.
  • “Quasi” reflects the fact we are no longer using a binomial model with true likelihood.
• The standard errors of the coefficients are \(SE_{Q}(\hat{\beta}_j) = \sqrt{\hat{\phi}} SE(\hat{\beta})\)
  • Inference is done using the \(t\) distribution to account for extra variability

Predicting ED wait and treatment times

Ataman and Sariyer (2021) use ordinal logistic regression to predict patient wait and treatment times in an emergency department (ED). The goal is to identify relevant factors that can be used to inform recommendations for reducing wait and treatment times, thus improving the quality of care in the ED.

Data: Daily records for ED arrivals in August 2018 at a public hospital in Izmir, Turkey.

Response variable: Wait time, a categorical variable with three levels: - Patients who wait less than 10 minutes - Patients whose waiting time is in the range of 10-60 minutes - Patients who wait more than 60 minutes

Ordinal logistic regression

Let \(Y\) be an ordinal response variable that takes levels \(1, 2, \ldots, J\) with associated probabilities \(p_1, p_2, \ldots, p_J\).

The proportional odds model can be written as the following:

\[\begin{aligned}&\log\Big(\frac{P(Y\leq 1)}{P(Y > 1)}\Big) = \beta_{01} + \beta_1x_1 + \dots + \beta_px_p \\ & \log\Big(\frac{P(Y\leq 2)}{P(Y > 2)}\Big) = \beta_{02} + \beta_1x_1 + \dots + \beta_px_p \\ & \dots \\ & \log\Big(\frac{P(Y\leq J-1)}{P(Y > J-1)}\Big) = \beta_{0{J-1}} + \beta_1x_1 + \dots + \beta_px_p\end{aligned}\]

Effect of arrival mode

Effect of triage level

Wrap up

• Covered fitting, interpreting, and drawing conclusions from GLMs

  • Looked at Poisson, Negative Binomial, and Logistic (binary, binomial, ordinal) in detail

• Used Pearson and deviance residuals to assess model fit and determine if new variables should be added to the model

• Addressed issues of overdispersion and zero-inflation

• Used the properties of the one-parameter exponential family to identify the best link function for any GLM