Linear versus logistic regression

We know that a linear regression is a convenient way to estimate the relationship between a quantitative continuous dependent variable, and one or more independent variables (of any type).

For instance, suppose we would like to estimate the relationship between two quantitative variables, \(X\) and \(Y\). Using the ordinary least squares method (the most common estimator used in linear regression), we obtain the following regression line:

Now suppose we are interested in estimating the impact of age on whether or not a patient has a certain disease. Age is considered as a quantitative continuous variable, while having the disease is binary (a patient is either ill or healthy).

Visually, we could have something like this:

If we fit a regression line (using the ordinary least square method) to the points, we obtain the following plot:

We see that the regression line goes below 0 and above 1 with respect to the \(y\)-axis. Since the dependent variable disease cannot take values below 0 (= healthy) nor above 1 (= ill), it is obvious that a linear regression is not appropriate for these data!

In addition to this limitation, the assumptions of normality and homoscedasticity, which are required in linear regression, are clearly not appropriate with these data since the dependent variable is binary and follows a Binomial distribution! R will not stop you from performing a linear regression on binary data, but this will produce a model of little interest.

This is where a logistic regression becomes handy as it takes into consideration these limitations.

Applied to our example, here is how the points are fitted using a binary logistic regression:

It is clear that this model is more appropriate.

The curve (known as a sigmoid) is obtained via a transformation of the predicted values. There are several possible choices for the link function, which aim is to constrain predicted values to be within the range of observed values.

The most widely used in practice is the logit function, which relates the probability of occurrence of an event (bounded between 0 and 1) to the linear combination of independent variables. The logit function also turns out to be the canonical link function for a Bernoulli or Binomial distribution. This transformation ensures that no matter in which range the \(X\) values are located, \(Y\) will only take numbers between 0 and 1.

One could say that the fitted values (= represented by the blue curve) taking values between 0 and 1 also does not seem to make sense since a patient can only be healthy or ill (and thus the dependent variable can only take the value 0 or 1, respectively). However, in a binary logistic regression it is not the outcome no disease/disease that is directly modeled, but the likelihood that a patient has the disease or not given his or her characteristics. This likelihood will be framed in terms of a probability to observe or not the disease in a patient, which is indeed included between 0 and 1, or 0% and 100%.

More generally, with a logistic regression we would like to model how the probability of success varies with the independent variables and determine whether or not these changes are statistically significant. We are actually going to model the logarithm of the odds, and the logistic regression model will be written as follows:

\[\begin{align} \log(odds(success)) &= logit(\pi) \\ &= \log\left(\frac{\pi}{1 - \pi}\right) \\ &= \beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p \end{align}\]

where \(\pi\) \((0 \le \pi \le 1)\) is the probability of an event happening (success) and denoted \(\pi = P(success)\). We find the values of \(\hat{\beta}_0\), \(\hat{\beta}_1\), \(\ldots\), \(\hat{\beta}_p\), which are used as estimates for \(\beta_0\), \(\beta_1\), \(\ldots\), \(\beta_p\), using the maximum likelihood method. This method is one of several methods used in statistics to estimate parameters of a mathematical model. The goal of the estimator is to estimate the parameters \(\beta_0\), \(\beta_1\), \(\ldots\), \(\beta_p\) which maximize the log likelihood function. Different algorithms have been established over the years for this non-linear optimization, but this is beyond the scope of the post.

Logistic regressions are very common in the medical field, for example to:

• estimate the risk factors associated with a disease or a harmful condition,
• predict the risk of developing a disease based on a patient’s characteristics, or
• determine the most important biological factors associated with a specific disease or condition.

However, logistic regressions are used in many other domains, for instance in:

• banking sector: estimate a debtor’s creditworthiness based on his or her profile (income, assets, liabilities, etc.),
• marketing: estimate a customer’s propensity to buy a product or service based on his or her profile (age, sex, salary, previous purchases, etc.),
• sports: estimate the probability of a player winning against another player as a function of the characteristics of the two opponents,
• politics: answer the question “Would a citizen vote for our political party at the next elections?”
• etc.

Univariable versus multivariable logistic regression

Now that it is more clear when a binary logistic regression should be used, we show how to perform one in R. We start by presenting univariable binary logistic regressions, and then multivariable binary logistic regressions.

