How can I prevent R from extrapolating the regression line outside the range of data points using a mixed model?

I am fitting a mixed model in R with an interaction term as follows: model <- lmer(y ~ x * z + (1|h), data = df). y and x are continuous variables and z is a categorical variable, and h is a random effect.

Data:

 n <- 50
df <- tibble::tibble(
  x = rnorm(n),
  z = sample(-1:1, n, replace = TRUE),
  h = sample(1:3, n, replace=TRUE), # random effect
  y = z*(0.5*x) + rnorm(n,0,0.2))
df$z <- as.factor(df$z)
df$h <- as.factor(df$h)

Then, I'm using the plot_predictions() function from the marginaleffects package, then I use facet_wrap(~ z) from the ggplot2 to visualize each regression separately of the different categories as follows:

plot_predictions(model, condition = c("x","z"), vcov = T, points= 0.3) + facet_wrap(~ z)

As you can notice the first and second plots extrapolate the regression line outside the range of their data according to the limits of data points of the third plot.

So, I try to control this by adding newdata = df:

plot_predictions(model, new data = df, by = c("x","z"), vcov = T, points= 0.3) + facet_wrap(~ z)

however, now I have a new problem; I lose the linear effect and it looks like a disruptive non-continuous line. Interestingly, it only happens when I fit a mixed model with random effects (i.e., glmm, hgam) but this is not the case for a simple linear model.

See plot here

Any idea how I can prevent this from happening, constraining my regression line (+- CIs) to the range of the data points for each panel and keeping the linear model line?

The code below seems to work as expected using the latest version of marginaleffects.

(PS: Next time, it would be convenient for people who try to give answers if you wrote a complete replicable example with libraries and a continuous code block that we can jut cut and paste.)

library(lme4)
library(ggplot2)
library(marginaleffects)
n <- 50
df <- tibble::tibble(
  x = rnorm(n),
  z = sample(-1:1, n, replace = TRUE),
  h = sample(1:3, n, replace=TRUE), # random effect
  y = z*(0.5*x) + rnorm(n,0,0.2))
df$z <- as.factor(df$z)
df$h <- as.factor(df$h)
model <- lmer(y ~ x * z + (1|h), data = df)

plot_predictions(model, newdata = df, by = c("x","z"), vcov = TRUE, points= 0.3) +
    facet_wrap(~ z, scales = "free")