• We will use functions from packages base, utils, and stats (pre-installed and pre-loaded)
• We may also use the packages below (specifying package::function for clarity).

# Load pckgs for this R session
options(scipen = 999)
# --- General 
library(here)     # tools find your project's files, based on working directory
library(dplyr)    # A Grammar of Data Manipulation
library(skimr)    # Compact and Flexible Summaries of Data
library(magrittr) # A Forward-Pipe Operator for R 
library(readr)    # A Forward-Pipe Operator for R 
library(tidyr)    # Tidy Messy Data
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax

# ---Plotting & data visualization
library(ggplot2)      # Create Elegant Data Visualisations Using the Grammar of Graphics
library(ggfortify)     # Data Visualization Tools for Statistical Analysis Results
library(scatterplot3d) # 3D Scatter Plot

# --- Statistics
library(MASS)       # Support Functions and Datasets for Venables and Ripley's MASS
library(factoextra) # Extract and Visualize the Results of Multivariate Data Analyses
library(FactoMineR) # Multivariate Exploratory Data Analysis and Data Mining
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests
library(car)        # Companion to Applied Regression
library(ROCR)       # Visualizing the Performance of Scoring Classifiers

# --- Tidymodels (meta package)
library(rsample)    # General Resampling Infrastructure  
library(broom)      # Convert Statistical Objects into Tidy Tibbles

• Logistic regression is a classification model used with a binary response variable, e.g.:
  • yes|no, or 0|1, or True|False in a survey question;
  • success|failure in a clinical trial experiment;
  • benign|malignant in a biopsy experiment.
• Logistic regression is a type of Generalized Linear Model (GLM): a more flexible version of linear regression that can work also for categorical response variables or count data (e.g. poisson regression).
• When logistic regression fits the coefficients \(\beta_0\), \(\beta_1\), …, \(\beta_k\) to the data, it minimizes errors using the Maximum Likelihood Estimation method (as opposed to linear regression’s Least Squares Error Estimation method).

If we have predictor variables like \(X_{1}\), \(X_{2}\), …, \(X_{k}\) and a binary response variable \(Y\) (where \(y_i = 0\) or \(y_i = 1\)), we need a “special” function to transform the expected value of the response variable into the \([0,1]\) outcome we’re trying to predict.

The logit function (of the GLM family) helps us determine the coefficients \(\beta_0\), \(\beta_1\), …, \(\beta_k\) that best fit this sort of data and it is defined as: \[ logit (p_i) = \ln\left( \frac{p_i}{1-p_i} \right)= \beta_0 + \beta_1 x_{1,i} + \cdots + \beta_k x_{k,i} \] where \(p_i\) is the probability that \(y_i = 1\)

Via the inverse-logit function (logistic), we solve for the probability \(p_i\), given values of the predictor variables, like so:

\[ P(y_i = 1 | x_{1,i}, \ldots, x_{k,i} ) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_{1,i} + \cdots + \beta_k x_{k,i})}} \]

Let’s see also the relationship between the binary outcome variable heart_disease and the binary predictor variable coffee_drinker.

• we use the handy count function from dplyr to count occurrences in categorical variable(s) combinations.

# Dataset manipulation
heart_train %>% 
  # count the unique values per each group from 2 categorical variables' combinations
  dplyr::count(heart_disease, coffee_drinker, name = "count_by_group") %>% 
  dplyr::group_by(coffee_drinker) %>% 
  dplyr::mutate(
    total_coffee_class = sum(count_by_group),
    proportion = count_by_group / total_coffee_class) %>% 
  dplyr::ungroup() %>% 
  # filter only those with heart disease
  dplyr::filter(heart_disease == "Yes_HD") %>% 
# Plot  
  ggplot(aes(x = coffee_drinker, y = proportion, fill = coffee_drinker)) + 
  geom_bar(stat = "identity") + 
  scale_fill_manual(values = c("#57b7ff", "#e07689")) + 
  theme_minimal() +
  ylab("Percent with Heart Disease") +
  xlab("Coffee Drinker (Y/N)") +
  ggtitle("Figure 3: Percent of Coffee Drinkers with Heart Disease") +
  labs(fill = "Coffee Drinker") + 
  scale_y_continuous(labels = scales::percent)

