Page **1** of **1**

### Interaction plots for mlwinfitMCMC fixed effects

Posted: **Wed Apr 07, 2021 4:05 pm**

by **wfleming**

Any recommendation for the best way to plot fixed interaction effects from mlwinfitMCMC objects?

Model is simple like:

Code: Select all

`m <- runmlwin(y ~ 1 + lev1x*lev1x + (1 | lev2id) + (1 | lev1id), estoptions = list(EstM = 1))`

and want to plot the lev1x*lev1x

Would usually use plot_model from 'sjPlot' for lm objects but naturally not written for mlwin objects. Would also do it manually extracting predicted values but predict() doesn't give the intervals for mlwin objects.

Any ideas greatly appreciated.

### Re: Interaction plots for mlwinfitMCMC fixed effects

Posted: **Fri Apr 09, 2021 3:43 pm**

by **ChrisCharlton**

You could manually create the data required for this plot by generating a chain of predictions and then calculating the mean and credible intervals from this. An example of this would be as follows:

Code: Select all

```
library(R2MLwiN)
library(coda)
data(tutorial, package = "R2MLwiN")
(model1 <- runMLwiN(normexam ~ 1 + standlrt * vrband + (1 | school) + (1 | student), estoptions = list(EstM = 1), data = tutorial))
# Data frame for variables used in fixed part of the model
modeldata <- model.matrix(normexam ~ 1 + standlrt * vrband, data = tutorial)
head(modeldata)
# Data frame containing parameter chains (specify variables in the order corresponding with the data frame above)
modelchains <- as.matrix(model1@chains)[, c("FP_Intercept", "FP_standlrt", "FP_vrbandvb2", "FP_vrbandvb3", "FP_standlrt:vrbandvb2", "FP_standlrt:vrbandvb3")]
# Prediction chain
predictchain <- as.mcmc(modelchains %*% t(modeldata))
# Mean prediction
prediction <- drop(apply(predictchain, 2, mean))
# 95% credible intervals
quantiles <- t(apply(predictchain, 2, quantile, c(0.025, 0.975)))
```

### Re: Interaction plots for mlwinfitMCMC fixed effects

Posted: **Tue Apr 20, 2021 2:05 pm**

by **wfleming**

Thanks Chris this worked well. Wouldn't have worked it out without you.

Only change I made was to use bayestestR::hdi() to get credible intervals instead of the usual quantiles.