Implenting Average Marginal Effects with Bayesian Ordinal Regression

The R package margins provides average marginal effects for any model developed with lm or glm, among others. However, margins does not work with brms. Jack O’Bailey developed two functions for brms here. There are two functions in this text file. The first function gets the average marginal effects for continuous independent variables. The second function get the average marginal effects for categorical variables. The link compares the results of the two functions with the margins command from the package margins. If you run the example, you can see that the results are remarkably similar.

Using the code above, this post attempts to recreate the functionality of the margins command from the margins package with ordinal logistic regression models.

The adapted function is below:

bayes_dydx_ordinal.factor <- function(model, data = NULL, variable, re_formula = NULL){
  
  # Get data from model where data = NULL
  if(is.null(data) == T){
    d <- model$data
  } else {
    d <- data
  }

  
  # Get outcome from model
  resp <- model$formula$resp
  resp_ord <- model$formula$resp
  # Omit outcome from data
  d <-
    d %>%
    dplyr::select(-resp)
  
  d <-
    d %>% 
    group_by_all() %>% 
    count(name = "w") %>% 
    ungroup()
  

  # Get factor levels
  levs <- levels(as.factor(d[[variable]]))
  base <- levs[1L]
  cont <- levs[-1L]
  
  
  # Create empty list for fitted draws
  f <- list()
  
  # For each list add fitted draws
  for (i in seq_along(levs)){
    
    # Fix variable in each list to factor level
    d[[variable]] <- levs[i]
    
    f[[i]] <- 
      d %>% 
      add_fitted_draws(model = model,
                       re_formula = NULL,
                       value = "eff") %>% 
      ##### Need to include ordered within response variable
      group_by_at(vars(".draw", ".category")) %>% 
      #group_by_at(vars(".draw")) %>% 
      summarise(
        # eff = sum(eff),
        # w = sum(w),
        # eff_w_a = sum(eff*w),
        eff_w = sum(eff * w)/sum(w)) %>% 
      #dplyr::select(eff_w) %>% 
      #dplyr::select(eff_w, w, eff, eff_w_a) %>% 
      ungroup() %>% 
      select(-.draw)
    
    # Compute contrast if not base level
    if (i > 1){
      f[[i]]$eff_w <- f[[i]]$eff_w - f[[1]][[2]] #### Change here because also included the order of response
    }
    
    # # Rename column
    names(f[[i]]) <- c(paste0("Category", levs[i]), levs[i])
    #names(f[[i]]) <- levs[i]
    
  }

  # Remove data frame
  d <- NULL
  
  # Create output object
  out <- do.call(cbind, f)

  # Return AMEs
  if (length(cont) == 1){
    
    out <- out %>% tibble() %>% select(1,cont)
    names(out)[2] <- "est"
    
    out %>% 
      mutate(
        var = variable,
        resp = cont
      ) %>% 
      dplyr::select(Category=1,var, resp, est)
    
  } else {
    
    out %>% 
      select(1, cont) %>% 
      tidyr::pivot_longer(cols = -1 ,
                          names_to = "variable", values_to = "value") %>%
      rename(resp = variable, est = value,
             Category = 1)
    
    
  }
  
} ### Close Function

The example uses the nhanes2 dataset, which is from the webuse package.

We next run the ordinal Bayesian model using brms. The code is below:

# Fit Bayesian model
bayes_ordinal <- brm(formula = health ~ female + black + age,
             family=cumulative("logit"),
             prior = c(prior(normal( 0, 1), class = b)),
             data = nhanes2_a,
            warmup = 1000, 
            iter = 2000, 
             chains = 2,
             cores = 2)

AMEByVar = function(x,y) {
  
 ParedEst = bayes_dydx_ordinal.factor(bayes_ordinal , nhanes2_a, variable = x ) %>% 
    group_by(Category, resp) %>% 
    summarise(p50 = median(est)) %>% 
    ungroup() %>% 
    mutate(Var = x)

}
 
BayesAME = purrr::map(as.list(c("black", "female")), AMEByVar) %>% 
  bind_rows(.)

As a first step, we can compare the log odds ratios between the Bayesian model with normally distributed priors and the previous model with flat priors:

As we can see, the coefficients are of similar value.

We next use the adapted function bayes_dydx_ordinal.factor to compare the marginal effects.

Below we compare the marginal effects from the bayes_dydx_ordinal.factor function to the margins function. Importantly, there is no documentation on how margins works with ordinal regressions. There is great documentation here and a great technical paper, but neither covers ordinal regression.

As we can see the marginal effects from the margins package do not match bayes_dydx_ordinal.factor. However, the margins estimates do not vary by outcome level. They are very similar at the first level, where health status is equal to poor compared to the estimates from bayes_dydx_ordinal.factor, but vary considerably from the other outcomes.

 f[[i]] <- 
      d %>% 
      add_fitted_draws(model = model,
                       re_formula = NULL,
                       value = "eff") %>% 
      ##### Need to include ordered within response variable
      group_by_at(vars(".draw", ".category")) %>% 
      #group_by_at(vars(".draw")) %>% 
      summarise(
        # eff = sum(eff),
        # w = sum(w),
        # eff_w_a = sum(eff*w),
        eff_w = sum(eff * w)/sum(w)) %>% 
      #dplyr::select(eff_w) %>% 
      #dplyr::select(eff_w, w, eff, eff_w_a) %>% 
      ungroup() %>% 
      select(-.draw)
    
    # Compute contrast if not base level
    if (i > 1){
      f[[i]]$eff_w <- f[[i]]$eff_w - f[[1]][[2]] #### Change here because also included the order of response
    }