使用权重的 GAM (mgcv) 的二项式模型

Question

我在尝试将二项式 GAM 拟合到数据时遇到了问题。可以通过两种方式对这些模型进行编码，(i) 提供一个比例作为响应变量，并将试验次数作为权重； (ii) 提供两栏，分别列出成功和失败。我有理由想要对我的数据点进行加权（与样本数量无关）。但是，我注意到如果我使用方法 (ii) 并添加权重（使用 weights 参数），我确实会得到非常奇怪的结果。此外，如果我提供相同的相对权重（但绝对值不同），我会得到非常不同的输出。当使用等效的 GLM 模型时（或者实际上，当使用 gam 包时）不会发生这种情况。如何为数据点提供一组权重？

这是一个 MRE：

library('mgcv')

# Random data.
x = 1:100
y_binom = cbind(rpois(100, 5 + x/2), rpois(100, 100))
w = sample(seq_len(100), 100, replace = TRUE)

# GAM models.
m1 = gam(y_binom ~ s(x), family = 'binomial')
m2 = gam(y_binom ~ s(x), weights = w / mean(w), family = 'binomial')
m3 = gam(y_binom ~ s(x), weights = w / sum(w), family = 'binomial')
m4 = gam(y_binom ~ s(x), weights = w * 100, family = 'binomial')

ms = list(m1, m2, m3, m4)

# Different RMSEs.
lapply(X = ms, FUN = function(x) return(sqrt(mean(x$residuals^2))))

# Different predictions, e.g.
plot(predict(m2), predict(m3))


# This does not happen with GLMs.
m1 = glm(y_binom ~ x, family = 'binomial')
m2 = glm(y_binom ~ x, weights = w / mean(w), family = 'binomial')
m3 = glm(y_binom ~ x, weights = w / sum(w), family = 'binomial')
m4 = glm(y_binom ~ x, weights = w * 100, family = 'binomial')

ms = list(m1, m2, m3, m4)

# Same RMSEs (for m2-m4).
lapply(X = ms, FUN = function(x) return(sqrt(mean(x$residuals^2))))

# Same predictions, e.g.
plot(predict(m2), predict(m3))

Answer 1

我认为您看到的差异是因为平滑有困难，而不是模型的 GLM 部分存在任何固有问题；您选择的权重会改变对数似然的大小，这会导致返回的模型略有不同。

我很快就会回到那个话题。首先，如果您只使用 gam():

安装普通或花园 GLM，“问题”就会消失

library('mgcv')

# Random data
set.seed(1)
x <- 1:100
y_binom <- cbind(rpois(100, 5 + x/2), rpois(100, 100))
w <- sample(seq_len(100), 100, replace = TRUE)

gam_m <- gam(y_binom ~ x, weights = w / mean(w), family = 'binomial')
glm_m <- glm(y_binom ~ x, weights = w / mean(w), family = 'binomial')

安装了完全相同的模型

> logLik(gam_m)
'log Lik.' -295.6122 (df=2)
> logLik(glm_m)
'log Lik.' -295.6122 (df=2)
> coef(gam_m)
(Intercept)           x 
 -2.1698127   0.0174864 
> coef(glm_m)
(Intercept)           x 
 -2.1698127   0.0174864

即使您通过使用不同的权重归一化来更改对数似然的大小，即使对数+似然不同，您也会得到相同的拟合模型：

gam_other <- gam(y_binom ~ x, weights = w / sum(w), family = 'binomial')

> logLik(gam_other)
'log Lik.' -2.956122 (df=2)
> coef(gam_other)
(Intercept)           x 
 -2.1698127   0.0174864

glm() 的行为在这方面是相同的：

> logLik(glm(y_binom ~ x, weights = w / sum(w), family = 'binomial'))
'log Lik.' -2.956122 (df=2)

# compare with logLik(gam_other)

在优化更边缘的情况下，这可能会崩溃，这就是 gam() 的情况。使用我的 gratia 包，我们可以轻松比较上面安装的两个 GAM：

# using your GAM m2 and m3 as examples
library(gratia)
comp <- compare_smooths(m2, m3)
draw(comp)

产生

请注意，默认情况下，这些图中的平滑包括与估计为线性的平滑时引入的偏差相关的校正。

如您所见，这两种配合是不同的；一种优化惩罚平滑一直回到线性函数，而另一种优化则没有受到太大惩罚。有了更多的数据，在 GLM 上拟合这个模型所涉及的额外复杂性（在 GAM 中我们必须 select 平滑参数）将被克服，我希望对数似然的变化不会有这么戏剧化的效果。

在这种情况下，关于 GAM 的一些理论开始变得有点松散，人们正在努力纠正或解决这些问题，但通常很难区分线性事物之间的区别或者在 link 函数的尺度上略微非线性。这里的真实函数在 link 函数的规模上略微非线性，但 m3 无法识别这一点，我认为部分原因是权重主导了似然计算。

使用权重的 GAM (mgcv) 的二项式模型

Binomial model with GAM (mgcv) using weights

r

gam

mgcv