# Marginal effects for nonlinear multilevel models

**URL:** <https://discourse.mc-stan.org/t/marginal-effects-for-nonlinear-multilevel-models/36050>\
**Category:** Modeling\
**Tags:** cmdstanpy\
**Created:** [July 29, 2024, 7:09pm UTC](https://discourse.mc-stan.org/t/marginal-effects-for-nonlinear-multilevel-models/36050 "2024-07-29T19:09:35Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![amillb](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/amillb/32/18914_2.png) [@amillb](https://discourse.mc-stan.org/u/amillb)\
**Post date:** [July 29, 2024, 7:09pm UTC](https://discourse.mc-stan.org/t/marginal-effects-for-nonlinear-multilevel-models/36050/1 "2024-07-29T19:09:35Z")

</div>

I’m trying to calculate marginal effects for a multilevel nonlinear (beta) model using cmdstanpy. I have continuous predictors at both the top and lower levels of the model, [roughly following this example](https://mc-stan.org/docs/2_18/stan-users-guide/multivariate-hierarchical-priors-section.html).

Marginal effects of the lower-level coefficients seem pretty straightforward, either through feeding fake data to `model.generate_quantities()`, or computing them directly in the generated quantities block (see code below).

For marginal effects of top-level coefficients, what is the best way to do this? Fake data doesn’t work, as the generated quantities block takes the lower-level coefficients (betas) as fixed, but the top-level predictors work through the betas. All the examples I see ([e.g. here](https://www.andrewheiss.com/blog/2021/11/10/ame-bayes-re-guide/)) assume that there are no top-level predictors (only groups).

My attempt (somewhat convoluted) is to recompute the betas in generated quantities. See below. Is there a better way to do it?

Also, is there any equivalent to `brms` for Python, or another package that allows the easy computation of marginal effects?

Thank you!

```no-highlight
data {
  int<lower=0> J; 
  int<lower=0> N;
  int<lower=0> L; // number of lower-level covariates (excl. intercept)
  int<lower=0> M; // number of top-level covariates (excl. intercept)
  array[N] int<lower=1,upper=J> group_num;
  vector[J] Z1; // use length J rather than N
  vector[N] X1;
  vector[N] y;

} 

transformed data {
    // roughly following https://mc-stan.org/docs/stan-users-guide/multivariate-hierarchical-priors.html
    // individual predictor matrix
    matrix[N, L+1] X;
    for (i in 1:N) {
        X[i,1] = 1;
        X[i,2] = X1[i];          
    } 
    // group-level predictors
    array[J] row_vector[M+1] Z;
    for (j in 1:J) {
        Z[j][1] = 1.0;
        Z[j][2] = Z1[j];
    } 
}

parameters {
  corr_matrix[(L+1)] Omega; // decomposed correlation matrix for individual-level predictors. 
  vector<lower=0>[L+1] tau; // scale of diagonals
  matrix[M+1,L+1] gamma; // top-level coefficients (including the intercept)
  array[J] vector[L+1] beta; // lower-level level coefficients (including the intercept)
  real<lower = 0> phi; // dispersion

} 

transformed parameters {
  vector[N] mu; 
  vector[N] xb;
  
  for (n in 1:N) {
      xb[n] = X[n] * beta[group_num[n]];
  }

  mu = inv_logit(xb);

}

model {
  // Weakly informative priors for covariance matrix and second-level parameters
  tau ~ cauchy(0,0.5);
  Omega ~ lkj_corr(2.0);
  to_vector(gamma) ~ cauchy(0, 1);
  phi ~ cauchy(0, 30);           
      
  {
    array[J] row_vector[L+1] z_gamma;
    for (j in 1:J) {
      z_gamma[j] = Z[j] * gamma;
    }
    beta ~ multi_normal(z_gamma, quad_form_diag(Omega, tau));
  }

  y ~ beta_proportion(mu, phi);

  }

generated quantities {
  vector[N] y_hat;        

  // marginal effects
  array[N] vector[L+1] me;   
  array[N] vector[L+1] me_group;
  
  for (n in 1:N) { 
    y_hat[n] = inv_logit(X[n] * beta[group_num[n]]);
   }

  // marginal effects of lower-level variables
  for (i in 1:L+1) {
      real yh2;
  
      for (n in 1:N) {
         yh2 = inv_logit(X[n] * beta[group_num[n]] + 0.0001*beta[group_num[n]][i] );
         me[n][i] = (yh2 - y_hat[n] ) / 0.0001;
      }
   }

   // marginal effects of top-level variables. 

   { real yh2;

     for (m in 1:M+1) { // loop over top-level covariates
     	 array[J] vector[L+1] me_beta = beta;          
         for (i in 1:L+1) {
	         for (j in 1:J) {  
               me_beta[j][i] = me_beta[j][i] + 0.0001 * gamma[m, i]; // add epsilon to top-level variables, and propagate that through the betas; equivalent to adding 0.0001 to Z[j]
            }
          }

          for (n in 1:N) {
             yh2 = inv_logit(X[n] * me_beta[group_num[n]]); 
	         me_group[n][m] = (yh2 - y_hat[n] ) / 0.0001;
          }

        
     }
   }

}

```

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 30, 2024, 6:11pm UTC](https://discourse.mc-stan.org/t/marginal-effects-for-nonlinear-multilevel-models/36050/2 "2024-07-30T18:11:53Z")

</div>

> [@amillb](#):
>
> marginal effects

Hi, @amillb, and welcome to the Stan forums.

Can you say a bit more about what you mean by “marginal effects”? I’m not familiar with the term, but it looks like you’re generating the variable called `me` with

```stan
 me[n][i] = (yh2 - y_hat[n] ) / 0.0001;

```

where it looks like `y_hat` is some kind of expectation for a logistic regression outcome. I’m also not sure why there is rescaling by 0.0001. Usually you want to eliminate that and work on natural scales.

Also, something like a `cauchy(0, 30)` isn’t giving you much of a prior compared to just a flat prior. You’ll find things sample better with weakly informative priors that don’t allow huge outliers like the Cauchy.

You can also make a lot of this code more efficient with vectorization, e.g., by defining `yh2` for `inv_logit(X * me_beta[group_num])` to turn this into a much more efficient matrix multiply. Similarly, you can set `me_group[, m] = (y2 - y_hat) / 0.0001`, and so on. It’s more efficient because it reduces the size of the autodiff expression graph, not because it eliminates a loop.

---

<div class="post-metadata">

**Author:** ![amillb](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/amillb/32/18914_2.png) [@amillb](https://discourse.mc-stan.org/u/amillb)\
**Post date:** [August 5, 2024, 9:16pm UTC](https://discourse.mc-stan.org/t/marginal-effects-for-nonlinear-multilevel-models/36050/3 "2024-08-05T21:16:35Z")

</div>

Hi @Bob_Carpenter,

Many thanks for replying; the suggestions for priors and for code efficiency are also appreciated.

By marginal effect, I mean the change in the dependent variable that results from a “small" change in a specified independent variable (i.e., the derivative). [See, for example, the discussion here.](https://www.andrewheiss.com/blog/2021/11/10/ame-bayes-re-guide/)

The rescaling by 0.0001 is because 0.0001 is added to the independent variable in the following line, so the division then approximates the derivative. Because the effect of the top-level variables is channeled through the lower-level coefficients `beta`, I call these adjusted coefficients `me_beta`.

```
    me_beta[j][i] = me_beta[j][i] + 0.0001 * gamma[m, i];

```

`yhat` and `yh2` are both the predicted outcome (y), based on a beta functional form. `yhat` is the predicted outcome for the actual values of `X` and `Z`. `yh2` is the predicted outcome with 0.0001 added to one of the group-level predictors `Z[m]` for every observation. The marginal effect is the difference between the two, scaled by 0.0001. Does that make sense?

Adam

PS I have lurked in the Stan forums for many years and found the discussions very helpful. Until now, my questions have been asked and answered already.