CODEWITHSUNDEEP
rregressioninteractionordinal
Rhys Maredudd Davies
Rhys Maredudd Davies

Reputation: 33

How to plot an interaction effect within an ordinal regression model?

I have an ordinal model for predicting anxiety severity, using the clm() function. Within the model there is significant interaction effect between two of the variables. I am looking for a way to demonstrate the interaction effect visually, so I can interpret the effect with greater clarity. I have tried looking for different options, but I can't find any applications that work with the clm() function of modelling ordinal regressions. The error I get for the example below is "Error in UseMethod("family") : no applicable method for 'family' applied to an object of class 'clm'". Modelling the ordinal regression with the polr function gives the same response.

m <- clm(anxiety_levels ~ pred_a * pred_b + pred_c, data, link = "logit")
interact_plot(m, pred = pred_a, modx = pred_b)

Any suggestions on how to plot the interaction from an ordinal regression would be greatly appreciated.

Upvotes: 1

Views: 1656

Answers (1)

Vincent
Vincent

Reputation: 17823

The marginaleffects package supports clm models. (Disclaimer: I am the author.)

However, note that the development version includes an important bug fix. You can install it like this:

library(remotes)
install_github("vincentarelbundock/marginaleffects")

Make sure you restart R completely, then we can estimate, summarize and plot the results. Note that we use facet_wrap(~group) often to present results for each outcome level separately in the plots:

library(ordinal)
library(ggplot2)
library(marginaleffects)

dat <- data.frame(
    y = factor(sample(letters[1:5], 100, replace = TRUE)),
    x1 = rnorm(100),
    x2 = rnorm(100))
mod <- clm(y ~ x1 * x2, data = dat)

Adjusted predictions

plot_cap(mod, condition = "x1") + facet_wrap(~group)



plot_cap(mod, condition = c("x1", "x2")) + facet_wrap(~group)

Marginal effects (slope)

mfx <- marginaleffects(mod)
summary(mfx)
#> Average marginal effects 
#>    Group Term     Effect Std. Error  z value Pr(>|z|)    2.5 %   97.5 %
#> 1      a   x1 -0.0066902   0.021233 -0.31508  0.75270 -0.04831 0.034927
#> 2      a   x2 -0.0018825   0.021849 -0.08616  0.93134 -0.04471 0.040941
#> 3      b   x1 -0.0058553   0.017050 -0.34343  0.73128 -0.03927 0.027561
#> 4      b   x2 -0.0022098   0.017486 -0.12637  0.89944 -0.03648 0.032062
#> 5      c   x1 -0.0017640   0.004803 -0.36729  0.71341 -0.01118 0.007649
#> 6      c   x2 -0.0015586   0.004723 -0.33002  0.74138 -0.01081 0.007698
#> 7      d   x1  0.0031011   0.010273  0.30186  0.76276 -0.01703 0.023236
#> 8      d   x2  0.0006885   0.010593  0.06499  0.94818 -0.02007 0.021449
#> 9      e   x1  0.0112084   0.031057  0.36089  0.71818 -0.04966 0.072080
#> 10     e   x2  0.0049625   0.031719  0.15645  0.87568 -0.05721 0.067130
#> 
#> Model type:  clm 
#> Prediction type:  response

plot(mfx) + facet_wrap(~group)



plot_cme(mod, effect = "x1", condition = "x2") + facet_wrap(~group)

Upvotes: 3

Related Questions