Logistic Regression in R

1 Logistics Regression Model

Remember from the lecture that we are fitting a regression model with a binary outcome.
As such, the model is as follows:

\[\begin{eqnarray*} \log \Big(\frac{\hat{p}}{1 - \hat{p}} \Big) & = & a + b_{1}x_1 + \ldots + b_{q}x_{q} \end{eqnarray*}\]

The left hand side of this model is the logarithm of the odds of success. It is a monotonously increasing function, which means when \(b_1\) is positive, if \(x_1\) increases, the odds of success increases. Try plot this in Desmos!
Thereby, the probability of success of \(\hat{p}\) can be written as follows:

\[\begin{eqnarray*} \hat{p} & = & \frac{\exp(a + b_{1}x_1 + \ldots + b_{q}x_{q})}{1 + \exp(a + b_{1}x_1 + \ldots + b_{q}x_{q})} \end{eqnarray*}\]

The above means once we estimate the coefficients of the model, we can estimate the probability of success of the outcome of interest.

2 Import data

Let’s revisit Alzheimer’s data set, and consider the task of building a logistic regression model with diagnosis as its response variable and variables age, education, hippocampus volume, and female as its predictors.

Code

library(tidyverse) 

alzheimer_data <- read_csv("https://raw.githubusercontent.com/COSMOS-DataScience/slides/main/data/alzheimer_data.csv")

3 Data pre-processing

3.1 Clean up the columns

Running a logistic regression model in R is pretty straightforward. Before we do that, let’s use select() function to select the interested variable (diagnosis, age, education, lhippo, rhippo, and female).

Code

alzheimer_data <- alzheimer_data %>% 
  select(diagnosis, age, educ, female, lhippo,rhippo)

glimpse(alzheimer_data)

Rows: 2,700
Columns: 6
$ diagnosis <dbl> 0, 0, 0, 0, 1, 0, 0, 2, 0, 2, 0, 0, 0, 1, 0, 1, 2, 2, 2, 1, …
$ age       <dbl> 74, 56, 77, 74, 75, 72, 64, 78, 73, 81, 66, 65, 66, 73, 78, …
$ educ      <dbl> 12, 16, 18, 20, 14, 16, 16, 17, 18, 13, 16, 16, 17, 20, 13, …
$ female    <dbl> 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, …
$ lhippo    <dbl> 2.2900, 3.2606, 2.6990, 3.0600, 2.9342, 3.2100, 3.6800, 1.73…
$ rhippo    <dbl> 2.9200, 3.3321, 2.5028, 3.0000, 3.2890, 3.0600, 3.8200, 2.30…

3.2 Data wrangling

Let’s begin by transforming the response to a new feature with two categories: no symptoms (0) versus mild or strong symptoms (1).
It is noticeable that female is a numeric variable. As such, we should make sure R recognizes that feature as a factor variable.
Create a new variable hippo as the sum of left hippocampus volume (lhippo) and right hippocampus volume (rhippo).

Code

alzheimer_data <- alzheimer_data %>% 
  mutate(diag = ifelse(diagnosis %in% c(1, 2), "1", "0"),
         diag = as.factor(diag),
         female = as.factor(female),
         hippo = lhippo + rhippo)

3.3 Data visualization

Let’s explore these relationships using visualizations.

What is the relationship between diagnosis v.s. hippocampus?
What is the relationship between diagnosis v.s. gender?
What is the relationship between diagnosis v.s. education years?

Code

alzheimer_data %>% 
  ggplot(aes(x = hippo, y = diag, color = diag))+
  geom_boxplot()+
  theme_classic() +
  labs(x = 'Hippocampus volume (cc)', y = 'Diagnosis group')+
  theme(aspect.ratio = 0.5)

Code

alzheimer_data %>% 
  mutate(female = fct_recode(female, Male="0",Female="1")) %>% 
  ggplot(aes(x = female, fill = diag))+
  geom_bar(position = 'fill') +
  theme_classic() +
  labs(x = 'Female', y = 'Diagnosis group')+
  theme(aspect.ratio = 0.5)

