从具有多个“by”平滑项的 gam (mgcv) 进行预测

Question

我安装了一个 gam (mgcv)，其中 y 通过因子 x1 和因子 x2 以单独的平滑项建模为时间 (t) 的函数。我还包括了 r，它应该模仿随机因子。

library(mgcv)
dat <- data.frame(y = rnorm(1500,500,100),
                  t = rep(c(1:5),100),
                  x1 = as.factor(rep(letters[1:3],length.out=1500)),
                  x2 = as.factor(rep(LETTERS[1:5],length.out=1500)),
                  r = as.factor(rep(letters[4:26],length.out=1500)))


fit <- gam(y~s(t,by=x1,k=4)+s(x1,bs="re")+s(t,by=x2,k=4)+s(x2,bs="re")+s(r,bs="re"),data=dat,
           family="gaussian",method="REML")

我现在想要预测 x2 的所有水平上 x1 水平的 y。然而，x2 中的新水平将导致预测中出现 NA。

pred <- data.frame(predict(fit,type="response",newdata = data.frame(t=seq(0,5,0.1),
                                                                    x1="a",
                                                                    x2="new_level",
                                                                    r="new_level")
                           ,exclude = c("x2","r"),se.fit = T))

pred$fit

如果我将 x2 设置为 x2 中包含的级别之一，即使我指定将其排除，它也会预测将根据我使用的值/级别而变化的值。有趣的是，r 达到了新的水平，但没有导致 NA。所以，我认为这与“by”平滑术语有关。但是，我无法想出解决方案。 一种可能是拟合两个模型，其中仅包含其中一项平滑项 w 每个？这样可以吗？

Answer 1

听起来您只需要正确指定

exclude

即可。例如，

library("gratia")
library("mgcv")
su_eg4 <- data_sim("eg4", n = 400,  dist = "normal", scale = 2, seed = 1)
m <- gam(y ~ fac + s(x2, by = fac) + s(x0, by = fac),
         data = su_eg4, method = "REML")

该模型具有以下平滑项

> smooths(m)
[1] "s(x2):fac1" "s(x2):fac2" "s(x2):fac3" "s(x0):fac1" "s(x0):fac2" "s(x0):fac3"

因此，如果您想在此处排除所有包含

s(x0)

的平滑，则需要使用它们的全名。因子

by

项按照此模板命名：

s(<x>):<f><level>

，其中

• <x>

  是平滑协变量，

• <f>

  是因子的名称

• <level>

  是此平滑所针对的

  <f>

  的级别。

如果你有很多关卡，把它们全部写出来可能会很乏味，所以我们可以搜索相关的字符串

sms <- smooths(m)
excl <- grepl("x0", sms)

然后相应地调整

exclude

参数值

exclude = sms[excl]

作为

sms[excl]

评估为

> sms[excl]
[1] "s(x0):fac1" "s(x0):fac2" "s(x0):fac3"

在你的情况下，你似乎想使用

sms <- smooths(fit)
excl <- grepl("x2", sms) | grepl("s\\(r\\)", sms)

捕获

r

的随机截距以及涉及

x2

的任何平滑。

---

我实际上认为你的模型会更简单地指定使用：

m <- <- gam(y ~
    s(t) + # "average" smooth of t
    s(t, x1, bs = "sz") + # smooth difference from `s(t)` for each level of x1
    s(t, x2, bs = "sz") + # smooth difference from `s(t)` for each level of x2
    s(r, bs = "re"), # random intercept term for `r`
  family = "gaussian", method = "REML")

然后你就可以使用

exclude = c("s(t,x2)", "s(r)")

排除您想要排除的平滑。

请注意，使用

sz

基础时，组均值包含 于基础中，因此所示模型中不包含

x1

或

x2

或其 ranef 等效项的参数效应。