Chapter 9 Multinomial Logistic Regression

9.1 Objectives

At the end of this chapter, readers should be able:

  • to understand the concept of logistic regression model to analyze data with polychotomous (multinomial) outcome
  • to estimate parameters of interest in a logistic regression model from data with polychotomous (multinomial) outcome
  • to make inference based on a logistic regression model from data with polychotomous (multinomial) outcome
  • to predict the outcome based on a logistic regression model from data with polychotomous (multinomial) outcome
  • to perform model checking of logistic regression model from data with with polychotomous (multinomial) outcomes

9.2 Introduction

Some data come with multinomial outcomes in which case, the outcome variable is a nominal or polychotomous variable with more than two levels. In multinomial outcome data, the outcome has no natural ordering. if it has, then it is best treated as a ordinal outcome data. For these data, we can modify the binary logistic model to make estimation and inference (D. Kleinbaum 2010 ; Hosmer, Lemeshow, and Sturdivant 2013).

Variables with more than two levels are known as either:

  1. multinomial data
  2. polychotomous data
  3. polytomous data

If we employ logistic regression to such data, then the analysis is known as polytomous or multinomial logistic regression. Again, the polychotomous outcome does not have any natural order. When the categories of the outcome variable do have a natural order, ordinal logistic regression is preferred.

9.3 Examples of multinomial outcome variables

Examples of data with polychotomous outcome variables (with more than two levels) include:

  • disease symptoms that have been classified by subjects as being absent, mild, moderate, or severe,
  • tumour invasiveness classified as in situ, locally invasive, ormetastatic, or
  • patient preferred treatment regimen, selected from among three or more options for example oral medication only, oral medication plus injection medication or injection only.

A numerical outcome can be categorized based on different cut-off points. The newly created categorical variable is now either treated as a nominal or polychotomous or multinomial outcome or as a ordinal outcome. That justifies the use of multinomial logistic regression.

9.4 Models for multinomial outcome data

With a multinomial outcome data, an extension of logistic regression know as multinomial logit or multinomial logistic regression can be performed. In multinomial logistic regression, one of the categories of the outcome variable is designated as the reference category and each of the other levels is compared with this reference. The choice of reference category can be arbitrary and is at the discretion of the researcher. Most software set the first category as the reference (also known as the baseline category) by default.

Other models that can analyze data with polychotomous outcome include:

  1. Stereotype logistic regression - each independent variable has one value for each individual
  2. Alternative-specific variables

9.5 Estimation for Multinomial logit model

Remember, interpreting and assessing the significance of the estimated coefficients are the main objectives in regression analysis. in multinomial logistic regression, we would like to model the relationship between covariates with the outcome variable that has more than two categories but without ordering or ranking.

The actual values taken by the dependent variable are irrelevant. In Stata, the exponentiated beta \(\exp^{\beta}\) will generate the so-called the relative-risk ratio. The dependent variable again, is a discrete variable and the model is estimated using maximum likelihood.

In multinomial logistic regression for example in data with three categories of the outcome, the sum of the probabilities for the three outcome categories must be equal to 1 (the total probability). The comparison is done between two categories each time. Because of that, each comparison considers only two probabilities, and the probabilities in the ratio do not sum to 1. Thus, the two odds-like expressions in multinomial logistic regression are not true odds.

Multinomial logistic regression can be thought as of simultenously fitting binary logits for all comparisons among the alternatives.

9.5.1 Log Odds and Odds Ratios

In the special case where the covariate is binary, coded 0 or 1, we simplify the notation to \(OR_j = OR_j(1,0)\). Let use an example where data have 3 categories of outcome; 0,1 and 2.

Let’s say we have a dataset with the outcome variable, \(Y\), and is coded as \(0, 1,\) or \(2\). In practice one should check that the software package that is going to be used allows a \(0\) code as the smallest category. We have used packages that require that the codes begin with \(1\).

SO the logit functions (log odds) when the outcome is a D variable with (0,1 and 2 values) are as below

  1. the log odds for comparison between 0 and 1 is

