04-Logistic.Rmd

# 4 Logistic Regression {#Logistic}

## 4.1 The Logistic Regression Model

### 4.1.1 The Logistic Regression Model

\begin{equation} 
  \mathrm{logit}[\pi(x)] = \mathrm{log}\left[\frac{\pi(x)}{1=\pi(x)}\right]= \alpha + \beta x.
  (\#eq:logitOfSuccess)
\end{equation}

\begin{equation} 
  \pi(x) = \frac{e^{\alpha + \beta x}}{1+e^{\alpha + \beta x}}.
  (\#eq:exponentialFunction)
\end{equation} 

### 4.1.2 Odds Ratio and Linear Approximaiton Interpreetations {#4.1.2}

$$\frac{\pi(x)}{1-\pi(x)}=\mathrm{exp}(\alpha + \beta x) = e^\alpha(e^\beta)^x.$$  

```{r}
Crabs <- read.table("http://users.stat.ufl.edu/~aa/cat/data/Crabs.dat",
  header = TRUE
)
library(RColorBrewer)
colors <- brewer.pal(n = 4, name = "Dark2")

fit <- glm(y ~ width, family = binomial, data = Crabs)
summary(fit)

at28 <- predict(fit, data.frame(width = 28), type = "resp")
at29 <- predict(fit, data.frame(width = 29), type = "resp")

x <- seq(20, 34, by = .01)

df <- tibble(x = x) %>%
  mutate(y = predict(fit, data.frame(width = x), type = "resp"))


spl <- smooth.spline(df$x, df$y, spar = 0.3)
newx <- seq(min(df$x), max(df$x), 0.1)
pred <- predict(spl, x = newx, deriv = 0)

# solve for tangent at a given x
newx <- 28
pred0 <- predict(spl, x = newx, deriv = 0)
pred1 <- predict(spl, x = newx, deriv = 1)
yint <- pred0$y - (pred1$y * newx)
xint <- -yint / pred1$y

tang <- tibble(x = df$x) %>%
  mutate(y = yint + pred1$y * x) %>%
  filter(x > 25.5 & x < 30.5)

library(ggthemes)

aLabel <- as.vector(expression(paste(beta, pi, "(1-", pi, ")")))

ggplot() +
  geom_line(data = df, aes(x = x, y = y)) +
  geom_line(data = tang, aes(x = x, y = y)) +
  annotate(
    geom = "segment", x = 28, y = at28, xend = 28, yend = at29,
    color = colors[1]
  ) +
  annotate(
    geom = "segment", x = 28, y = at29, xend = 29, yend = at29,
    color = colors[1]
  ) +
  annotate(
    geom = "text", x = 26.35, y = at28 + .02,
    label = 'paste(beta, pi, "(1-" , pi, ")")', parse = TRUE,
    hjust = "left",
    color = colors[1]
  ) +
  annotate(
    geom = "text", x = 28.3, y = at29 + .03,
    label = "1",
    hjust = "left",
    color = colors[1]
  ) +
  theme_few() +
  theme(axis.title.y = element_text(angle = 0, vjust = .5)) +
  scale_x_continuous(breaks = seq(20, 34, by = 2)) +
  scale_y_continuous(breaks = seq(0, 1, by = .2)) +
  ggtitle("Figure 4.1") +
  xlab("Width") +
  ylab("Probability")
```      

### 4.1.3 Example: Whethr a Female Horsehoe Crab Has Satelites {#x4.1.3}

```{r}
Crabs <- read.table("http://users.stat.ufl.edu/~aa/cat/data/Crabs.dat",
  header = TRUE
)

Crabs %>%
  filter(row_number() %in% c(1, 2, n()))

library(RColorBrewer)
colors <- brewer.pal(n = 4, name = "Dark2")

library(gam)
gam.fit <- gam(y ~ s(width), family = binomial, data = Crabs)
fit <- glm(y ~ width, family = binomial, data = Crabs)

plot(jitter(y, 0.08) ~ width, data = Crabs, xlab = "Width", ylab = "Satellites")
curve(predict(gam.fit, data.frame(width = x), type = "resp"),
  col = colors[1], add = TRUE
)
curve(predict(fit, data.frame(width = x), type = "resp"),
  col = colors[2], add = TRUE
)

summary(fit)
# estimated probability of satellite at width = 21.0
predict(fit, data.frame(width = 21.0), type = "resp")

predict(fit, data.frame(width = mean(Crabs$width)), type = "resp")
```

$$\mathrm{logit}[\hat\pi(x)] = -12.351 + 0.497 x.$$  

$$\hat\pi(x) =\frac{\mathrm{exp}(-12.351 + 0.497 x)}{1+\mathrm{exp}(-12.351 + 0.497 x)}.$$ 

$$\frac{\mathrm{exp}[-12.351 + 0.497 (21.0)]}{1+\mathrm{exp}[-12.351 + 0.497 (21.0)]} = 0.129.$$

### 4.1.4 Logistic Regression with Retrospective Studies

### 4.1.5 Normally Distributed X Implies Logistic Regression for Y

## 4.2 Statistical Inference for Logistic Regression {#x4.2}

### 4.2.1 Confidence Intervals for Effects {#x4.2.1}

```{r}
Crabs <- read.table("http://users.stat.ufl.edu/~aa/cat/data/Crabs.dat",
  header = TRUE
)

fit <- glm(y ~ width, family = binomial, data = Crabs)
summary(fit) # z value & p-value Wald test

suppressMessages(
  confint(fit) # profile likelihood confidence interval
)

suppressPackageStartupMessages(library(car))
Anova(fit) # likelihood-ratio test of width effect
```

### 4.2.2 Significance Testing {#x4.2.2}

### 4.2.3 Fitted Values and Confidence Intervals for Probabilities
$$\hat P(Y=1)=\mathrm{exp}(\hat\alpha + \hat\beta x)/[1 + (\hat\alpha + \hat\beta x)],$$ 

```{r}
Crabs <- read.table("http://users.stat.ufl.edu/~aa/cat/data/Crabs.dat",
  header = TRUE, stringsAsFactors = TRUE
)

fit <- glm(y ~ width, family = binomial, data = Crabs)
pred.prob <- fitted(fit) # ML fitted value estimate of P(Y = 1)
lp <- predict(fit, se.fit = TRUE) # linear predictor
LB <- lp$fit - 1.96 * lp$se.fit # confidence bounds for linear predictor
LB <- lp$fit + qnorm(0.025) * lp$se.fit # better confidence bound
UB <- lp$fit + qnorm(0.975) * lp$se.fit
LB.p <- exp(LB) / (1 + exp(LB)) # confidence bounds for P(Y = 1)
UB.p <- exp(UB) / (1 + exp(UB))


library(dplyr) # bind_cols, filter, row_number
predictions <-
  bind_cols(
    Width = Crabs$width,
    `Predicted probaility` = pred.prob,
    `Lower CB` = LB.p,
    `Upper CB` = UB.p
  )

predictions %>%
  filter(row_number() %in% c(1, 7, n()))
```

```{r}
fit <- glm(y ~ width, family = binomial, data = Crabs)

library(tidyverse)
data.plot <- tibble(width = 18:34)
lp <- predict(fit, newdata = data.plot, se.fit = TRUE)

library(ggthemes)
data.plot <- data.plot %>%
  mutate(
    pred.prob = exp(lp$fit) / (1 + exp(lp$fit)),
    LB = lp$fit + qnorm(0.025) * lp$se.fit,
    UB = lp$fit + qnorm(0.975) * lp$se.fit,
    LB.p = exp(LB) / (1 + exp(LB)),
    UB.p = exp(UB) / (1 + exp(UB))
  )

data.plot %>%
  ggplot(aes(x = width)) +
  geom_point(data = Crabs, aes(x = width, y = jitter(y, .1))) +
  geom_line(aes(y = pred.prob)) +
  geom_ribbon(aes(ymin = LB.p, ymax = UB.p), fill = "gray", alpha = 0.5) +
  theme_few() +
  ylab("Probability (satellites)") +
  scale_y_continuous(breaks = seq(0, 1, by = .2)) +
  xlab("Width") +
  ggtitle("Figure 4.3")
```

### 4.2.4 Why Use a Model to Estimate Probabilities?

## 4.3 Logistic Regression with Categorical Predictors

### 4.3.1 Indicator Variables Represent Categories of Predictors
$$\mathrm{logit}[P(Y = 1)] = \alpha + \beta_1 x + \beta_2 z$$

| x | z | Logit |
| :---: | :---: | :---: |
| 0 | 0 | $\alpha$ |
| 1 | 0 | $\alpha + \beta_1$ |
| 0 | 0 | $\alpha + \beta_2$ |
| 0 | 0 | $\alpha + \beta_1 + \beta_2$ |

$$ = [\alpha + \beta_1 (1) + \beta_2 z] - [\alpha + \beta_1 (0) + \beta_2 z] = \beta_a.$$

### 4.3.2 Example: Survey about Marijuana Use

```{r}
library(tidyverse)
Marijuana <- read.table("http://users.stat.ufl.edu/~aa/cat/data/Marijuana.dat",
  header = TRUE, stringsAsFactors = TRUE
)

Marijuana

fit <- glm(yes / (yes + no) ~ gender + race,
  weights = yes + no, family = binomial,
  data = Marijuana
)

theFit <- summary(fit)
theFit

# pull the 2nd and 3rd element from the "Estimate" column, exponentiate & round
round(exp(theFit$coefficients[2:3, "Estimate"]), 2)

# easier to follow syntax
theCoef <- as.data.frame(coef(theFit))
theCoef %>%
  rownames_to_column(var = "Effect") %>%
  select("Effect", "Estimate") %>%
  filter(Effect != "(Intercept)") %>%
  mutate(`Odds Ratio` = exp(Estimate)) %>%
  mutate_if(is.numeric, round, digits = 2)
```

```{r}
# library(car)
car::Anova(fit) # likelihood-ratio test for individual explanatory variables
```

### 4.3.3 ANOVA-Type Model Representation of Factors
$$\mathrm{logit}[P(Y=1)] = \alpha + \beta_i^x + \beta_k^z$$

### 4.3.4 Tests of Conditoinal Independence and of Homogeneity for Three-Way Contingency Tables

$$\mathrm{logit}[P(Y = 1)] = \alpha + \beta x + \beta_k^z,$$

## 4.4 Multiple Logistic Regression

$$\mathrm{logit}[P(Y=1)] = \alpha + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p.$$

### 4.4.1 Example: Horseshoe Crabs with Color and Width Predictors {#x4.4.1}

\begin{equation} 
  \mathrm{logit}[P(Y=1)] = \alpha + \beta_1 x + \beta_2 c_2 + \beta_3 c_3 + \beta_4 c_4.
  (\#eq:eq4dot3)
\end{equation} 

```{r}
Crabs <- read.table("http://users.stat.ufl.edu/~aa/cat/data/Crabs.dat",
  header = TRUE, stringsAsFactors = TRUE
)

fit <- glm(y ~ width + factor(color), family = binomial, data = Crabs)
summary(fit)
```

$$\mathrm{exp}[12.715 + 0.468(26.3)]/\{1 + \mathrm{exp}[12.715 + 0.468(26.3)]\}=0.399.$$  
$$\mathrm{exp}[11.385 + 0.468(26.3)]/\{1 + \mathrm{exp}[11.385 + 0.468(26.3)]\} = 0.715.$$

### 4.4.2 Model Comparison to Check Whether a Term is Needed

```{r}
summary(glm(y ~ width, family = binomial, data = Crabs))

# compare deviance = 194.45 vs 187.46 with color also in model
library(car)
Anova(glm(y ~ width + factor(color), family = binomial, data = Crabs))
```

### 4.4.3 Example: Treating Color as Quantitative or Binary

\begin{equation} 
  \mathrm{logit}[P(Y=1)] = \alpha + \beta_1 x + \beta_2 c.
  (\#eq:eq4dot4)
\end{equation} 


```{r}
fit2 <- (glm(y ~ width + color, family = binomial, data = Crabs))
summary(fit2)

anova(fit2, fit, test = "LRT") # likelihood-ratio test comparing models
```

```{r}
Crabs$c4 <- ifelse(Crabs$color == 4, 1, 0)
# Crabs$c4 <- I(Crabs$color == 4)

fit3 <- glm(y ~ width + c4, family = binomial, data = Crabs)
summary(fit3)

anova(fit3, fit, test = "LRT")
```

### 4.4.4 Allowing Interaction between Explanatory Variables

```{r}
(interaction <- glm(y ~ width + c4 + width:c4, family = binomial, data = Crabs))
summary(interaction)
```

```{r}
# P-value 0.28 from last paragraph of 4.4.4
car::Anova(interaction)
```

### 4.4.5 Effects Depend on Other Explanatory Variables in Model

## 4.5 Summarizing Effects in Logistic Regression

### 4.5.1 Probability-Based Interpretations

```{r}

fit3 <- glm(y ~ width + c4, family = binomial, data = Crabs)

round(
  predict(fit3,
    data.frame(c4 = 1, width = mean(Crabs$width)),
    type = "response"
  ),
  3
)

round(
  predict(fit3,
    data.frame(c4 = 0, width = mean(Crabs$width)),
    type = "response"
  ),
  3
)

round(
  predict(fit3, data.frame(c4 = mean(Crabs$c4), width = quantile(Crabs$width)),
    type = "response"
  ),
  3
)
```

### 4.5.2 Marginal Effects and Their Average {#x4.5.2}

```{r}
fit3 <- glm(y ~ width + c4, family = binomial, data = Crabs)

suppressPackageStartupMessages(library(mfx))

# with atmean = TRUE, finds effect only at the mean
logitmfx(fit3, atmean = FALSE, data = Crabs)
```

### 4.5.3 Standardized Interpretations {#x4.5.3}

## 4.6 Summzarizing Predictive Power: Classification Tables, ROC Curves and Multiple Correlation

### 4.6.1 Summarizing Predictive Power: Classification Tables {#x4.6.1}

```{r}
Crabs <- read.table("http://users.stat.ufl.edu/~aa/cat/data/Crabs.dat",
  header = TRUE, stringsAsFactors = TRUE
)

prop <- sum(Crabs$y) / nrow(Crabs)
prop

fit <- glm(y ~ width + factor(color), family = binomial, data = Crabs)

predicted <- as.numeric(fitted(fit) > prop) # predict y=1 when est. > 0.6416

xtabs(~ Crabs$y + predicted)
```

```{r, fig.cap="Table 4.4"}
# Make Table 4.4

predicted50 <- as.numeric(fitted(fit) > 0.50) # predict y=1 when est. > 0.50
xtabs(~ Crabs$y + predicted50)

# reorder levels to flip rows and columns
t1 <- xtabs(~ relevel(as.factor(Crabs$y), "1") + relevel(as.factor(predicted), "1"))
t1df <- as.data.frame.matrix(t1) %>% 
  rename("t1_1" = `1`) %>% 
  rename("t1_0" = `0`) 
  
t2 <- xtabs(~ relevel(as.factor(Crabs$y), "1") + relevel(as.factor(predicted50), "1"))
t2df <- as.data.frame.matrix(t2) %>% 
  rename("t2_1" = `1`) %>% 
  rename("t2_0" = `0`) 

Actual <-tibble(Actual = c("y = 1", "y = 0")) 
Total <- tibble(Total = c(111, 62))  
  

library(officer) # for fp_border() this needs to load before flextable

library(flextable)
suppressMessages(conflict_prefer("compose", "flextable"))

library(dplyr)  # for bind_cols()


# Make analysis table
`Table 4.4` <- bind_cols(Actual, t1df, t2df, Total)

# The header needs blank columns for spaces in actual table.
# The wide column labels use Unicode characters for pi. Those details are 
#   replaced later with the compose() function.
theHeader <- data.frame(
  col_keys = c("Actual", "blank", "t1_1", "t1_0",
               "blank2", "t2_1", "t2_0", 
               "blank3", "Total"),
  line1 = c("Actual", "", 
            rep("Prediction, \U1D70B = 0.6416", 2), "", 
            rep("Prediction, \U1D70B = 0.50", 2), "",
            "Total"),
  line2 = c("Actual", "", "t1_1", "t1_0", "", "t2_1", "t2_0", "", "Total"))

# Border lines
big_border <- fp_border(color="black", width = 2)

# Make the table - compose uses Unicode character 
# https://stackoverflow.com/questions/64088118/in-r-flextable-can-complex-symbols-appear-in-column-headings
flextable(`Table 4.4`, col_keys = theHeader$col_keys) %>% 
  set_header_df(mapping = theHeader, key = "col_keys") %>% 
  theme_booktabs() %>% 
  merge_v(part = "header") %>% 
  merge_h(part = "header") %>% 
  align(align = "center", part = "all") %>% 
  empty_blanks() %>% 
  fix_border_issues()  %>%
  set_caption(caption = "Classification tables for horseshoe crab data with width and factor color predictors.") %>% 
  set_table_properties(width = 1, layout = "autofit") %>% 
  hline_top(part="header", border = big_border) %>% 
  hline_bottom(part="body", border = big_border)  %>% 
  # i = 1 refers to top row of column headerings
  compose(i = 1, j = "t1_1", part = "header",
          value = as_paragraph("Prediction, \U1D70B",  
                               as_sub("0"), "= 0.6416")) %>% 
  compose(i = 1, j = "t2_1", part = "header",
          value = as_paragraph("Prediction, \U1D70B",  
                               as_sub("0"), "= 0.50")) %>% 
  # i = 2 refers to second row of column headings
  compose(part = "header", i = 2, j = "t1_1",  
          value = as_paragraph("y\U0302 = 1")) %>% 
  compose(part = "header", i = 2, j = "t1_0", 
          value = as_paragraph("y\U0302 = 0")) %>% 
  compose(part = "header", i = 2, j = "t2_1", 
          value = as_paragraph("y\U0302 = 1")) %>% 
  compose(part = "header", i = 2, j = "t2_0", 
          value = as_paragraph("y\U0302 = 0"))
```

$$\mathrm{Sensitivity} = P(\hat y = 1 | y = 1), \ \ \ \ \  \mathrm{Sensitivity} = P(\hat y = 0 | y = 0)$$

$$P(\mathrm{correct\ classif.}) = \mathrm{sensitivity}[P(y = 1)] + \mathrm{specificity}[1 - P(y = 1)],$$ 

### 4.6.2 Summarizing Predictive Power: ROC Curves

```{r}
fit <- glm(y ~ width + factor(color), family = binomial, data = Crabs)

suppressMessages(library(pROC))
suppressMessages(
  rocPlot <- roc(y ~ fitted(fit), data = Crabs)  
)

x <- data.frame(rocPlot$sensitivities,rocPlot$specificities, rocPlot$thresholds)

thePar <- par(pty = "s")
plot.roc(rocPlot, legacy.axes = TRUE, asp = F)
par(thePar)

auc(rocPlot)

# Make a prettier ROC plot including the probability cut points
model <- data.frame(outcome = Crabs$y, 
                    predicted = predict(fit, Crabs, type = "response"))  

suppressMessages(library(plotROC))
# `d` (disease) holds the known truth
# `m` (marker) holds the predictor values 
ggplot(model, aes(d = outcome, m = predicted)) + 
  geom_roc() + 
  style_roc(xlab = "1-Specificity",
            ylab ="Sensitivity",
            minor.breaks = c(seq(0, 0.1, by = 0.02), seq(0.9, 1, by = 0.02))) + 
  ggtitle("ROC curve with probabilty cutpoints") +
  theme(aspect.ratio=1) +
  geom_abline(slope = 1, color = "grey92")
  

```

### 4.6.3 Summarizing Predictive Power: Multiple Correlation {#x4.6.3}