• In principle, we could use a linear regression model to study likelihood of having Heart Disease in relation to risk factors: \[ Y = \beta_0 + \beta_1\text{(US\$ Spent Fast Food)} + \beta_2\text{(Coffee Drinker = YES)} \]

• But we’ll see why the logistic regression model is a better option.

\[ \log \left( \frac{p}{1-p} \right) = \beta_0 + \beta_1\text{(US\$ Spent Fast Food)} + \beta_2\text{(Coffee Drinker = YES)} \]

• Let’s compare the 2 models (for simplicity, we ignore the coffee_drinker variable for now.):

# --- 1) Linear regression model
linear_mod <- lm(heart_disease ~ fast_food_spend# + coffee_drinker
                 , data = heart_data)
# --- 2) Logistic regression model
logit_mod <- glm(heart_disease ~ fast_food_spend# + coffee_drinker
                 , data = heart_data, family = "binomial")

We can now extract the coefficients from the linear regression and logistic models, then estimate the predicted outcomes from these models: lin_pred, logit_pred, and logistic_pred (i.e. the conversion form log(odds) to probability).

# --- 1) Extract coefficients from linear regression model
intercept_lin <- coef(linear_mod)[1]
fast_food_spend_lin <- coef(linear_mod)[2]
coffee_drinker_lin <- coef(linear_mod)[3]

# --- 2) Extract coefficients from logit regression model
intercept_logit <- coef(logit_mod)[1]
fast_food_spend_logit <- coef(logit_mod)[2]
coffee_drinker_logit <- coef(logit_mod)[3]

# --- Estimate predicted data from different models   
heart_data <- heart_data %>%
  mutate(
    # Convert outcome variable to factor
    heart_disease_factor = factor(heart_disease, 
                                  labels = c("No Disease (Y=0)", "Disease (Y=1)")),
    # 1) Linear model prediction
    lin_pred = intercept_lin + fast_food_spend_lin * fast_food_spend, 
    # 2) Logit model prediction
    logit_pred = intercept_logit + fast_food_spend_logit * fast_food_spend,  coffee_drinker,
    # 3) Convert logit to probability (logistic model prediction)
    logistic_pred = 1 / (1 + exp(-logit_pred)) 
    ) %>%
  arrange(fast_food_spend)   

• Besides the poor fit, recall linear regression models implies certain assumptions:
  • Linear relationship between Y and the predictors X1 and X2
    
  • Residuals must be 1) approximately normally distributed and 2) uncorrelated
  • Homoscedasticity: residuals should have constant variance
  • Non collinearity: predictor variables should not be highly correlated with each other

• Assumptions that are not met in this case

# Set graphical parameters to plot side by side
par(mfrow = c(1, 2))

# Diagnostic plots
plot(linear_mod, which = 1)
plot(linear_mod, which = 2)

# Reset graphical parameters (optional, avoids affecting later plots)
par(mfrow = c(1, 1))

• Table 1 (next) shows the model output, with the coefficient estimate for each predictor.
  • The broom::tidy function converts the model summary into a more readable data frame.
  • The odds ratio (= exponentiated coefficient estimate) is more interpretable than the coefficient itself.

# Convert model's output summary into data frame
heart_model_coef <- broom::tidy(heart_model) %>% 
  # improve readability of significance levels
  dplyr::mutate('signif. lev.' = case_when(
    `p.value` < 0.001 ~ "***",
    `p.value` < 0.01 ~ "**",
    `p.value` < 0.05 ~ "*",
    TRUE ~ ""))%>%
  # add odds ratio column
  dplyr::mutate(odds_ratio = exp(estimate)) %>%
  dplyr::relocate(odds_ratio, .after = estimate) %>%
  dplyr::mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  # format as table
  knitr::kable() %>% 
  # reduce font size
  kable_styling(font_size = 20) %>% 
  # add table title
  kableExtra::add_header_above(c("Logistic Regression Analysis of Heart Disease Risk Factors"= 7))

heart_model_coef

• Intercept: gives the log-odds of heart disease when all predictor variables are zero. This is not generally interpreted, but the highly negative value suggests very low probability of heart disease in the sample of reference.
  • (If interpreted) the intercept term -11.0554 would translate as probability \(P = e^{-11.0554} / (1 + e^{-11.0554}) = 0.00002\) or \(0.002\%\).
  • …which means: when \(\text{coffee_drinker} = NO\), and \(\text{fast_food_spend} = 0\), and \(\text{income} = 0\), the probability of heart disease is as low as \(0.002\%\).

• Income: Based on a p-value = 0.8182, we conclude that income is not significantly associated with heart disease. (Anyhow, the odds ratio of ≈1 would suggest no change in odds based on income.)

• The positive coefficient of Fast Food Annual Spending (US$), 0.0024 (highly statistically significant as p-value = 0.0000), suggests a positive association with heart disease.

• This means that for each additional \(\Delta X_1 = +1 \; US\$\) spent on fast food annually:

  • the log(odds) of heart disease increases by: \(\Delta \log(\text{odds}) = \beta_{1} \times \Delta X_1 = 0.0024 \times 1 = 0.0024\)

  • the odds of heart disease (with spending) increase by a factor of: \(OR = e^{0.0024} \approx 1.0024\) (compared to without spending).

• The probability of heart disease can be computed using the logistic function: \(P_{HD=1} = \frac{e^{\beta_0 + (\beta_1 \times X_1)}}{1 + e^{\beta_0 + (\beta_1 \times X_1)}}\)
  • Solving for \(\beta_0 = -11.0554\), \(\beta_1 = 0.0024\) and \(X_1 = 1\) gives:
  \[P_{HD=1} = \frac{e^{-11.0554 + (0.0024 \times 1)}}{1 + e^{-11.0554 + (0.0024 \times 1)}} = 0.00001583981\]

which indicates a probability of heart disease is as low as \(0.00159\%\) when fast_food_spend \(= 1\; US\$\).

• However, considering that this is annual spending, we should use a more adequate scale!

• For example, for an additional \(\Delta X_1 = +100 \; US\$\) spent on fast food annually:

  • the log(odds) of heart disease increases by: \(\Delta \log(\text{odds}) = \beta_{1} \times \Delta X_1 = 0.0024 \times 100 = 0.24\)

  • the odds of heart disease (with spending) increase by a factor of: \(OR = e^{0.24} \approx 1.271\) (compared to without spending).

• The probability of heart disease can be computed using the logistic function: \(P_{HD=1} = \frac{e^{\beta_0 + (\beta_1 \times X_1)}}{1 + e^{\beta_0 + (\beta_1 \times X_1)}}\)

• Solving for \(\beta_0 = -11.0554\), \(\beta_1 = 0.0024\) and and \(X_1 = 100\) gives:
  \[P_{HD=1} = \frac{e^{-11.0554 + (0.0024 \times 100)}}{1 + e^{-11.0554 + (0.0024 \times 100)}} \approx 0.00002008\]

  which (this time) indicates a (still very low probability) of heart disease at \(0.002008\%\) when fast_food_spend \(= 100\; US\$\).
  

• The negative coefficient of Coffee Drinker (=YES), -0.7296 expressed in log-odds, means that coffee drinkers have lower odds of having heart disease (relative to non-coffee drinkers).
• This effect is statistically significant as p-value = 0.0122.
• Transforming the estimate into odds ratio, we obtain \(\text{O.R.} = e^{-0.7296} = 0.48\) , meaning coffee-drinkers have 0.48 (or 48%) the odds that non-coffee drinkers have to experience heart disease (holding all other predictors fixed). \[ \text{Odds Ratio} = e^{-0.7296} = 0.48 \]
  • How much is the drop of the odds? Since O.R. is < 1, I have a percentage decrease:
    \((1 - 0.48) = 0.52\) or \(52\%\) lower odds of having heart disease compared to non-coffee drinkers.

• Odds represent the ratio of the probability of an event occurring to the probability of it not occurring:
  • Example: If the probability of heart disease is 0.25 (25%), then the odds are:

\[ \text{Odds} = \frac{P(\text{event})}{1 - P(\text{event})} \]

\[ \frac{0.25}{1 - 0.25} = \frac{0.25}{0.75} = 0.33 \]

This means that for every 1 person without heart disease, there are 0.33 people with it.
(Or equivalently: for every 3 people without heart disease, 1 person has it.)

• Odds Ratio (O.R.) compares the odds of an event occurring in one group (exposed) relative to another group (reference).

  • Example: If the odds of heart disease for coffee drinkers are 0.48, and for non-coffee drinkers are 1, then the odds ratio is:

\[\text{Odds Ratio} = \frac{\text{Odds}_{\text{exposed}}}{\text{Odds}_{\text{reference}}} = \frac{\text{Odds}_{\text{coffee drinkers}}}{\text{Odds}_{\text{non-coffee drinkers}}} = \frac{0.48}{1} = 0.48\]

• Interpretation:
  • If O.R. > 1, the event is more likely in the exposed group.
  • If O.R. < 1, the event is less likely in the exposed group.
  • If O.R. = 1, there is no difference between groups.

If the odds of heart disease are:

• Coffee drinkers (exposed group): \(0.48\)
• Non-coffee drinkers (reference group): \(1\)

Then the Odds Ratio is:

\[\text{Odds Ratio} = \frac{\text{Odds}_{\text{coffee drinkers}}}{\text{Odds}_{\text{non-coffee drinkers}}} = \frac{0.48}{1} = 0.48\]

This means coffee drinkers have 52% lower odds of heart disease than non-coffee drinkers.

Let’s use the logistic regression model to make predictions on the training data. First, we’ll do it manually:

• Consider an individual that drinks coffee, spends US$5,000 on fast food per year, and has an annual income of US$50,000.

• We can calculate the probability of this individual having heart disease using our logistic regression model:

\[p_i(HD=YES) = \frac{1}{1 + e^{-(\beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3)}} = \]

\[ = \frac{1}{1 + e^{-(-11.0554 -0.7296 \times 1 + 0.0024 \times 5000 -0.0000023 \times 50,000 )}} = 0.4951735\]

(decimals count!)

Now, let’s use the predict() function to make a prediction for the same individual.

individual <- data.frame(
  coffee_drinker = "Yes_Coffee", # factor 
  fast_food_spend = 5000, 
  income = 50000)

# Make a prediction
predict(heart_model, individual, type = "response")

        1 
0.4951735 

In fact, we can use the predict() function to make predictions for all the individuals observed in the training data heart_train.

• heart_train$heart_disease_pred = added column with the predicted probabilities of heart disease for each individual.

# Add calculated predictions to the training data
heart_train$heart_disease_pred <- predict(heart_model, type = "response")

# Subset  of 3 with heart_disease = Yes_HD
heart_train[heart_train$heart_disease == "Yes_HD", ][1:3,] 

    heart_disease coffee_drinker fast_food_spend   income heart_disease_pred
207        Yes_HD      No_Coffee        4723.998 48956.17         0.51420472
210        Yes_HD     Yes_Coffee        4748.477 20655.20         0.36601782
242        Yes_HD     Yes_Coffee        3932.141 14930.18         0.07756531

# Subset  of 3 with heart_disease = No_HD
heart_train[heart_train$heart_disease == "No_HD", ][1:3,]

  heart_disease coffee_drinker fast_food_spend    income heart_disease_pred
1         No_HD      No_Coffee        1823.816 44361.625        0.001086288
3         No_HD      No_Coffee        2683.873 31767.139        0.008566989
6         No_HD     Yes_Coffee        2298.971  7491.559        0.001762176

A 3D scatterplot of observations shows the first 3 principal components’ scores.

• For this one, we need the scatterplot3d() function of the scatterplot3d package;
• The color argument assigned to the Label variable;
• To add a legend, we use the legend() function and specify its coordinates via the xyz.convert() function.

# 3D scatterplot ...
plot_3d <- with(PC_scores, 
                scatterplot3d::scatterplot3d(PC_scores$PC1, 
                                             PC_scores$PC2, 
                                             PC_scores$PC3, 
                                             color = as.numeric(Label), 
                                             pch = 19, 
                                             main ="Figure 3: 3D Scatter Plot", 
                                             xlab="PC1",
                                             ylab="PC2",
                                             zlab="PC3"))

# ... + legend
legend(plot_3d$xyz.convert(0.5, 0.7, 0.5), 
       pch = 19, 
       yjust=-0.6,
       xjust=-0.9,
       legend = levels(PC_scores$Label), 
       col = seq_along(levels(PC_scores$Label)))

Biplots have two key elements: scores (the 2 axes) and loadings (the vectors). As in the scores plot, each point represents an observation projected in the space of principal components where:

• Biopsies of the same class are located closer to each other, which indicates that they have similar scores referred to the 2 main principal components;
• The loading vectors show strength and direction of association of original variables with new PC variables.

As expected from PCA, the single PC1 accounts for variance in almost all original variables, while V9 has the major projection along PC2.