mgcv：修复 GAMM 中的平滑参数和模型嵌套的有效性
mgcv: Fix smoothing parameter in GAMM and validity of model nesting
我目前正在尝试使用 mgcv 包中的 gamm 对协变量 x（年龄）和性别 0-1 变量之间的相互作用进行建模。在为每个性别指定一个具有平滑项的主模型（我们称之为 M0）之后，我想检验更简单的假设，即性别之间的差异是线性的（而不是任意平滑的）。我有以下两个问题：
- 当尝试正确嵌套模型时，我想从
M0 中取出 0 性别平滑的平滑参数，并在更简单的模型规范中使用它。但是我收到以下错误消息：
Error in gamm.setup(gp, pterms = pTerms, data = mf, knots = knots, parametric.only = FALSE, :
gamm can not handle linked smoothing parameters (probably from use of `id' or adaptive smooths)
- 第二个问题比较蠢。当我从每个性别的平滑到 0 性别平滑和线性差异到 1 性别时，模型是否甚至嵌套？
下面是一个例子。我模拟了一些随机数据，所以数据没有显示我实际数据的行为，但问题仍然存在。
library(mgcv)
### Simulate random data
x <- rnorm(100, mean = 10, sd = 1.5)
y <- rnorm(100, mean = 1, sd = 0.025)
id <- sample(1:10, size = 100, replace = T)
id <- as.factor(id)
gender <- sample(c(0,1), size = 100, replace = T)
### Specify main model, M0
ctrl <- list(niterEM=0,optimMethod="L-BFGS-B", msMaxIter = 100)
M0 <- gamm(y~s(x, by = as.factor(gender)) + gender,
random=list(id=~1+x), control=ctrl)
### Specify model with linear difference between gender0 and gender1
M1 <- gamm(y~s(x) + gender:x + gender,
random=list(id=~1+x), control=ctrl)
### Testing
anova(M0$lme, M1$lme)
### Problems:
sp0 <- M0$gam$sp[1]
M1 <- gamm(y~s(x, sp = sp0) + gender:x + gender,
random=list(id=~1+x), control=ctrl)
Error in gamm.setup(gp, pterms = pTerms, data = mf, knots = knots, parametric.only = FALSE, :
gamm can not handle linked smoothing parameters (probably from use of `id' or adaptive smooths)
有什么想法吗？提前致谢。
关于 gamm 错误
这是一件很有意思的事情……嗯，还是先说明一下逻辑吧
原则上固定gamm中的平滑参数是非法的，因为gamm会将平滑的摆动分量视为随机效果， lme 估计的方差（因为你有高斯数据）。如果您尝试修复平滑参数，这实际上是说您想要修复随机效应的方差。嗯，lme 不允许你这样做（而且我怀疑这种尝试在混合建模中是否合法。）
因此 gamm 将禁用任何可能的平滑参数约束，包括：
- 平滑参数的下限，通过
min.sp；
- 链接平滑共享相同的平滑参数，通过
s() 中的 id；
- 通过
sp、s()直接指定平滑参数。
前两个完美检查。 gamm 没有像 gam 这样的 min.sp 参数；即使你通过 ... 指定它，它也没有机会被使用（因为后来它是 NULL 在 gamm.setup 期间传递给 gam.setup，所以你指定的 min.sp 被忽略）。 id 的规范也将通过您看到的错误消息检测到，但是当然您没有指定 id 所以上面的错误在这里没有报告正确的问题，因此是一个错误。
第三个，其实还没有通过gamm直接查过。理想情况下，一旦 gamm / gam 公式被解释（通过 interpret.gam），sp 应该重置为 -1 如果它不容易是, 并应发出有关此的警告消息。但是，缺少这一部分。因此，目前你能做的最好的事情就是不指定 sp.
你有有效的模型嵌套吗？
现在让我们谈谈您对嵌套的关注。 嵌套与基础设置有关，而不是与平滑参数的选择有关。只要您有相同的基础集（相同的基础类型、相同的“节点”数量和/或位置） ，模型矩阵将相同。例如，在您的模型 M0 和 M1 中，您具有与 mgcv 默认 bs = 'tp', k = 10 相同的 s() 配置。因此 s() 的设计矩阵在您的两个模型中是相同的。 by = factor(gender) 只是将此 s() 复制到主模型 M0 中 gender 的所有级别。可能不容易看到，其实你的M1是嵌套在M0.
让我们考虑一个小例子来验证这一点。为了简单起见，我不会使用 mgcv 中的 s(x)，而是使用 poly(x, degree = 2)（假设是 s(x)）。让我们先生成一些玩具数据：
x <- 1:10
f <- gl(2, 5, labels = c("M", "F"))
由于 f 不是有序因子，s(x, by = factor(f)) 通过为 f 的所有级别复制 s(x) 生成设计矩阵：
## original design matrix for `s(x)`
X0 <- poly(x, 2)
## design matrix for `f`, without contrasting
Xf <- model.matrix(~ f + 0)
## design matrix for level `M`
X1 <- X0 * Xf[, 1]
## design matrix for level `F`
X2 <- X0 * Xf[, 2]
## design matrix for `s(x, by = f)` "please, imagine it as `poly`"
X <- cbind(X1, X2)
# 1 2 1 2
# [1,] -0.49543369 0.52223297 0.00000000 0.00000000
# [2,] -0.38533732 0.17407766 0.00000000 0.00000000
# [3,] -0.27524094 -0.08703883 0.00000000 0.00000000
# [4,] -0.16514456 -0.26111648 0.00000000 0.00000000
# [5,] -0.05504819 -0.34815531 0.00000000 0.00000000
# [6,] 0.00000000 0.00000000 0.05504819 -0.34815531
# [7,] 0.00000000 0.00000000 0.16514456 -0.26111648
# [8,] 0.00000000 0.00000000 0.27524094 -0.08703883
# [9,] 0.00000000 0.00000000 0.38533732 0.17407766
#[10,] 0.00000000 0.00000000 0.49543369 0.52223297
你的第二个模型M1，只有一个平滑项s(x)，它的设计矩阵是X0。
我们可以通过以下方式看到您的 M1 嵌套在 M0 中：
- 由于
X1 + X2 = X0，s(x)和s(x, by = f)的设计矩阵具有相同的列跨度，因此s(x)嵌套在s(x, by = f)；
- 因为
x 嵌套在 s(x) 中，所以 x:f 嵌套在 s(x, by = f) 中。
我会做什么
虽然您的模型已经很好地嵌套，但主模型 M0 与您的 M1 没有可比较的解释。您的主模型 M0 最终将对每个级别进行独立的平滑处理，而您的 M1 则专注于两组之间的差异。
如果能控制M0承认一种“参考平滑+差异平滑”的形式就好了。然后，如果差异平滑变成一条线，而不是实际拟合 M1 我们已经知道没有证据表明组之间存在非线性差异。
在mgcv中，如果你的因子是有序的，将构建差异平滑。所以我建议你通过以下方式拟合你的主要模型：
gender1 <- ordered(gender) ## create an ordered factor
s(x) + s(x, by = gender1) + gender
如果估计结果显示差异平滑 s(x, by = gender1) 为一条线，您知道您可以改为拟合更简单的模型
s(x) + gender:x + gender
即使不使用 F-test.
注意，有一个有序的因子by是非常重要的，这样才能平滑地构造“差异”。如果你这样做
s(x) + s(x, by = gender) + gender ## note, it is "gender" in "by"
s(x) 和 s(x, by = gender) 完全线性相关。生成的模型矩阵将是 rank-deficient.
一些跟进
I forgot to include in my example that I at first compared the same model parametrized as s(x, by = as.factor(gender)) and s(x) + s(x, by = gender) by AIC (recall gender is 0-1 numerical variable). These models are statistically equivalent, but the smoothing parameters are obviously estimated differently in the cases, and the AIC thus differ a bit.
哦，是的。您的 gender 是二进制的，因此数值 by 也是构建差异平滑的好主意。但是要小心。数值 by 不会产生居中平滑。因此，s(x) + s(x, by = gender) 将隐含一个截距列，与模型截距混淆。你应该选择 s(x) + s(x, by = gender) - 1.
我目前正在尝试使用 mgcv 包中的 gamm 对协变量 x（年龄）和性别 0-1 变量之间的相互作用进行建模。在为每个性别指定一个具有平滑项的主模型（我们称之为 M0）之后，我想检验更简单的假设，即性别之间的差异是线性的（而不是任意平滑的）。我有以下两个问题：
- 当尝试正确嵌套模型时，我想从
M0中取出 0 性别平滑的平滑参数，并在更简单的模型规范中使用它。但是我收到以下错误消息：
Error in gamm.setup(gp, pterms = pTerms, data = mf, knots = knots, parametric.only = FALSE, :
gamm can not handle linked smoothing parameters (probably from use of `id' or adaptive smooths)
- 第二个问题比较蠢。当我从每个性别的平滑到 0 性别平滑和线性差异到 1 性别时，模型是否甚至嵌套？
下面是一个例子。我模拟了一些随机数据，所以数据没有显示我实际数据的行为，但问题仍然存在。
library(mgcv)
### Simulate random data
x <- rnorm(100, mean = 10, sd = 1.5)
y <- rnorm(100, mean = 1, sd = 0.025)
id <- sample(1:10, size = 100, replace = T)
id <- as.factor(id)
gender <- sample(c(0,1), size = 100, replace = T)
### Specify main model, M0
ctrl <- list(niterEM=0,optimMethod="L-BFGS-B", msMaxIter = 100)
M0 <- gamm(y~s(x, by = as.factor(gender)) + gender,
random=list(id=~1+x), control=ctrl)
### Specify model with linear difference between gender0 and gender1
M1 <- gamm(y~s(x) + gender:x + gender,
random=list(id=~1+x), control=ctrl)
### Testing
anova(M0$lme, M1$lme)
### Problems:
sp0 <- M0$gam$sp[1]
M1 <- gamm(y~s(x, sp = sp0) + gender:x + gender,
random=list(id=~1+x), control=ctrl)
Error in gamm.setup(gp, pterms = pTerms, data = mf, knots = knots, parametric.only = FALSE, :
gamm can not handle linked smoothing parameters (probably from use of `id' or adaptive smooths)
有什么想法吗？提前致谢。
关于 gamm 错误
这是一件很有意思的事情……嗯，还是先说明一下逻辑吧
原则上固定gamm中的平滑参数是非法的，因为gamm会将平滑的摆动分量视为随机效果， lme 估计的方差（因为你有高斯数据）。如果您尝试修复平滑参数，这实际上是说您想要修复随机效应的方差。嗯，lme 不允许你这样做（而且我怀疑这种尝试在混合建模中是否合法。）
因此 gamm 将禁用任何可能的平滑参数约束，包括：
- 平滑参数的下限，通过
min.sp； - 链接平滑共享相同的平滑参数，通过
s()中的id； - 通过
sp、s()直接指定平滑参数。
前两个完美检查。 gamm 没有像 gam 这样的 min.sp 参数；即使你通过 ... 指定它，它也没有机会被使用（因为后来它是 NULL 在 gamm.setup 期间传递给 gam.setup，所以你指定的 min.sp 被忽略）。 id 的规范也将通过您看到的错误消息检测到，但是当然您没有指定 id 所以上面的错误在这里没有报告正确的问题，因此是一个错误。
第三个，其实还没有通过gamm直接查过。理想情况下，一旦 gamm / gam 公式被解释（通过 interpret.gam），sp 应该重置为 -1 如果它不容易是, 并应发出有关此的警告消息。但是，缺少这一部分。因此，目前你能做的最好的事情就是不指定 sp.
你有有效的模型嵌套吗？
现在让我们谈谈您对嵌套的关注。 嵌套与基础设置有关，而不是与平滑参数的选择有关。只要您有相同的基础集（相同的基础类型、相同的“节点”数量和/或位置） ，模型矩阵将相同。例如，在您的模型 M0 和 M1 中，您具有与 mgcv 默认 bs = 'tp', k = 10 相同的 s() 配置。因此 s() 的设计矩阵在您的两个模型中是相同的。 by = factor(gender) 只是将此 s() 复制到主模型 M0 中 gender 的所有级别。可能不容易看到，其实你的M1是嵌套在M0.
让我们考虑一个小例子来验证这一点。为了简单起见，我不会使用 mgcv 中的 s(x)，而是使用 poly(x, degree = 2)（假设是 s(x)）。让我们先生成一些玩具数据：
x <- 1:10
f <- gl(2, 5, labels = c("M", "F"))
由于 f 不是有序因子，s(x, by = factor(f)) 通过为 f 的所有级别复制 s(x) 生成设计矩阵：
## original design matrix for `s(x)`
X0 <- poly(x, 2)
## design matrix for `f`, without contrasting
Xf <- model.matrix(~ f + 0)
## design matrix for level `M`
X1 <- X0 * Xf[, 1]
## design matrix for level `F`
X2 <- X0 * Xf[, 2]
## design matrix for `s(x, by = f)` "please, imagine it as `poly`"
X <- cbind(X1, X2)
# 1 2 1 2
# [1,] -0.49543369 0.52223297 0.00000000 0.00000000
# [2,] -0.38533732 0.17407766 0.00000000 0.00000000
# [3,] -0.27524094 -0.08703883 0.00000000 0.00000000
# [4,] -0.16514456 -0.26111648 0.00000000 0.00000000
# [5,] -0.05504819 -0.34815531 0.00000000 0.00000000
# [6,] 0.00000000 0.00000000 0.05504819 -0.34815531
# [7,] 0.00000000 0.00000000 0.16514456 -0.26111648
# [8,] 0.00000000 0.00000000 0.27524094 -0.08703883
# [9,] 0.00000000 0.00000000 0.38533732 0.17407766
#[10,] 0.00000000 0.00000000 0.49543369 0.52223297
你的第二个模型M1，只有一个平滑项s(x)，它的设计矩阵是X0。
我们可以通过以下方式看到您的 M1 嵌套在 M0 中：
- 由于
X1 + X2 = X0，s(x)和s(x, by = f)的设计矩阵具有相同的列跨度，因此s(x)嵌套在s(x, by = f)； - 因为
x嵌套在s(x)中，所以x:f嵌套在s(x, by = f)中。
我会做什么
虽然您的模型已经很好地嵌套，但主模型 M0 与您的 M1 没有可比较的解释。您的主模型 M0 最终将对每个级别进行独立的平滑处理，而您的 M1 则专注于两组之间的差异。
如果能控制M0承认一种“参考平滑+差异平滑”的形式就好了。然后，如果差异平滑变成一条线，而不是实际拟合 M1 我们已经知道没有证据表明组之间存在非线性差异。
在mgcv中，如果你的因子是有序的，将构建差异平滑。所以我建议你通过以下方式拟合你的主要模型：
gender1 <- ordered(gender) ## create an ordered factor
s(x) + s(x, by = gender1) + gender
如果估计结果显示差异平滑 s(x, by = gender1) 为一条线，您知道您可以改为拟合更简单的模型
s(x) + gender:x + gender
即使不使用 F-test.
注意，有一个有序的因子by是非常重要的，这样才能平滑地构造“差异”。如果你这样做
s(x) + s(x, by = gender) + gender ## note, it is "gender" in "by"
s(x) 和 s(x, by = gender) 完全线性相关。生成的模型矩阵将是 rank-deficient.
一些跟进
I forgot to include in my example that I at first compared the same model parametrized as
s(x, by = as.factor(gender))ands(x) + s(x, by = gender)by AIC (recallgenderis 0-1 numerical variable). These models are statistically equivalent, but the smoothing parameters are obviously estimated differently in the cases, and the AIC thus differ a bit.
哦，是的。您的 gender 是二进制的，因此数值 by 也是构建差异平滑的好主意。但是要小心。数值 by 不会产生居中平滑。因此，s(x) + s(x, by = gender) 将隐含一个截距列，与模型截距混淆。你应该选择 s(x) + s(x, by = gender) - 1.