Remember that in both cases, the dependent variable must be a qualitative variable with two outcomes (hence the name binary logistic regression). The difference between a univariable and multivariable binary logistic regression lies in the fact that:

• for a univariable binary logistic regression, there is only one independent variable, while
• for a multivariable binary logistic regression, there are two ore more independent variables.

It is true that the term “univariable” may be confusing here because there are two variables in the model (i.e., one dependent variable and one independent variable). However, it is called univariable binary logistic regression to indicate that only one independent variable is considered in the model, as opposed to multivariable binary logistic regression where several independent variables are considered in the model.

To draw a parallel with linear regression:

• a univariable binary logistic regression is the equivalent of a simple linear regression, whereas
• a multivariable binary logistic regression is the equivalent of a multiple linear regression

when the dependent variable is binary instead of quantitative continuous. This is the reason a univariable binary logistic regression is sometimes called a simple binary logistic regression and a multivariable binary logistic regression sometimes called a multiple binary logistic regression.

The terms univariable/multivariable should not be confused with univariate/multivariate. The number of dependent variables characterizes the model as univariate or multivariate; univariate refers to a model with only one dependent variable, while multivariate refers to a model that simultaneously predicts more than one dependent variable. Usually, the intent is to differentiate models based on the number of independent variables. This distinction is made thanks to the terms univariable and multivariable. Multivariable refers to a model relating multiple predictor variables to a dependent variable, whereas univariable refers to a model relating one single independent variable to a dependent variable.

Univariable binary logistic regression

As mentioned above, we start with a univariable binary logistic regression, that is, a binary logistic regression with only one independent variable.

In R, a binary logistic regression can be done with the glm() function and the family = "binomial" argument. Similar to linear regression, the formula used inside the function must be written as dependent variable ~ independent variable (in this order!).

While the dependent variable must be categorical with two levels, the independent variable can be of any type. However, interpretations differ depending on whether the independent variable is qualitative or quantitative.

For completeness, we illustrate this type of regression with both a quantitative and a qualitative independent variable, starting with a quantitative independent variable.

# save model
m1 <- glm(heart_disease ~ age,
  data = dat,
  family = "binomial"
)

# print results
summary(m1)

## 
## Call:
## glm(formula = heart_disease ~ age, family = "binomial", data = dat)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.05122    0.76862  -3.970  7.2e-05 ***
## age          0.05291    0.01382   3.829 0.000128 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 409.95  on 296  degrees of freedom
## Residual deviance: 394.25  on 295  degrees of freedom
## AIC: 398.25
## 
## Number of Fisher Scoring iterations: 4

# OR for age
exp(coef(m1)["age"])

##      age 
## 1.054331

# prob(heart disease) for age = 0
exp(coef(m1)[1]) / (1 + exp(coef(m1)[1]))

## (Intercept) 
##  0.04516478

# 95% CI for the OR for age
exp(confint(m1,
  parm = "age"
))

##    2.5 %   97.5 % 
## 1.026699 1.083987

# predict probability to develop heart disease
pred <- predict(m1,
  newdata = data.frame(age = c(30)),
  type = "response"
)

# print prediction
pred

##         1 
## 0.1878525

# predict probability to develop heart disease
pred <- predict(m1,
  newdata = data.frame(age = c(30)),
  type = "response",
  se = TRUE
)

# print prediction
pred$fit

##         1 
## 0.1878525

# 95% confidence interval for the prediction
lower <- pred$fit - (qnorm(0.975) * pred$se.fit)
upper <- pred$fit + (qnorm(0.975) * pred$se.fit)
c(lower, upper)

##          1          1 
## 0.07873357 0.29697138

# load package
library(sjPlot)

# plot
plot_model(m1,
  type = "pred",
  terms = "age"
) +
  labs(y = "Prob(heart disease)")

# levels for sex
levels(dat$sex)

## [1] "female" "male"

# save model
m2 <- glm(heart_disease ~ sex,
  data = dat,
  family = "binomial"
)

# print results
summary(m2)

## 
## Call:
## glm(formula = heart_disease ~ sex, family = "binomial", data = dat)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0438     0.2326  -4.488 7.18e-06 ***
## sexmale       1.2737     0.2725   4.674 2.95e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 409.95  on 296  degrees of freedom
## Residual deviance: 386.12  on 295  degrees of freedom
## AIC: 390.12
## 
## Number of Fisher Scoring iterations: 4

# OR for sex
exp(coef(m2)["sexmale"])

##  sexmale 
## 3.573933

# prob(disease) for sex = female
exp(coef(m2)[1]) / (1 + exp(coef(m2)[1]))

## (Intercept) 
##   0.2604167

chisq.test(table(dat$heart_disease, dat$sex))

## 
## 	Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(dat$heart_disease, dat$sex)
## X-squared = 21.852, df = 1, p-value = 2.946e-06

# predict probability to develop heart disease
pred <- predict(m2,
  newdata = data.frame(sex = c("male")),
  type = "response"
)

# print prediction
pred

##         1 
## 0.5572139

# plot
plot_model(m2,
  type = "pred",
  terms = "sex"
) +
  labs(y = "Prob(heart disease)")

Multivariable binary logistic regression

The interpretation of the coefficients in multivariable logistic regression is similar to the interpretation in univariable regression, except that this time it estimates the multiplicative change in the odds in favor of \(Y = 1\) when \(X\) increases by 1 unit, while the other independent variables remain unchanged.

This is similar to multiple linear regression, where a coefficient gives the expected change of \(Y\) for an increase of 1 unit of \(X\), while keeping all other variables constant.

The main advantages of using a multivariable logistic regression compared to a univariable logistic regression are to:

• consider the simultaneous (rather than isolated) effect of independent variables,
• take into account potential confounding and/or interaction effects, and
• improve predictions.

For this illustration, suppose we would like to estimate the relationship between heart disease and all variables present in the data frame, that is, age, sex, chest pain type and maximum heart rate achieved:

# save model
m3 <- glm(heart_disease ~ .,
  data = dat,
  family = "binomial"
)

# print results
summary(m3)

## 
## Call:
## glm(formula = heart_disease ~ ., family = "binomial", data = dat)
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -0.060150   1.962091  -0.031 0.975544    
## age                         0.042814   0.019009   2.252 0.024302 *  
## sexmale                     1.686330   0.349352   4.827 1.39e-06 ***
## chest_painatypical angina  -0.120481   0.641396  -0.188 0.851000    
## chest_painnon-anginal pain -0.124331   0.571093  -0.218 0.827658    
## chest_painasymptomatic      1.963723   0.548877   3.578 0.000347 ***
## max_heartrate              -0.030326   0.007975  -3.802 0.000143 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 409.95  on 296  degrees of freedom
## Residual deviance: 275.26  on 290  degrees of freedom
## AIC: 289.26
## 
## Number of Fisher Scoring iterations: 5

Note that the formula heart_disease ~ . is a shortcut to include all variables present in the data frame in the model as independent variables, except heart_disease.

Based on the \(p\)-values displayed in the last column of the coefficients table, we conclude that, at the 5% significance level, age, sex and maximum heart rate achieved are all significantly associated with heart disease (\(p\)-values < 0.05).

For the variables age, sex and max_heartrate, there is only one \(p\)-value (the result of the test on the nullity of the coefficient).

For the variable chest_pain, 3 \(p\)-values are displayed. This is normal: similar to linear regression when a categorical variable with more than two levels is included in the model, one test is performed for each comparison between the reference level and the other levels.

In our case, the reference level for the variable chest_pain is typical angina as it comes first:

levels(dat$chest_pain)

## [1] "typical angina"   "atypical angina"  "non-anginal pain" "asymptomatic"

Therefore, a test is performed for the comparison between:

1. typical angina and atypical angina,
2. typical angina and non-anginal pain, and
3. typical angina and asymptomatic.

But here we are not interested in comparing levels of the variable chest pain, we would like to test the overall effect of chest pain on heart disease. For this, we are going to compare two models via a likelihood ratio test (LRT):

1. a model which includes all the variables of interest and the variable chest_pain, and
2. the exact same model but which excludes the variable chest_pain.

The first one is referred as the full or complete model, whereas the second one is referred as the reduced model.

We compare these two models with the anova() function:

# save reduced model
m3_reduced <- glm(heart_disease ~ age + sex + max_heartrate,
  data = dat,
  family = "binomial"
)

# compare reduced with full model
anova(m3_reduced, m3,
  test = "LRT"
)

## Analysis of Deviance Table
## 
## Model 1: heart_disease ~ age + sex + max_heartrate
## Model 2: heart_disease ~ age + sex + chest_pain + max_heartrate
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       293     325.12                          
## 2       290     275.26  3    49.86 8.558e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notice that the reduced model must come before the complete model in the anova() function. The null hypothesis of this test is that the two models are equivalent.

At the 5% significance level (see the \(p\)-value at the right of the R output), we reject the null hypothesis and we conclude that the complete model is significantly better than the reduced model at explaining the presence of heart disease. This means that chest pain is significantly associated with heart disease (which was expected since the comparison between typical angina and asymptomatic was found to be significant).

A comparison of two models via the LRT will be shown again later, when discussing about interactions.

Now that we have shown that all four independent variables were significantly associated with the presence of heart disease, we can interpret the coefficients in order to know the direction of the relationships and most importantly, quantify the strength of these relationships.

Like univariable binary logistic regression, it is easier to interpret these relationships through OR. But this time, we also print the 95% CI of the OR in addition to the OR (rounded to 3 decimals) so that we can easily see which ones are significantly different from 1:

# OR and 95% CI
round(exp(cbind(OR = coef(m3), confint(m3))), 3)

##                               OR 2.5 % 97.5 %
## (Intercept)                0.942 0.020 44.353
## age                        1.044 1.006  1.084
## sexmale                    5.400 2.776 10.971
## chest_painatypical angina  0.886 0.252  3.191
## chest_painnon-anginal pain 0.883 0.293  2.814
## chest_painasymptomatic     7.126 2.509 22.030
## max_heartrate              0.970 0.955  0.985

From the OR and their 95% CI computed above, we conclude that:

• Age: the odds of having a heart disease are multiplied by a factor of 1.04 for each one-unit increase in age, all else being equal.
• Sex: the odds of developing a heart disease for males are 5.4 times the odds for females, all else being equal.
• Chest pain: the odds of developing a heart disease for people suffering from chest pain of the type “asymptomatic” are 7.13 times the odds for people suffering from chest pain of the type “typical angina”, all else being equal. We refrain from interpreting the other comparisons as they are not significant at the 5% significance level (1 is included in their 95% CI).
• Maximum heart rate achieved: the odds of having a heart disease are multiplied by a factor of 0.97 for each one-unit increase in maximum heart rate achieved, all else being equal.
• Intercept: we also refrain from interpreting the intercept as it is not significantly different from 0 at the 5% significance level (1 is included in the 95% CI).

If you are interested in printing only the OR for which the \(p\)-value of the coefficient is < 0.05, here is the code:

exp(coef(m3))[coef(summary(m3))[, "Pr(>|z|)"] < 0.05]

##                    age                sexmale chest_painasymptomatic 
##              1.0437437              5.3996293              7.1258054 
##          max_heartrate 
##              0.9701289

Remember that we can always write the interpretations in terms of probabilities with the formula \((OR - 1) \times 100\), where OR corresponds to the odds ratio.

For instance, for the maximum heart rate achieved, the OR = 0.97, so the interpretation becomes: the chance of developing a heart disease increases by (0.97 \(- 1) \times 100 =\) -3% for each one-unit increase in maximum heart rate achieved, which is equivalent to say that the chance of developing a heart disease decreases by 3% for each one-unit increase in maximum heart rate achieved.

For illustrative purposes, suppose now that we would like to predict the probability that a new patient develops a heart disease. Suppose that this patient is a 32-year-old woman, suffering from chest pain of the type non-anginal and she achieved a maximum heart rate of 150. The probability that she develops a heart disease is:

# create data frame of new patient
new_patient <- data.frame(
  age = 32,
  sex = "female",
  chest_pain = "non-anginal pain",
  max_heartrate = 150
)

# predict probability to develop heart disease
pred <- predict(m3,
  newdata = new_patient,
  type = "response"
)

# print prediction
pred

##          1 
## 0.03345948

