GAM 的年度计数指数按站点查看长期趋势

Annual count index from GAM looking at long-term trends by site

我有兴趣使用广义加性模型 (gams) 来估计在几个不同站点监测的计数随时间的共享全球趋势。我读过 Pederson et al. (2019) 对分层游戏 (hgams) 的精彩介绍,我相信我可以按如下方式设置模型(Pederson 等人 (2019) GS 模型),

fit_model = gam(count ~ s(year, m = 2) + s(year, site, bs = 'fs', m = 2),
                data = count_df,
                family = nb(link = 'log'),
                method = 'REML')

我可以绘制局部效应平滑图,查看拟合诊断,一切看起来都很合理。我的问题是如何提取非中心的年度相对计数指数?我的第一个想法是将估计截距(时间序列开始时跨站点的平均计数)添加到 s(year) 平滑(共享全局平滑)。但是我不确定围绕平滑的不确定性是否已经包含了估计截距中的不确定性?或者如果我需要添加它?多亏了令人惊叹的 R 库 mgcvgratiadplyr.

,所有这一切才成为可能

你的方法不包括常数项中的不确定性,它只是改变了周围的一切。

如果你想这样做,使用 constant 参数到 gratia:::draw.gam() 会更容易:

draw(fit_model, select = "s(year)", constant = coef(fit_model)[1L])

它可以完成您的代码所做的工作,而无需(您)付出太多努力。

一个更好的方法——使用 {gratia},因为你已经在使用它了——创建一个数据框,其中包含一系列超过 year 的值,然后使用 gratia::fitted_values() 从模型中为 year 的这些值生成估计值。为了得到你想要的东西(这似乎是为了排除拟合的随机平滑分量,这样你就将随机分量设置为等于 link 尺度上的 0)你需要将该平滑传递给 exclude 参数:

## data to predict at
new_year <- with(count_df,
                 tibble(year = gratia::seq_min_max(year, n = 100),
                        site = factor(levels(site)[1], levels = levels(site)))

## predict
fv <- fitted_values(fit_model, data = new_year, exclude = "s(year,site)")

如果您想了解 exclude,请参阅 ?predict.gam