\[g_1(x) = ln\frac{P(D=1|X_1)} {P(D=0|X_1)}=\alpha_1 + (\beta_{11}X_1)\]

  1. and, the log odds for comparison between 0 and 2 is

\[g_2(x) = ln\frac{P(D=2|X_1)} {P(D=0|X_1)}=\alpha_2 + (\beta_{21}X_1)\]

If for example, we assume that the outcome labelled with \(Y=0\) is the reference outcome. The subscript on the odds ratio indicates which outcome is being compared to the reference category outcome. The odds ratio of outcome \(Y = j\) versus \(Y = 0\) for the covariate values of \(x = a\) versus \(x = b\) is:

\[OR_{j}(a,b)= \frac{Pr(Y=j|x=a)/(Pr(Y=0|x=a)} {Pr(Y=j|x=b)/Pr(Y=0|x=b)}\]

Each odds ratio is calculated in a manner similar to that used in standard logistic regression. That is:

\[OR_1(X=1,X=0)= \frac{Pr(D=1|X=1)/(Pr(D=0|X=1)} {Pr(Y=1|X=0)/Pr(D=0|X=0)}=\exp^{\beta_{11}}\]

\[OR_2(X=2,X=0)= \frac{Pr(D=2|X=1)/(Pr(D=0|X=1)} {Pr(D=2|X=0)/Pr(D=0|X=0)}=\exp^{\beta_{21}}\]

9.5.2 Conditional probabilities

The conditional probabilities for the multinomial logit model are:

\[Pr(D = 0 | x) = \frac{1}{1+e^{g_1(x)} + e^{g_2(x)}}\]

\[Pr(D = 1 | x) = \frac{e^{g_1(x)}}{1+e^{g_1(x)} + e^{g_2(x)}}\]

\[Pr(D = 2 | x) = \frac{e^{g_2(x)}}{1+e^{g_1(x)} + e^{g_2(x)}}\]

9.6 Prepare environment

Start brand new analysis with a new project in RStudio. To do so,

  • Go to RStudio Menu
  • Click File
  • Click New Project
  • Choose either New Directory or Existing Directory (to point to existing folder usually with the dataset)

9.6.1 Load libraries

We will load

  • here package that enables easy file referencing
  • tidyverse for data wrangling and plotting
  • haven to read data in various statistical formats
  • gtsummary to produce statistical tables
  • VGAM package where we will be using the vglm() function to perform multinomial logistic regression
  • kableExtra to produce nice tables for the results
library(here)
## here() starts at D:/Data_Analysis_CRC_multivar_data_analysis_codes
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(haven)
library(gtsummary)
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

9.6.2 Dataset

The datasets contains all variables of our interest. For the purpose of this assignment, we want to explore association of hypertension status, weight and total cholesterol with the result of screening FBS. The variables in the datasets as follow:

  1. fbs : Fasting Blood Sugar (mmol/L). The data ranges from 2.51 mmol/L to 28.01 mmol/L.
  2. totchol : Total Cholesterol (mmol/L). The data ranges from 0.18 mmol/L to 23.14 mmol/L.
  3. hpt : Hypertension Status. Coded as Yes and No.
  4. weight : Body weight measures in kilogram. The data ranges between 30kg to 187.8kg.

9.6.3 Read data

We will use haven::read_dta() function to read stata .dta data into the memory. Remember, we can also use foreign::read.dta() function to read stata .dta data. It is your choice.

We use here() to indicate that the folder data contains the dataset. And then we specify metabolic_syndrome.dta as the dataset to be read.

ms <- read_dta(here('data','metabolic_syndrome.dta'))

We can then quickly glance at the number of observations and the type of variables in the ms dataset.

glimpse(ms)
## Rows: 4,340
## Columns: 21
## $ codesub  <chr> "R-S615112", "MAA615089", "M-M616372", "MFM615361", "R-A61578…
## $ age      <dbl> 70, 20, 29, 25, 37, 43, 26, 28, 48, 20, 56, 55, 26, 39, 18, 2…
## $ hpt      <chr> "yes", "no", "no", "no", "no", "no", "no", "no", "no", "no", …
## $ smoking  <chr> "never smoked", "still smoking", "never smoked", "still smoki…
## $ dmdx     <chr> "no", "no", "no", "no", "no", "no", "no", "no", "no", "no", "…
## $ height   <dbl> 1.54, 1.74, 1.54, 1.60, 1.44, 1.46, 1.47, 1.61, 1.68, 1.55, 1…
## $ weight   <dbl> 40.0, 54.6, 37.0, 48.4, 44.5, 45.5, 48.0, 40.0, 48.4, 41.0, 5…
## $ waist    <dbl> 76.0, 83.0, 83.0, 83.5, 85.0, 90.0, 91.0, 80.0, 82.0, 79.0, 9…
## $ hip      <dbl> 61, 62, 63, 64, 64, 64, 64, 64, 65, 66, 67, 67, 67, 67, 67, 6…
## $ msbpr    <dbl> 135.0, 105.0, 91.0, 117.0, 102.0, 124.0, 120.0, 85.0, 112.0, …
## $ mdbpr    <dbl> 80.0, 58.0, 60.0, 68.5, 78.0, 65.5, 77.0, 60.0, 74.0, 52.0, 1…
## $ hba1c    <dbl> 5.2, 5.3, 4.8, 4.8, 5.1, 5.1, 4.8, 4.9, 5.6, 4.2, 5.1, 5.3, 4…
## $ fbs      <dbl> 3.99, 4.26, 4.94, 4.60, 4.60, 4.42, 3.82, 4.40, 4.80, 3.68, 6…
## $ mogtt1h  <dbl> 7.06, 8.63, 6.26, 4.31, 9.49, 6.29, NA, 6.43, 9.23, 6.70, 8.7…
## $ mogtt2h  <dbl> 3.22, 6.49, 5.15, 3.85, 7.71, 5.65, 5.88, 4.89, 4.29, 2.59, 7…
## $ totchol  <dbl> 5.43, 5.13, 5.55, 4.01, 5.21, 6.19, 4.33, 5.84, 6.14, 6.02, 6…
## $ ftrigliz <dbl> 1.06, 1.17, 0.72, 1.12, 0.78, 1.11, 0.73, 0.79, 1.63, 0.81, 1…
## $ hdl      <dbl> 1.65, 1.59, 2.24, 1.21, 1.43, 2.18, 0.98, 1.81, 1.63, 1.47, 1…
## $ ldl      <dbl> 2.69, 2.79, 2.55, 1.83, 2.40, 2.93, 1.82, 3.43, 3.71, 2.77, 3…
## $ gender   <chr> "female", "male", "female", "male", "female", "female", "fema…
## $ crural   <chr> "rural", "rural", "rural", "rural", "rural", "rural", "rural"…

9.6.4 Data Wrangling

We will

  • select only variable fbs, totchol, hpt and weight
  • convert all character variables to factor variables
  • then take a glance of the updated ms dataset again
ms <- ms %>% 
  select(fbs, totchol, hpt , weight) %>% 
  mutate(across(where(is.labelled), as_factor))
glimpse(ms)
## Rows: 4,340
## Columns: 4
## $ fbs     <dbl> 3.99, 4.26, 4.94, 4.60, 4.60, 4.42, 3.82, 4.40, 4.80, 3.68, 6.…
## $ totchol <dbl> 5.43, 5.13, 5.55, 4.01, 5.21, 6.19, 4.33, 5.84, 6.14, 6.02, 6.…
## $ hpt     <chr> "yes", "no", "no", "no", "no", "no", "no", "no", "no", "no", "…
## $ weight  <dbl> 40.0, 54.6, 37.0, 48.4, 44.5, 45.5, 48.0, 40.0, 48.4, 41.0, 52…

9.6.5 Create new categorical variable from fbs

Let us create a categorical (also known as a factor variable) based on this classification:

  1. Normal if fbs is less than 6.1 mmol/L
  2. Impaired Fasting Glucose (IFG) if fbs is between 6.1 mmol/L to 6.9 mmol/L
  3. Diabetis Mellitus (DM) if fbs is 7.00 mmol/L or higher
ms <- ms %>%
  mutate(cat_fbs = cut(fbs, 
                       breaks = c(2.50 , 6.10 , 6.90 , 28.01 ),
                       labels = c("Normal","IFG", "DM")))
ms %>% 
  count(cat_fbs)
## # A tibble: 4 × 2
##   cat_fbs     n
##   <fct>   <int>
## 1 Normal   3130
## 2 IFG       364
## 3 DM        595
## 4 <NA>      251

We notice that there were 250 data has no values recorded as NA. Thus, we decide to remove observations when there are missing values for variable cat_fbs.

ms <- ms %>%
  filter(!is.na(cat_fbs)) 
ms %>%
  count(cat_fbs)
## # A tibble: 3 × 2
##   cat_fbs     n
##   <fct>   <int>
## 1 Normal   3130
## 2 IFG       364
## 3 DM        595

9.6.6 Exploratory data analysis

Next, is to return the table of summary statistics

ms %>%
  tbl_summary(by = cat_fbs,
              statistic = list(all_continuous() ~ "{mean} ({sd})", 
                              all_categorical() ~ "{n} ({p}%)"),
              type = list(where(is.logical) ~ "categorical"),
              label = list(fbs ~ "Fasting Blood Sugar (mmol/L)", 
                           totchol ~ "Total Cholesterol (mmol/L)", hpt~"Hypertension", weight ~ "Weight (kg)"),
              missing_text = "Missing") %>% 
  modify_caption("**Table 1. Survey Participant Characteristic**")  %>%
  modify_header(label ~ "**Variable**") %>%
  modify_spanning_header(c("stat_1", "stat_2", "stat_3") ~ "**Glycemic Control Status**") %>%
  modify_footnote(all_stat_cols() ~ "Mean (SD) or Frequency (%)") %>%
  bold_labels() %>%
  as_gt()
TABLE 6.3: Table 1. Survey Participant Characteristic
Variable Glycemic Control Status
Normal
N = 3,1301		 IFG
N = 3641		 DM
N = 5951
Fasting Blood Sugar (mmol/L) 4.76 (0.82) 6.44 (0.22) 10.41 (3.63)
Total Cholesterol (mmol/L) 5.69 (1.24) 6.08 (1.37) 6.16 (1.31)
    Missing 4 0 1
Hypertension 281 (9.0%) 75 (21%) 126 (21%)
Weight (kg) 63 (14) 68 (14) 68 (15)
    Missing 1 0 0
1 Mean (SD) or Frequency (%)

9.6.7 Confirm the order of cat_fbs

levels(ms$cat_fbs)
## [1] "Normal" "IFG"    "DM"

However, we would like the DM as the smallest category. To do that we will use the fct_relevel() function.

ms <- ms %>% 
  mutate(cat_fbs = fct_relevel(cat_fbs, 
                               c("DM", 'IFG', 'Normal')))
levels(ms$cat_fbs)
## [1] "DM"     "IFG"    "Normal"

9.7 Estimation

Our intention to investigate the relationship between totchol, hpt and weight with the outcome variables cat_fbs. Thus, we will perform multinomial logistic regression model to estimate the relation for 2 models:

  • Model 1: DM vs Normal
  • Model 2: IFG vs Normal

In both models, the reference group is Normal

9.7.1 Single Independent Variable

For independent variable totchol

log_chol <- vglm(cat_fbs ~ totchol, 
                 multinomial, data = ms)
summary(log_chol)
## 
## Call:
## vglm(formula = cat_fbs ~ totchol, family = multinomial, data = ms)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -3.33584    0.21306 -15.656  < 2e-16 ***
## (Intercept):2 -3.59935    0.25764 -13.970  < 2e-16 ***
## totchol:1      0.28357    0.03446   8.229  < 2e-16 ***
## totchol:2      0.24661    0.04176   5.905 3.53e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
## 
## Residual deviance: 5634.468 on 8164 degrees of freedom
## 
## Log-likelihood: -2817.234 on 8164 degrees of freedom
## 
## Number of Fisher scoring iterations: 4 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):2'
## 
## 
## Reference group is level  3  of the response

This is the model where hpt is the independent variable

log_hpt <- vglm(cat_fbs ~ hpt, 
                 multinomial, data = ms)
summary(log_hpt)
## 
## Call:
## vglm(formula = cat_fbs ~ hpt, family = multinomial, data = ms)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -1.80412    0.04983 -36.205  < 2e-16 ***
## (Intercept):2 -2.28830    0.06173 -37.067  < 2e-16 ***
## hptyes:1       1.00205    0.11823   8.475  < 2e-16 ***
## hptyes:2       0.96743    0.14389   6.724 1.77e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
## 
## Residual deviance: 5637.16 on 8174 degrees of freedom
## 
## Log-likelihood: -2818.58 on 8174 degrees of freedom
## 
## Number of Fisher scoring iterations: 4 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Reference group is level  3  of the response

And lastly, the independent variable is weight

log_wt <- vglm(cat_fbs ~ weight, 
                 multinomial, data = ms)
summary(log_wt)
## 
## Call:
## vglm(formula = cat_fbs ~ weight, family = multinomial, data = ms)
## 
## Coefficients: 
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -3.214303   0.204303 -15.733  < 2e-16 ***
## (Intercept):2 -3.672632   0.249121 -14.742  < 2e-16 ***
## weight:1       0.023860   0.002988   7.984 1.42e-15 ***
## weight:2       0.023372   0.003627   6.444 1.16e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
## 
## Residual deviance: 5638.998 on 8172 degrees of freedom
## 
## Log-likelihood: -2819.499 on 8172 degrees of freedom
## 
## Number of Fisher scoring iterations: 4 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):2'
## 
## 
## Reference group is level  3  of the response

9.7.2 Multiple Independent Variables

We feel that totchol, hpt and weight are all important independent variables. Hence, we want to fit a model with the three independent variables as the covariates.

mlog <- vglm(cat_fbs ~ totchol + hpt + weight, 
             multinomial, data = ms)
summary(mlog)
## 
## Call:
## vglm(formula = cat_fbs ~ totchol + hpt + weight, family = multinomial, 
##     data = ms)
## 
## Coefficients: 
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -4.907874   0.303768 -16.157  < 2e-16 ***
## (Intercept):2 -5.112990   0.366325 -13.958  < 2e-16 ***
## totchol:1      0.277217   0.035055   7.908 2.62e-15 ***
## totchol:2      0.239413   0.042391   5.648 1.63e-08 ***
## hptyes:1       0.899987   0.120451   7.472 7.91e-14 ***
## hptyes:2       0.867136   0.145569   5.957 2.57e-09 ***
## weight:1       0.022678   0.003074   7.378 1.61e-13 ***
## weight:2       0.021994   0.003710   5.928 3.07e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
## 
## Residual deviance: 5476.397 on 8158 degrees of freedom
## 
## Log-likelihood: -2738.199 on 8158 degrees of freedom
## 
## Number of Fisher scoring iterations: 5 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):1', '(Intercept):2'
## 
## 
## Reference group is level  3  of the response

9.7.3 Model with interaction term between independent variables

Then, we hypothesize that there could be a significant interaction between totchol and weight. And to test the hypothesis, we extend the multivariable logistic regression model by adding an interaction term. This interaction term is a product between variable weight and totchol.

mlogi <- vglm(cat_fbs ~ totchol + hpt + weight + totchol*weight, 
              multinomial, data = ms)
summary(mlogi)
## 
## Call:
## vglm(formula = cat_fbs ~ totchol + hpt + weight + totchol * weight, 
##     family = multinomial, data = ms)
## 
## Coefficients: 
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1    -5.1253325  0.9671480  -5.299 1.16e-07 ***
## (Intercept):2    -7.0506874  1.1195704  -6.298 3.02e-10 ***
## totchol:1         0.3155740  0.1605324   1.966  0.04932 *  
## totchol:2         0.5747241  0.1878902   3.059  0.00222 ** 
## hptyes:1          0.8984459  0.1204877   7.457 8.87e-14 ***
## hptyes:2          0.8643701  0.1455478   5.939 2.87e-09 ***
## weight:1          0.0260899  0.0142834   1.827  0.06776 .  
## weight:2          0.0514503  0.0165054   3.117  0.00183 ** 
## totchol:weight:1 -0.0006015  0.0023804  -0.253  0.80052    
## totchol:weight:2 -0.0051091  0.0028023  -1.823  0.06828 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
## 
## Residual deviance: 5473.004 on 8156 degrees of freedom
## 
## Log-likelihood: -2736.502 on 8156 degrees of freedom
## 
## Number of Fisher scoring iterations: 5 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):1', '(Intercept):2'
## 
## 
## Reference group is level  3  of the response

The interaction term in our model showed p-values above 0.05 (p-values of 0.80 and 0.07, respectively). As the p-value is bigger than the level of significance at \(5\%\) and the value of regression parameters for the interaction terms are likely not clinically meaningful, we have decided not to use the model with an interaction term.

9.8 Inferences

For the inference, we will

  • calculate the \(95\%\) CI (interval estimates)
  • calculate the p-values (hypothesis testing)

There is no facility inside thebroom::tidy() function to generate confidence intervals for object with class vglm. Because of that we will use the coef(), confint() and cind(0 functions to produce a rather nice table of inferences.

We are going to follow these steps:

  • set the number of digits equal to 2 to limit the decimal numbers
  • return the regression coefficents for all \(\hat\beta\) as an object named b_fitmlog2
  • return the the confidence intervals for all \(\hat\beta\) as an object named ci_fitmlog2
  • combine the \(\hat\beta\) and the corresponding \(95\%\) CIs
b_mlog <- coef(mlog)
ci_mlog <- confint(mlog) 
b_ci_mlog <- data.frame(b_mlog,ci_mlog) %>%
  rename("log odds" = b_mlog, "Lower CI" = X2.5.., "Upper CI" = X97.5..)
b_ci_mlog %>% 
  kbl(digits = 2, booktabs = T, caption = "Log odds from multinomial logistic regression") %>%
  kable_styling(position = "center")
TABLE 9.1: Log odds from multinomial logistic regression
log odds Lower CI Upper CI
(Intercept):1 -4.91 -5.50 -4.31
(Intercept):2 -5.11 -5.83 -4.40
totchol:1 0.28 0.21 0.35
totchol:2 0.24 0.16 0.32
hptyes:1 0.90 0.66 1.14
hptyes:2 0.87 0.58 1.15
weight:1 0.02 0.02 0.03
weight:2 0.02 0.01 0.03

Afterwards, we will exponentiate the coefficients to obtain the relative-risk ratio. We then combine the results to the previous table. Finally, we will name the columns of the object tab_fitmlog2.

rrr_mlog <- exp(b_ci_mlog)
tab_mlog <- cbind(b_ci_mlog, rrr_mlog)
colnames(tab_mlog) <- c('b', 'lower b', 'upper b',
                        'RRR', 'lower RRR', 'upper RRR')
tab_mlog %>%
  kbl(digits = 2, booktabs = T, caption = "Log odds and RRR from multinomial logistic regression") %>%
  kable_styling(position = "center")
TABLE 9.2: Log odds and RRR from multinomial logistic regression
b lower b upper b RRR lower RRR upper RRR
(Intercept):1 -4.91 -5.50 -4.31 0.01 0.00 0.01
(Intercept):2 -5.11 -5.83 -4.40 0.01 0.00 0.01
totchol:1 0.28 0.21 0.35 1.32 1.23 1.41
totchol:2 0.24 0.16 0.32 1.27 1.17 1.38
hptyes:1 0.90 0.66 1.14 2.46 1.94 3.11
hptyes:2 0.87 0.58 1.15 2.38 1.79 3.17
weight:1 0.02 0.02 0.03 1.02 1.02 1.03
weight:2 0.02 0.01 0.03 1.02 1.01 1.03

9.9 Interpretation

The result from our multivariable logistic regression models can be interpreted as below:

  • Every increment 1 mmol/L of totchol (Total cholesterol) when controlling for hypertension status and weight; a) the odds of being in DM group (in comparison to Normal) change by 0.277 (Adjusted RRR = 1.32, \(95\%\) CI : 1.232, 1.413, p-value <0.001) and b) the odds of being in IFG group (in comparison to being in Normal) change by 0.239 (Adjusted RRR = 1.27, \(95\%\) CI : 1.169,1.381, p-value <0.001).
  • Every increment 1 kg of weight when controlling for hypertension status and total cholesterol; a) the odds of being in DM group (in comparison to being in Normal) increase by 0.023 (Adjusted RRR = 1.02, \(95\%\) CI : 1.017,1.029, p-value <0.001), and b) the odds of being in IFG group (in comparison being in Normal) increase by 0.022 (Adjusted RRR = 1.02, \(95\%\) CI : 1.015,1.030, p-value <0.001).
  • In the population with hypertension (as compared with participant without hypertension) when controlling for weight and total cholesterol; a) their odds of being in DM group (in comparison to being in Normal) change by 0.900 (Adjusted RRR = 2.5, \(95\%\) CI : 1.942, 3.114, p-value <0.001). When repeated research on this population, the odds ranging between reduced the odds by 0.481 to 0.677, and b) and their odds of being in IFG group (in comparison to being Normal) change by 0.867 (Adjusted RRR = 2.38, 95% CI : 1.789, 3.166, p-value <0.001).

9.10 Prediction

We can calculate the predicted probability of each category of outcome by using the predict() function. Below are the result for the top 10 observations.

predict.vgam(mlog, type = 'response') %>% 
  head(10)
##            DM        IFG    Normal
## 1  0.15254733 0.09528434 0.7521683
## 2  0.09000250 0.05817369 0.8518238
## 3  0.07042058 0.04534241 0.8842370
## 4  0.06047226 0.04095046 0.8985773
## 5  0.07525337 0.04882982 0.8759168
## 6  0.09712040 0.06068516 0.8421944
## 7  0.06497274 0.04348092 0.8915463
## 8  0.08025639 0.05100727 0.8687363
## 9  0.10144228 0.06337974 0.8351780
## 10 0.08548801 0.05392687 0.8605851

9.11 Presentation of multinomial regression model

We can make a better table using knitr and kableExtra packages. Or we can save the results as a .csv file so we can edit it using spreadsheets.

tab_mlog %>%
  write_csv("results_multinomial.csv")

9.12 Summary

In this chapter, we extend binary logistic regression analysis to data with a categorical outcome variable but has three or more levels. We describe the multinomial logistic regression model, the log odds and also the conditional probabilities. When then show estimation for the single and multiple multinomial logistic regression including models with an interaction. Following that, we describe how to make inferences and perform predictions. Lastly, we make use of kableExtra package to produce a nice table for our model. We are big fans of two excellent books on logistic regression. The first one is Applied Logistic Regression (Hosmer, Lemeshow, and Sturdivant 2013) and the second one is Logistic Regression: Self Learning Text (D. Kleinbaum 2010). You will acquire strong understanding of the concepts and model building process from the two books. To see what packages and functions relevant to logistic regression, you can refer to A Handbook of Statistical Analyses Using R (Everitt and Hothorn 2017).

References

Everitt, Brian S., and Torsten Hothorn. 2017. HSAUR: A Handbook of Statistical Analyses Using r (1st Edition). https://CRAN.R-project.org/package=HSAUR.
Hosmer, D. W., S. Lemeshow, and R. X. Sturdivant. 2013. Applied Logistic Regression. Wiley Series in Probability and Statistics. Wiley. https://books.google.com.my/books?id=bRoxQBIZRd4C.
Kleinbaum, David. 2010. Logistic Regression : A Self-Learning Text. New York: Springer.