Overlaying Basis Functions for GAM Plot Part II: Random Effects

Question

I recently posted a question about overlaying basis function on a GAM plot. This was quite helpful for modeling fixed effects models. However, I want to do the same thing for GAMMs and I'm running into a bit of a road block. I've tried simulating some data to achieve this in a way as similar as possible to the original question. The main difference is that I used a typical crossed random effects design, adding in subjects and items to create a GAMM later.

#### Load Libraries and Set Theme ####
library(mgcv)
library(tidyverse)
library(gratia)
theme_set(theme_bw())

#### Simulate Data ####
dat <- data_sim("eg2", seed = 123, n = 1000)
cross.dat <- dat %>% 
  mutate(subject = factor(rbinom(n=1000,size=100,prob=.5)),
         item = factor(rbinom(n=1000,size=50,prob=.5)))
cross.dat

Which looks like this:

# A tibble: 1,000 × 6
         y      f      x      z subject item 
     <dbl>  <dbl>  <dbl>  <dbl> <fct>   <fct>
 1 -1.36   0.631  0.288  0.274  49      24   
 2 -1.75   0.326  0.788  0.594  55      30   
 3  0.346  0.382  0.409  0.160  55      23   
 4  0.0416 0.306  0.883  0.853  42      26   
 5 -4.86   0.234  0.940  0.848  55      18   
 6  2.51   0.428  0.0456 0.478  55      26   
 7  0.873  0.374  0.528  0.774  50      26   
 8  4.90   0.0642 0.892  0.295  52      30   
 9  1.50   0.134  0.551  0.0656 44      27   
10 -0.502  0.392  0.457  0.441  53      26   
# … with 990 more rows
# ℹ Use `print(n = ...)` to see more rows

Then I plotted the data in a similar way to the original question like so:

#### Fit GAMM ####
m <- gam(
  y 
  ~ s(x, bs = "cr")
  + s(z, bs = "cr")
  + ti(x, z, bs = "cr")
  + s(subject, bs = "re")
  + s(item, bs = "re"),
  data = cross.dat,
  method = "REML"
)
plot(m, select = 2)

#### Data Slice and Basis ####
ds <- data_slice(m, z = evenly(z, n = 200))
bs <- basis(m, term = "s(z)", data = ds)

#### Sum Basis Functions at Each Value ####
spl <- bs |>
  group_by(z) |>
  summarise(spline = sum(value))

#### Plot ####
bs |> 
  ggplot(aes(x = z,
             y = value, 
             colour = bf, 
             group = bf)) +
  geom_line(show.legend = FALSE) +
  geom_line(aes(x = z, 
                y = spline), 
            data = spl,
            linewidth = 1.5,
            inherit.aes = FALSE) +
  labs(y = expression(f(z)), 
       x = "z",
       title = "Simulated Basis Functions")

This looks to be close, but for some reason it looks off to me:

Is this correct? If not how do I fix it?