If we trust our model, the probability that this new patient will develop a heart disease is predicted to be 3.35%.

We can also visualize the results thanks to the plot_model() function, three effects at the same time:

1. effect of age, sex and chest pain type on the predicted probability of developing a heart disease, and
2. effect of maximum heart rate achieved, sex and chest pain type on the predicted probability of developing a heart disease.

# 1. age, sex and chest pain on prob of disease
plot_model(m3,
  type = "pred",
  terms = c("age", "chest_pain", "sex"),
  ci.lvl = NA # remove confidence bands
) +
  labs(y = "Prob(heart disease)")

# 2. max heart rate, chest pain and sex on prob of disease
plot_model(m3,
  type = "pred",
  terms = c("max_heartrate", "chest_pain", "sex"),
  ci.lvl = NA # remove confidence bands
) +
  labs(y = "Prob(heart disease)")

For more clarity in the plots, confidence bands are removed thanks to ci.lvl = NA.

These plots confirm results obtained above, that is:

• there is a positive relationship between age and the presence of heart disease,
• there is a negative relationship between maximum heart rate achieved and the presence of heart disease,
• the chance of developing a heart disease is higher for patients suffering from chest pain of the type asymptomatic and similar for the 3 other types of chest pain, and
• the chance of developing a heart disease is higher for males than for females.

Interaction

In the previous sections, potential interaction effects were omitted.

An interaction occurs when the relationship between an independent variable and the outcome variable depends on the value or the level taken by another independent variable. On the contrary, if the relationship between an independent variable and the dependent variable remains unchanged no matter the value taken by another independent variable, we cannot conclude that there is an interaction effect.

In our case, there would be an interaction if for example the relationship between age and heart disease depends on the sex. There would be an interaction, for instance, if the relationship between age and heart disease was positive for females, and negative for males, or vice versa. Or if the relationship between age and heart disease was much stronger or much weaker for females than for males.

Let’s see if there is an interaction between age and sex, and more importantly, whether or not this interaction is significant. For this, we need to build two models:

1. one model containing only the main effects, so without the interaction, and
2. one model containing the main effects and the interaction.

# save model without interaction
m4 <- glm(heart_disease ~ age + sex,
  data = dat,
  family = "binomial"
)

# save model with interaction
m4_inter <- glm(heart_disease ~ age * sex,
  data = dat,
  family = "binomial"
)

We first assess the interaction visually via the plot_model() function:

# plot
plot_model(m4_inter,
  type = "pred",
  terms = c("age", "sex"),
  ci.lvl = NA # remove confidence bands
) +
  labs(y = "Prob(heart disease)")

Since the two curves of the predicted probabilities are relatively similar and follow the same pattern, the relationship between age and the presence of heart disease does not seem to depend on the sex, indicating that there may indeed be no interaction. However, we would like to test it more formally via a statistical test.

For this, we can compare the two models (the one without compared to the one with interaction) with a likelihood ratio test (LRT), using the anova() function:

anova(m4, m4_inter,
  test = "LRT"
)

## Analysis of Deviance Table
## 
## Model 1: heart_disease ~ age + sex
## Model 2: heart_disease ~ age * sex
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       294     364.43                     
## 2       293     364.23  1  0.20741   0.6488

Remember that it is always the reduced model as the first argument in the anova() function, and then the more complex model as the second argument.

The test confirms what we supposed based on the plot: at the 5% significance level, we do not reject the null hypothesis that the two models are equivalent. Since the only difference between the two models is that an interaction term is added in the complete model, we do not reject the hypothesis that there is no interaction between the age and the sex (\(p\)-value = 0.649).

This conclusion could have also been obtained more simply with the drop1() function:

drop1(m4_inter,
  test = "LRT"
)

## Single term deletions
## 
## Model:
## heart_disease ~ age * sex
##         Df Deviance    AIC     LRT Pr(>Chi)
## <none>       364.23 372.23                 
## age:sex  1   364.43 370.43 0.20741   0.6488

As you can see, this gives the same \(p\)-value of 0.649.

In practice, a non-significant interaction is removed from the model before interpreting its results. In our case, only the main effects of age and sex would remain in the model.

This leads us to model selection, or which variables should be included in our final model. This is discussed in the next section.

Model selection

In practice, we often have several models, corresponding to the different combinations of independent variables and their interactions. Finding the best model is not easy.

In general, the best practice is to obtain a final model that is as parsimonious as possible, that is, with as few parameters as possible. A parsimonious model is easier to interpret and generalize, and also more powerful from a statistical point of view. On the other hand, it should not be too simple so that it still captures the variations or patterns in the data. In general, while more variables are often better than one, too many is often worse than a few.

The two most common approaches to obtain a final model are the following:

1. Adjust the model by removing the main effects and their interactions which are not significant with respect to their \(p\)-values, obtained by testing the nullity of the corresponding coefficients using a statistical test such as the likelihood ratio test or the Wald test. If there are several main effects or interactions which are not significant, interactions must be removed before removing any main effect. Moreover, it is recommended to remove interactions and independent variables one by one (starting with the one with the highest \(p\)-value), as removing a variable or an interaction may make another variable or interaction that was initially non-significant significant.
2. Adjust the model by using AIC (Akaike Information Criterion) or BIC (Bayesian Information Criterion). These procedures allow to select the best model (according to AIC or BIC) by finding an equilibrium between simplicity and complexity. These selection processes usually lead to a final model with as few parameters as possible, but which captures as much information in the data as possible. Note that only models with the same dependent variable can be compared using AIC or BIC. Models with different dependent variables cannot be compared using these criteria.

The first method requires that:

• the underlying assumptions are valid,
• the sample size is sufficiently large, and
• the models are nested (i.e., the complete model includes at least all the variables included in the reduced model).

Moreover, the second method can be used with widely used criteria when selecting variables, and more importantly, in a completely autonomous way in R.

For this reason, the second method is recommended and more often used in practice.

This second method, referred as the stepwise selection, is divided into 3 types:

1. backward selection: we start from the most complete model (containing all independent variables and usually also their interactions), and the interactions/main effects are deleted at each step until the model cannot be improved,
2. forward selection: we start from the most basic model containing only the intercept, and the independent variables/interactions are added at each step until the model cannot be improved, or
3. mixed selection: we apply both the backward and forward selection to determine the best model according to the desired criterion.

We show how to select the best model according to AIC using the mixed stepwise selection, illustrated with all variables present in the data frame as independent variables and all possible second order interactions:

# save initial model
m5 <- glm(heart_disease ~ (age + sex + chest_pain + max_heartrate)^2,
  data = dat,
  family = "binomial"
)

# select best model according to AIC using mixed selection
m5_final <- step(m5,
  direction = "both", # both = mixed selection
  trace = FALSE # do not display intermediate steps
)

# display results of final model
summary(m5_final)

## 
## Call:
## glm(formula = heart_disease ~ age + sex + chest_pain + max_heartrate + 
##     age:max_heartrate, family = "binomial", data = dat)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                19.4386591  8.1904201   2.373 0.017628 *  
## age                        -0.3050017  0.1414986  -2.156 0.031122 *  
## sexmale                     1.7055353  0.3507149   4.863 1.16e-06 ***
## chest_painatypical angina  -0.0086463  0.6547573  -0.013 0.989464    
## chest_painnon-anginal pain -0.0590333  0.5844000  -0.101 0.919538    
## chest_painasymptomatic      1.9724490  0.5649516   3.491 0.000481 ***
## max_heartrate              -0.1605658  0.0540933  -2.968 0.002994 ** 
## age:max_heartrate           0.0023314  0.0009433   2.471 0.013457 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 409.95  on 296  degrees of freedom
## Residual deviance: 268.60  on 289  degrees of freedom
## AIC: 284.6
## 
## Number of Fisher Scoring iterations: 5

According to the AIC (which is the default criterion when using the step() function), the best model is the one including:

• age,
• sex,
• chest_pain,
• max_heartrate, and
• the interaction between age and max_heartrate.

For your information, you can also easily compare models manually using AIC or the pseudo-\(R^2\) with the tab_model() function, also available in the {sjPlot} R package:

tab_model(m3, m4, m5_final,
  show.ci = FALSE, # remove CI
  show.aic = TRUE, # display AIC
  p.style = "numeric_stars" # display p-values and stars
)

Pseudo-\(R^2\) is a generalization of the coefficient of determination \(R^2\) often used in linear regression to judge the quality of a model. Like the \(R^2\) in linear regression, the pseudo-\(R^2\) varies from 0 to 1, and can be interpreted as the percentage of the null deviance explained by the independent variable(s). The higher the pseudo-\(R^2\) and the lower the AIC, the better the model.

Note that there are several pseudo-\(R^2\), such as:

• Likelihood ratio \(R^2_{L}\),
• Cox and Snell \(R^2_{CS}\),
• Nagelkerke \(R^2_{N}\),
• McFadden \(R^2_{McF}\), and
• Tjur \(R^2_{T}\).

The tab_model() function gives the Tjur \(R^2_{T}\) by default.

Based on the AIC and the Tjur \(R^2_{T}\), the last model is considered as the best one among the 3 considered.

Note that, even though a model is deemed the best one among the ones you have considered (based on one or several criteria), it does not necessarily mean that it fits the data well. There are several methods to check the quality of a model and to check if it is appropriate for the data at hand. This is the topic of the next section.

Quality of a model

Usually, the goal of building a model is to be able to predict, as precisely as possible, the response variable for new data.

In the next sections, we present some measures to judge the quality of a model, starting with the easiest and most intuitive one, followed by two widely used in the medical domain, and finally two other metrics common in the field of machine learning.

# create a vector of predicted probabilities
preds <- predict(m5_final,
  newdata = select(dat, -heart_disease), # remove real outcomes
  type = "response"
)

# if probability < threshold, patient is considered not to have the disease
preds_outcome <- ifelse(preds < 0.5,
  0,
  1
)

# transform predictions into factor and set labels
preds_outcome <- factor(preds_outcome,
  levels = c(0, 1),
  labels = c("no disease", "disease")
)

# compare observed vs. predicted outcome
tab <- table(dat$heart_disease, preds_outcome,
  dnn = c("observed", "predicted")
)

# print results
tab

##             predicted
## observed     no disease disease
##   no disease        132      28
##   disease            33     104

accuracy <- sum(diag(tab)) / sum(tab)
accuracy

## [1] 0.7946128

# sensitivity
sensitivity <- tab[2, 2] / (tab[2, 2] + tab[2, 1])
sensitivity

## [1] 0.7591241

# specificity
specificity <- tab[1, 1] / (tab[1, 1] + tab[1, 2])
specificity

## [1] 0.825

# load package
library(pROC)

# save roc object
res <- roc(heart_disease ~ fitted(m5_final),
  data = dat
)

# plot ROC curve
ggroc(res, legacy.axes = TRUE)

# print AUC
res$auc

## Area under the curve: 0.87

# plot ROC curve with AUC in title
ggroc(res, legacy.axes = TRUE) +
  labs(title = paste0("AUC = ", round(res$auc, 2)))

Reporting results

As we have seen before, odds ratios are useful when reporting results of binary logistic regressions.

Computing these odds ratios together with the confidence intervals is not particularly difficult. However, presenting them in a table for a publication or a report can quickly become time consuming, in particular if you have many models and many independent variables.

Luckily, there are two packages which saved me a lot of time and which I use almost every time I need to report results of a logistic regression.

The first package, called {gtsummary} is useful to report results of one regression at a time. The second one is the {finalfit} package.¹ This packages is more appropriate if you need to report results of several regressions at a time.

# load package
library(gtsummary)

# print table of results
tbl_regression(m5_final, exponentiate = TRUE)

# load packages
library(tidyverse)
library(gt)
library(finalfit)

# set dependent and independent variables
dependent <- "heart_disease"
independent <- c("age", "sex", "chest_pain", "max_heartrate")

# save results of univariable logistic regressions
glmuni <- dat |>
  glmuni(dependent, independent) |>
  fit2df(
    explanatory_name = "Variables",
    estimate_name = "Crude OR",
    estimate_suffix = " (95% CI)"
  )

# print results
glmuni |>
  gt()

# save results of full model
glmmulti_full <- dat |>
  glmmulti(dependent, independent) |>
  fit2df(
    explanatory_name = "Variables",
    estimate_name = "Adjusted OR - full model",
  )

# print results
glmmulti_full |>
  gt()

