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:
