p values / confidence intervals for smooth terms

[charliercooper]

Hello :) I am trying to include a smooth term in my model, but am unsure where I can find a p-value or confidence intervals for this? It does not seem to be in the tidy tibbles for fixed or random effects? I did come across another post showing how to extract p-values (summary(fit$sd_report, select = "fixed", p.value = TRUE), am I correct in thinking the final term "ln_smooth_sigma" shows the p-value for my smooth term? Many thanks in advance! Charlie

[seananderson]

The ln_smooth_sigma term represents the log SD of the random effects that represent penalties on the smoother basis functions. So, as exp(ln_smooth_sigma) gets bigger, the basis functions start contributing more and more to the model fit and you get something that is wiggly. However, I don't think that has the same interpretation as the p-value in mgcv for a smoother, which is described in this paper. The implementation in sdmTMB is similar to in brms and I think gamm4. It uses mgcv::smooth2random(). Short answer is that as that ln_smooth_sigma term gets small the relationship between the predictor and response becomes linear, but there isn't a good p-value we've worked out to work with. I'd work with it visually. The p-value on ln_smooth_sigma probably doesn't make much sense since I assume it's reflecting whether it's different from zero, which means whether exp(ln_smooth_sigma) is different from 1, which doesn't mean much here. In practice that SD will tend to collapse to zero if the data don't support a wiggly function and standard errors will blow up.

If you're curious, you can inspect those elements like this:

library(sdmTMB)
mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 20)
fit <- sdmTMB(
  density ~ s(depth), 
  data = pcod_2011, 
  mesh = mesh, 
  family = tweedie(link = "log")
)

# basis function matrix:

head(fit$smoothers$Zs[[1]])
#>              [,1]        [,2]       [,3]        [,4]         [,5]        [,6]
#> [1,] -0.006094308 -0.14024402 -0.1186314  0.09694187  0.016469037  0.01086803
#> [2,]  0.011140592  0.36464171  0.1361885 -0.16215870 -0.081179418  0.03501230
#> [3,]  0.009321785  0.38152610  0.1323766 -0.16949601 -0.080817456  0.03879639
#> [4,]  0.006858794  0.05201337  0.1046212  0.01779886  0.002230583 -0.02047334
#> [5,] -0.038154705 -0.34876616 -0.4735616 -0.21523507  0.126339775  0.09679300
#> [6,] -0.140193869 -0.29346273 -0.5554432 -0.40324786  0.139291519  0.26437792
#>              [,7]         [,8]
#> [1,]  0.016414029  0.010448578
#> [2,] -0.048794117 -0.006254895
#> [3,] -0.048518369 -0.008718771
#> [4,]  0.002920848  0.003360885
#> [5,]  0.060572322 -0.038893021
#> [6,]  0.086184861 -0.177432912
est <- as.list(fit$sd_report, "Est")
se <- as.list(fit$sd_report, "Std. Error")

# parameters that are multiplied with the basis function matrix:

est$b_smooth
#>            [,1]
#> [1,]   3.194859
#> [2,]  -8.348903
#> [3,]  11.084558
#> [4,]  12.161908
#> [5,]  -4.574198
#> [6,] -15.515514
#> [7,]  12.324199
#> [8,]  -5.279749

# log SD of those random effects:

est$ln_smooth_sigma
#>          [,1]
#> [1,] 2.570412

^{Created on 2023-05-24 with reprex v2.0.2}