# select the variables to be included in the final model
independent_final <- c("age", "sex", "chest_pain")

# save results of final model
glmmulti_final <- dat |>
  glmmulti(dependent, independent_final) |>
  fit2df(
    explanatory_name = "Variables",
    estimate_name = "Adjusted OR - final model",
    estimate_suffix = " (95% CI)"
  )

# print results
glmmulti_final |>
  gt()

# save descriptive statistics
summary <- dat |>
  summary_factorlist(dependent, independent, fit_id = TRUE)

# save results of regressions
output <- summary |>
  finalfit_merge(glmuni) |>
  finalfit_merge(glmmulti_full) |>
  finalfit_merge(glmmulti_final)

# print all results
output |>
  dplyr::select(-fit_id, -index) |>
  dplyr::rename(
    Variables = label,
    " " = levels
  ) |>
  gt()

# load packages
library(tidyverse)
library(gt)
library(finalfit)

# set variables
dependent <- "heart_disease"
independent <- c("age", "sex", "chest_pain", "max_heartrate")
independent_final <- c("age", "sex", "chest_pain")

# save descriptive statistics
summary <- dat |>
  summary_factorlist(dependent, independent, fit_id = TRUE)

# save results of univariable logistic regressions
glmuni <- dat |>
  glmuni(dependent, independent) |>
  fit2df(
    explanatory_name = "Variables",
    estimate_name = "Crude OR",
    estimate_suffix = " (95% CI)"
  )

# save results of full model
glmmulti_full <- dat |>
  glmmulti(dependent, independent) |>
  fit2df(
    explanatory_name = "Variables",
    estimate_name = "Adjusted OR - full model",
  )

# save results of final model
glmmulti_final <- dat |>
  glmmulti(dependent, independent_final) |>
  fit2df(
    explanatory_name = "Variables",
    estimate_name = "Adjusted OR - final model",
    estimate_suffix = " (95% CI)"
  )

# save merged results
output <- summary |>
  finalfit_merge(glmuni) |>
  finalfit_merge(glmmulti_full) |>
  finalfit_merge(glmmulti_final)

# print all results
output |>
  dplyr::select(-fit_id, -index) |>
  dplyr::rename(
    Variables = label,
    " " = levels
  ) |>
  gt()

dat |> or_plot(dependent, independent,
  table_text_size = 3.5 # reduce text size
)

Conditions of application

For results to be valid and interpretable, a binary logistic regression requires:

1. the dependent variable to be binary,
2. independence of the observations: no repeated measurements or matched data, otherwise generalize linear mixed effect models (GLMM) should be used,
3. linearity of continuous independent variables and the log-odds outcome: take age and heart disease as an example. If heart disease is more frequent or less frequent as age rises, the model will work well. However, if children and the elderly are at high risk of having a heart disease, but those in middle years are not, then the relationship is not linear, or not monotonic, meaning that the response does not only go in one direction,
4. a sufficiently large sample size (for confidence intervals and hypothesis tests to be valid), and
5. no multicollinearity: independent variables should not be highly correlated with each other, otherwise coefficients and OR can become unstable.

Here is how to verify each of them:

1. this is obvious; check if the dependent variable has indeed only two levels,
2. this is often not tested formally, but verified through the design of the experiment,
3. quantitative independent variables should have a linear relationship between their log-odds and their observed values. A visual check is sufficient, see below with age, maximum heart rate achieved and model m3 as example:

# linearity to the log-odds?
dat |>
  dplyr::select(age, max_heartrate) |>
  mutate(log_odds = predict(m3)) |>
  pivot_longer(-log_odds) |>
  ggplot(aes(log_odds, value)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~name)

4. in practice, it is recommended to have at least 10 times as many events as parameters in the model, and
5. the variance inflation factors (VIF) is a well known measure of multicollinearity. It should be below 10 or 5, depending on the field of research. VIF can be computed with the vif() function, available in the {car} package:

# load package
library(car)

# compute VIF for model m3
vif(m3)

##                   GVIF Df GVIF^(1/(2*Df))
## age           1.205246  1        1.097837
## sex           1.155071  1        1.074742
## chest_pain    1.113010  3        1.018005
## max_heartrate 1.143125  1        1.069170