3

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:

enter image description here

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

Shawn Hemelstrand
  • 2,676
  • 4
  • 17
  • 30
  • Why do you think it is off? – Gavin Simpson Mar 20 '23 at 12:41
  • I'm just having trouble eyeballing what the polynomials should look like compared to the spline. In the first question it seemed more obvious. The parabolas that touch seem quite clear, and when they snake around the curves of the middle region that also seems right, but the left and right region don't seem immediately apparent that they are shaped by the basis functions. Perhaps that is just me being naive. – Shawn Hemelstrand Mar 20 '23 at 13:07
  • For example, what basis function helps form the far left side of the spline? Is it just the orange looking basis function? – Shawn Hemelstrand Mar 20 '23 at 13:10
  • The far left is defined by all those non-zero values of the basis functions; quite a few are ~ -1, so it doesn't need that many to sum to give the spline at -3 – Gavin Simpson Mar 20 '23 at 17:52
  • I felt like I already responded to this awhile ago, but thanks for clarifying. That seems to make sense. – Shawn Hemelstrand Apr 12 '23 at 22:04

0 Answers0