Mediation when using Cumulative Probit and Negative Binomial Models

Modeling
,
1

Hi there,

I am carrying out a mediation analysis in brms that contains an ordinal mediator, modelled with the cumulative(probit) family, and a zero-inflated outcome, modelled with the negative binomial family. For example:

bf_med <- bf(m ~ x) + cumulative(probit)
bf_out <- bf(y ~ m + x) + negbinomial()
m1 <- brm(bf_med + bf_out + set_rescor(FALSE), data=df, cores=4)

Based on this post I understand that the different paths must be standardized to calculate the mediation, direct, and indirect effects.

My questions are:

  1. To calculate the direct effect of the treatment on the outcome, the x → y pathway, would this be standardized in the same way as the mediator on outcome effect, the m → y pathway? As they are both from the second sub-model.

  2. When standardizing the x → y and m → y pathways, how should the negative binomial family be accounted for? I’m guessing exp() should be used somewhere, to inverse the log link function, but should this be used directly on the posterior, after standardization, or is something different required?

The code below is a modified version the simulation/ mediation code from the previous post, but now to simulate y as zero-inflated count data. This is my first attempt at simulation/ reprex, so please let me know if there are any issues.

Thanks,
Josh

library(brms)
library(tidyverse)


set.seed(100)


###Data simulation
n_obs <- 100
x_alpha <- 0
x_sigma <- 1

mx_beta <- 0.5
m_tau <- c(-1, 0, 1)

#regression coefficients
beta0 <- 1
beta1 <- 0.5


# Ordered probit random draws
rop <- function(x, tau){
  op <- vector(mode = 'numeric', length = length(x))
  n_levels <- length(tau) + 1
  probs <- rep(NA, times = n_levels + 1)
  probs[1] <- 0
  probs[n_levels + 1] <- 1
  
  for(n in 1:length(x)){
    probs[2:n_levels] <- pnorm(tau - x[n])
    
    op[n] <- sample(1:n_levels, size = 1, prob = probs[2:(n_levels + 1)] - probs[1:n_levels])
  }
  
  op
}




df <- data.frame(row = 1:n_obs) %>%
  mutate(x = rnorm(n(), x_alpha, x_sigma),
         m = rop(mx_beta * x, m_tau),
         y = rnbinom(n_obs, mu = exp(beta0 + beta1*x + m), size = 0.8))



###brms model

bf_med <- bf(m ~ x) + cumulative(probit)
bf_out <- bf(y ~ m + x) + negbinomial()
m1 <- brm(bf_med + bf_out + set_rescor(FALSE), data= df, cores=4)

summary(m1)


###posterior samples and compute the indirect effect (mediation effect) a*b, as per Andrew Johnson approach

s1 <- as_draws_df(m1) 

sd_x <- sd(df$x)
sd_y <- sd(df$y)
sd_m <- sd(df$m)


#calculate standardized a, the standardized x -> m pathway

a <- s1$b_m_x * sd_x / ( sqrt( s1$b_m_x^2 + sd_x^2 + 1 ) ) 


#calculate standardized b, the standardized m -> y pathway - mediator effect

## Q2. exp() in here somewhere?

b <- s1$b_y_m * sd_m / sd_y


#calculate a*b, the indirect effect (mediation effect)

ab <- a*b


## calculate the standardized direct effect, x -> y pathway 

## Q1. standardize the same way as for the mediator effect, b?

## Q2. exp() in here somewhere?

direct <- s1$b_y_x * sd_x / sd_y


median(a) ## standardized a

median(b) ## standardized b, mediator effect

median(ab) ## standardized indirect effect

median(direct) ## standardized direct effect

median(ab+direct) ## standardized total effect
2

Hi, @joshb, and welcome to the Stan forums. Sorry nobody has responded. Your question is fine, it’s just that we’re overwhelmed with brms questions because very few of our developers use brms. @paul.buerkner will know enough stats and brms to answer, but I’m not sure who else would.

3

I don’t think that the coefficients have to be standardized for the multiplication in mediation effects to work. However, I think that a lot of the math of mediation models relies on linear modeling. This may be extended to cumulative models because they have a linear latent variable, which can be used for this purpose. But I honestly don’t know what math is correct for your case. What you could do it to extend the simulation code in the other post to cover your case and to think about what data generating process makes most sense. Then, we could see if it can be represented in brms and how to extract the relevant parameters of interest, if possible.

4

Hi both, and thank you for your replies.

I’ve edited the simulation code to more closely match my case. x is a five level Likert-type measure (though treated as numeric), m is again a five level Likert-type item and y is a zero inflated count variable. I guess my question is whether it is possible for the effect of x and m on the count variable y to be standardized and directly compared to the effect of x on the latent variable for m, m_hat?

Thanks,
Josh


set.seed(100)


###Data simulation

n_obs <-1000

## x = likert-type data with 5 levels and a bimodal distribution

x <- sample(1:5, n_obs, replace = TRUE, prob = c(0.1, 0.35, 0.2, 0.35, 0.1))

norm_distr <- rnorm(1000, 2.5, 1) ## normal distribution for ordered probit random draws


## effect size of x on m and cumulative probability for levels of m 

mx_beta <- 0.5
m_tau <- c(-1.1, 1, 1.5, 3)


## effects of x and m on y

yx_beta <- 0.2
ym_beta <- 0.8


# Ordered probit random draws

rop <- function(norm_distr, tau){
  op <- vector(mode = 'numeric', length = length(norm_distr))
  n_levels <- length(tau) + 1
  probs <- rep(NA, times = n_levels + 1)
  probs[1] <- 0
  probs[n_levels + 1] <- 1
  
  for(n in 1:length(norm_distr)){
    probs[2:n_levels] <- pnorm(tau - norm_distr[n])
    
    op[n] <- sample(1:n_levels, size = 1, prob = probs[2:(n_levels + 1)] - probs[1:n_levels])
  }
  
  op
}


## simulate dataframe

sim_df <- data.frame(row = 1:n_obs) %>%
  mutate(x = x,
         m = rop(mx_beta * x,
                 m_tau),
         y = rnbinom(1000, mu = exp(0 + yx_beta * x + ym_beta * m), size = 0.6))


## brms model

f1 <- bf(m ~ x, family = cumulative(probit))
f2 <- bf(y ~ m + x, family = negbinomial())
m1 <- brm(f1 + f2 + set_rescor(FALSE),
          data = sim_df,
          cores = 4,
          iter = 4000)

summary(m1)


## standardization from previous code:

###posterior samples and compute the indirect effect (mediation effect) a*b, as per Andrew Johnson approach

s1 <- as_draws_df(m1) 

sd_x <- sd(sim_df$x)
sd_y <- sd(sim_df$y)
sd_m <- sd(sim_df$m)


#calculate standardized a, the standardized x -> m pathway

a <- s1$b_m_x * sd_x / ( sqrt( s1$b_m_x^2 + sd_x^2 + 1 ) ) 


#calculate standardized b, the standardized m -> y pathway - mediator effect

## Q2. exp() in here somewhere?

b <- s1$b_y_m * sd_m / sd_y


#calculate a*b, the indirect effect (mediation effect)

ab <- a*b


## calculate the standardized direct effect, x -> y pathway 

## Q1. standardize the same way as for the mediator effect, b?

## Q2. exp() in here somewhere?

direct <- s1$b_y_x * sd_x / sd_y


median(a) ## standardized a

median(b) ## standardized b, mediator effect

median(ab) ## standardized indirect effect

median(direct) ## standardized direct effect

median(ab+direct) ## standardized total effect