Code

alzheimer_data %>% 
  ggplot(aes(x = educ, y = diag, color = diag))+
  geom_boxplot()+
  theme_classic() +
  labs(x = 'Education in years', y = 'Diagnosis group')+
  theme(aspect.ratio = 0.5)

4 Model fitting

We will use glm() function to fit logistics regression.

Code

logistic_model <- glm(
  diag ~ educ + age + female + hippo, 
  family = "binomial", 
  data = alzheimer_data)

summary(logistic_model)


Call:
glm(formula = diag ~ educ + age + female + hippo, family = "binomial", 
    data = alzheimer_data)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  6.417112   0.595993  10.767  < 2e-16 ***
educ        -0.053793   0.013290  -4.048 5.17e-05 ***
age          0.017096   0.004262   4.011 6.04e-05 ***
female1     -1.347349   0.097097 -13.876  < 2e-16 ***
hippo       -1.037441   0.059413 -17.462  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3692.7  on 2699  degrees of freedom
Residual deviance: 3019.4  on 2695  degrees of freedom
AIC: 3029.4

Number of Fisher Scoring iterations: 4

5 Model interpretation:

If we plug in the estimated values from the model above, we have

\[\begin{eqnarray*} \log \Big(\frac{\hat{p}}{1- \hat{p}} \Big) & = & 6.42 -0.05\,x_{\text{educ}} + 0.02\,x_{\text{age}} -1.35\,x_{\text{female}} -1.04\,x_{\text{hippo}} \end{eqnarray*}\]

Let’s say we want to estimate the probability of having symptoms for Alzheimer’s disease for a person with specific characteristics.

What is the probability of having symptoms for Alzheimer’s disease for a 40 years old women, with 10 education year, and 6 cc of hippocampus volume?

newx_df <- data.frame(
  age    = 40,
  female = factor(0, levels = levels(alzheimer_data$female)),
  educ   = 10,
  hippo  = 6
) 
pred_prob = predict(logistic_model, newdata = newx_df, type = "response") %>%
  round(digits = 2)

\[\begin{eqnarray*} \hat{P}(\text{diag = 1}) & = & \frac{\exp(6.42 + -0.05\,\times 10 + 0.02\,\times 40 + -1.35 \times 1 -1.04 \times 6)}{1 + \exp(6.42 + -0.05\,\times 10 + 0.02\,\times 40 + -1.35 \times 1 -1.04 \times 6)} \end{eqnarray*} = 0.58\]

6 Format the output to a publication-ready table

Code

library(gtsummary) # for tbl_regression() function

logistic_model %>%
  tbl_regression(
    estimate_fun = function(x){style_number(x, digits = 3)}, 
    exponentiate = TRUE)

Characteristic	OR	95% CI	p-value
educ	0.948	0.923, 0.973	<0.001
age	1.017	1.009, 1.026	<0.001
female
0	—	—
1	0.260	0.215, 0.314	<0.001
hippo	0.354	0.315, 0.398	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

7 Model Evaluation:

To split the data into training and validation sets using the rsample package in R, you can use the initial_split() function. Here’s an example of how you can split the data:

Code

library(rsample)

set.seed(0)
data_split <- initial_split(alzheimer_data, prop = 0.7) 

train_data <- training(data_split)
test_data <- testing(data_split)

7.0.1 Fit the model with train data

As we saw, next step after splitting the data into train and test would be training the model using training data:

Code

logistic_model2 <- glm(diag ~ educ + age + hippo + female, family=binomial, data=train_data)
summary(logistic_model2)


Call:
glm(formula = diag ~ educ + age + hippo + female, family = binomial, 
    data = train_data)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  6.595232   0.722410   9.129  < 2e-16 ***
educ        -0.064738   0.015687  -4.127 3.68e-05 ***
age          0.016344   0.005117   3.194   0.0014 ** 
hippo       -1.032922   0.071457 -14.455  < 2e-16 ***
female1     -1.333364   0.116105 -11.484  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2578.5  on 1888  degrees of freedom
Residual deviance: 2112.9  on 1884  degrees of freedom
AIC: 2122.9

Number of Fisher Scoring iterations: 4

Code

logistic_model2 %>%
  tbl_regression(estimate_fun = function(x) style_number(x, digits = 3), exponentiate = TRUE)

Characteristic	OR	95% CI	p-value
educ	0.937	0.909, 0.966	<0.001
age	1.016	1.006, 1.027	0.001
hippo	0.356	0.309, 0.409	<0.001
female
0	—	—
1	0.264	0.209, 0.330	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

7.0.2 Validate on test data

Followed, by testing it via the validation set. This means to calculate the probability of success for each subject in the test set:

Code

pred_prob <- logistic_model2 %>% 
  predict(test_data,type="response")

We are now ready to calculate the accuracy of our trained model. To accomplish this, we translate all probabilities of success above to 0.5 to a 1 and otherwise to a 0, followed by tracking the number of correct predictions (1’s correctly predicted as 1’s and 0’s correctly predicted as 0’s).

Code

predicted.classes <- ifelse(pred_prob > 0.5, "1", "0")
acc = mean(predicted.classes == test_data$diag) %>% 
  round(2)

print(paste("This model yields a", acc, "accuracy rate!"))

[1] "This model yields a 0.71 accuracy rate!"

8 Let’s practice

Try to build a logistics model to predict AD (diagnosis = 2) or not AD (diagnosis = 1 or 0). In other words, transforming the response to a new feature with two categories: no or mild symptoms(0) versus strong symptoms (1).

Note: give the data/model a new name!

Question: What is your accuracy on validation set? Is it higher than chance?

alzheimer_data <- read_csv("https://raw.githubusercontent.com/COSMOS-DataScience/slides/main/data/alzheimer_data.csv")

Rows: 2700 Columns: 57
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): id
dbl (56): diagnosis, age, educ, female, height, weight, bpsys, bpdias, hrate...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

AD_data <- alzheimer_data %>% 
  mutate(diag = ifelse(diagnosis %in% c(2), "1", "0"),
         diag = as.factor(diag),
         female = as.factor(female),
         hippo = lhippo + rhippo)
 
# data splitting
set.seed(0)
data_split <- initial_split(AD_data, prop = 0.7) 

train_data <- training(data_split)
test_data <- testing(data_split)

# model fit
logistic_AD <- glm(diag ~ educ + age + hippo + female, family=binomial, data=train_data)
summary(logistic_AD)


Call:
glm(formula = diag ~ educ + age + hippo + female, family = binomial, 
    data = train_data)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.8440109  0.7767975   6.236 4.49e-10 ***
educ        -0.0482131  0.0170757  -2.823  0.00475 ** 
age         -0.0004922  0.0060436  -0.081  0.93509    
hippo       -0.8407063  0.0737020 -11.407  < 2e-16 ***
female1     -0.9188226  0.1287350  -7.137 9.52e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1874.3  on 1888  degrees of freedom
Residual deviance: 1663.3  on 1884  degrees of freedom
AIC: 1673.3

Number of Fisher Scoring iterations: 4

# prediction
pred_prob <- logistic_AD %>% 
  predict(test_data,type="response")

predicted.classes <- ifelse(pred_prob > 0.5, "1", "0")
acc = mean(predicted.classes == test_data$diag) %>% 
  round(2)

print(paste("This model yields a", acc, "accuracy rate!"))

[1] "This model yields a 0.8 accuracy rate!"

# chance level?
chance1 = mean(test_data$diag == sample(test_data$diag, size = nrow(test_data), replace = TRUE)) %>% 
  round(2)
chance2 = max(table(test_data$diag))/nrow(test_data) %>% 
  round(2)