One approach to this problem is to use the data manipulation tools in the merTools
package in order to create a data.frame that simulates the relationship you want to explore. Here is an example using the VerbAgg
data from the lme4
package, but you should be able to easily adopt it to your data and model.
data(VerbAgg)
fmVA <- glmer(r2 ~ (Anger + Gender + btype + situ)^2 +
(1|id) + (1|item), family = binomial,
data = VerbAgg)
We prep the data using the draw
function in merTools
. Here we
draw the average observation from the model frame. We then wiggle
the
data by expanding the dataframe to include the same observation repeated
but with different values of the variable specified by the var
parameter. Here, we expand the dataset to all values of btype
, situ
,
and Anger
subsequently.
# Select the average case
newData <- draw(fmVA, type = "average")
newData <- wiggle(newData, var = "btype", values = unique(VerbAgg$btype))
newData <- wiggle(newData, var = "situ", values = unique(VerbAgg$situ))
newData <- wiggle(newData, var = "Anger", values = unique(VerbAgg$Anger))
head(newData, 10)
r2 Anger Gender btype situ id item
1 N 20 F curse other 5 S3WantCurse
2 N 20 F scold other 5 S3WantCurse
3 N 20 F shout other 5 S3WantCurse
4 N 20 F curse self 5 S3WantCurse
5 N 20 F scold self 5 S3WantCurse
6 N 20 F shout self 5 S3WantCurse
7 N 11 F curse other 5 S3WantCurse
8 N 11 F scold other 5 S3WantCurse
9 N 11 F shout other 5 S3WantCurse
10 N 11 F curse self 5 S3WantCurse
In the next step we simply pass this new dataset to predictInterval
in order to generate predictions for these counterfactuals and we bind these predictions to the original data.
plotdf <- predictInterval(fmVA, newdata = newData, type = "probability",
stat = "median", n.sims = 1000)
plotdf <- cbind(plotdf, newData)
Then we plot the predicted values against the continuous variable, Anger
,
and facet and group on the two categorical variables situ
and btype
respectively.
ggplot(plotdf, aes(y = fit, x = Anger, color = btype, group = btype)) +
geom_point() + geom_smooth(aes(color = btype), method = "lm") +
facet_wrap(~situ) + theme_bw() +
labs(y = "Predicted Probability")
Which produces the following plot of fitted values:
