2

我目前正在尝试使用包中的协变量（年龄）和性别 0-1 变量之间的交互进行x建模。在为每个性别指定了一个平滑项的主要模型（让我们称之为）之后，我想测试一个更简单的假设，即性别之间的差异是线性的（而不是任意平滑的）。我有以下两个问题：gammmgcvM0

  1. 在尝试正确嵌套模型时，我想从中获取 0-gender smoothing 的平滑参数M0，并将其用于更简单的模型规范中。但我收到以下错误消息：

gamm.setup 中的错误（gp，pterms = pTerms，数据 = mf，节 = 节，parametric.only = FALSE，：

gamm 无法处理链接的平滑参数（可能来自使用“id”或自适应平滑）

  1. 第二个问题更愚蠢。当我从每个性别的平滑变为 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)

gamm.setup 中的错误（gp，pterms = pTerms，数据 = mf，节 = 节，parametric.only = FALSE，：

gamm 无法处理链接的平滑参数（可能来自使用“id”或自适应平滑）

有什么想法吗？提前致谢。

4

1 回答 1

3

关于gamm错误

这是一件非常有趣的事情……好吧，我应该先解释一下逻辑。

原则上，在 中固定平滑参数是非法的gamm，因为gamm会将平滑的摆动分量视为随机效应，其方差由lme（因为您有高斯数据）来估计。如果您尝试修复平滑参数，则本质上是说您要修复随机效应的方差。好吧，lme不允许你这样做（我怀疑这种尝试在混合建模中是否合法。）

因此gamm将禁用对平滑参数的任何可能约束，包括：

  1. 平滑参数的下限，通过min.sp；
  2. 链接平滑共享相同的平滑参数，通过idin s()；
  3. 直接指定平滑参数 via sp, via s()。

前两个已完美检查。gamm没有像;这样的min.sp论点 gam即使您通过 指定它...，也没有机会使用它（因为稍后它被NULL传递给gam.setupduring gamm.setup，因此您的指定min.sp被忽略）。id您看到的错误消息也会检测到规范，但您当然没有指定id，因此上述错误未在此处报告正确的问题，因此是一个错误。

第三个，其实并没有直接通过gamm。理想情况下，一旦(by )解释了gamm/gam公式，则应将其重置为（如果不是），并发出有关此的警告消息。然而，这部分是缺失的。因此，目前您能做的最好的事情就是不指定.interpret.gamsp-1sp

你有有效的模型嵌套吗？

现在让我们转向您对嵌套的关注。嵌套与基础设置有关，与平滑参数的选择无关。只要您具有相同的一组基础（相同的基础类型、相同的数量和/或“结”的位置），模型矩阵将是相同的。例如，在您的模型M0和中，您具有与defaultM1相同的配置。因此，您的两个模型中的设计矩阵相同。只是将此复制到您的主要模型中的所有级别。也许它不容易看到，但实际上你的嵌套在.s()mgcvbs = 'tp', k = 10s()by = factor(gender)s()genderM0M1M0

让我们考虑一个小例子来验证这一点。为简单起见，我不会使用s(x)from mgcv，而是使用poly(x, degree = 2)（想象它是s(x)）。让我们首先生成一些玩具数据：

x <- 1:10
f <- gl(2, 5, labels = c("M", "F"))

由于f不是有序因子，因此通过对 的所有级别进行s(x, by = factor(f))复制来生成设计矩阵：s(x)f

## 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：

  1. As X1 + X2 = X0，设计矩阵s(x)和s(x, by = f)具有相同的列跨度，因此s(x)嵌套在 中s(x, by = f)；
  2. 因为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。

一些后续

s(x, by = as.factor(gender))我忘了在我的示例中包含我首先比较了与s(x) + s(x, by = gender)AIC参数化的相同模型（召回gender是 0-1 数值变量）。这些模型在统计上是等效的，但是在不同的情况下，平滑参数的估计显然不同，因此 AIC 略有不同。

哦是的。你gender是二进制的，所以数字by也是构建差分平滑的好主意。但请小心执行此操作。数值by不会产生居中平滑。因此，s(x) + s(x, by = gender)将隐含一个截距列，与模型截距相混淆。你应该去s(x) + s(x, by = gender) - 1。

于 2016-11-09T12:41:40.150 回答