mgcv 语法中的嵌套示例设计-解网

问：

我无法将我的示例设计转换为正确的 mgcv 包语法以获得随机效果。

这是设置：我们拖网捕鱼，每个“收集”地点一次，每月在相同的 12 个地点（希望如此），连续 7 个月，每年（我认为这些是我的随机拦截）。然而，有时数据库中的站点要么被意外跳过（NA），要么被备份站点取代（这是一个附近的站点，当恶劣的条件阻碍我们在那里拖网时，可以替换预期的站点）。数据还有一个重复测量方面，其中在同一地点记录了许多鱼的长度（这只是一个协变量，但认为随机截距也用于重复测量）。如果我的目标是忠实地表现这个设计，我应该如何构建我的随机效果？像这样的东西（？

（fSite = 因子位点;fCYR = 因子日历年）

s（fSite， bs='re'） + s（fCYR，， bs='re'）

或

s（fSite， fCYR， bs='re'） <- 这是否仅在每个 fCYR 级别中存在相同的 12 个站点时才有效？

数据：

> dput(station)
structure(list(CYR = c(2009, 2009, 2009, 2009, 2009, 2009, 2010, 
2010, 2010, 2010, 2010, 2010, 2010, 2011, 2011, 2011, 2011, 2011, 
2011, 2011, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2013, 2013, 
2013, 2013, 2013, 2013, 2013, 2014, 2014, 2014, 2014, 2014, 2014, 
2015, 2015, 2015, 2015, 2015, 2015, 2016, 2016, 2016, 2016, 2016, 
2016, 2016, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2018, 
2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2019, 2019, 2019, 
2019, 2019, 2019, 2019, 2019, 2019), Month = c(6, 7, 8, 9, 10, 
11, 5, 6, 7, 8, 9, 10, 11, 5, 6, 7, 8, 9, 10, 11, 5, 6, 7, 8, 
9, 10, 11, 4, 5, 6, 7, 8, 9, 10, 5, 6, 7, 8, 9, 10, 5, 6, 7, 
8, 9, 10, 5, 6, 7, 8, 9, 10, 11, 5, 6, 7, 8, 9, 10, 11, 12, 1, 
3, 5, 6, 7, 8, 9, 10, 11, 2, 3, 5, 6, 7, 8, 9, 10, 11), `Coll. Site 1` = c(21, 
23, 22, 23, 22, 23, 22, 22, 21, 68, 23, 22, 22, 70, 20, 23, 22, 
20, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 
22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 21, 22, 22, 22, 
22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 
22, 22, NA, NA, 22, 22, 22, NA, NA, 22, 22, 22, 22, 22), `Coll. Site 2` = c(40, 
70, 70, 70, 68, 68, 68, 68, 70, 101, 112, 70, 70, 143, 70, 70, 
68, 67, 70, 67, 67, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 
70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 
70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 
70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70), 
    `Coll. Site 3` = c(70, 112, 112, 122, 123, 112, 122, 122, 
    111, 136, 134, 124, 135, 158, 122, 123, 124, 124, 112, 123, 
    123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 
    123, 123, 123, 123, 123, 123, 124, 123, 123, 123, 123, 123, 
    123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 
    123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 123, 
    123, 123, 123, NA, 123, 123, 123, 123, 123, 123, 123), `Coll. Site 4` = c(147, 
    147, 133, 134, 146, 159, 54, 134, 134, 159, 156, 145, 146, 
    167, 159, 146, 144, 146, 147, 146, 145, 146, 146, 146, 146, 
    146, 146, 146, 146, 146, 146, 146, 145, 146, 146, 146, 146, 
    146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 
    146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 
    146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 146, 
    146, 146, 146, 146, 146, 146), `Coll. Site 5` = c(176, 173, 
    168, 168, 169, 172, 168, 168, 168, 168, 168, 169, 167, 241, 
    169, 172, 146, 168, 168, 169, 167, 168, 168, 168, 168, 168, 
    168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 
    168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 
    168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 168, 
    168, 168, 168, 168, 168, 168, 168, 168, NA, 168, 168, 168, 
    168, 168, 168, 168, 168), `Coll. Site 6` = c(241, 258, 241, 
    241, 241, 241, 144, 241, 169, 241, 241, 241, 241, 266, 241, 
    241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 240, 241, 
    241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 
    241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 240, 
    241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 
    241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 
    241, 241, 241, 241), `Coll. Site 7` = c(280, 279, 266, 266, 
    253, 266, 24, 266, 241, 270, 266, 266, 266, 270, 270, 266, 
    266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 257, 266, 
    266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 
    266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 
    257, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 
    266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 266, 
    266, NA, 266), `Coll. Site 8` = c(281, 283, 270, 270, 257, 
    270, NA, 270, 266, 279, 270, 270, 270, 301, 301, 270, 270, 
    270, 270, 270, 270, NA, 270, 270, 270, 270, 265, 270, 270, 
    270, 270, 270, 270, 270, 270, NA, 270, 270, 270, 270, 270, 
    270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 266, 
    270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 
    270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 270, 
    270, 270), `Coll. Site 9` = c(314, 302, 301, 314, 301, 301, 
    NA, 301, 270, 301, 301, 301, 301, 402, 402, 301, 301, 301, 
    314, 301, 314, 312, 301, 301, 302, 301, 301, 301, 301, 301, 
    NA, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 
    301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 
    301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 
    301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 301, 
    301), `Coll. Site 10` = c(401, 314, 303, 609, 314, 402, NA, 
    314, 301, 314, 314, 314, 314, 618, 618, 314, 314, 314, 403, 
    402, 402, 314, 314, 314, 314, 314, 314, 314, 314, 314, 314, 
    314, 314, 314, 314, 314, 314, 314, 315, 314, 314, 314, 314, 
    314, 314, 314, 314, 314, 314, 314, 314, 314, 314, 314, 314, 
    314, 314, 314, 314, 314, 314, 314, 314, NA, 314, 314, 314, 
    314, 314, 314, 314, 314, 314, 314, 314, 314, 314, 314, 314
    ), `Coll. Site 11` = c(618, 402, 402, 618, 402, 618, NA, 
    402, 314, 402, 402, 402, 402, 314, 266, 402, 402, 402, 609, 
    618, 618, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 
    402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 
    402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 
    402, NA, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 
    402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402, 402
    ), `Coll. Site 12` = c(NA, 621, 618, NA, 621, NA, NA, 618, 
    402, 618, 619, 618, 618, NA, NA, 621, 618, 618, 616, 314, 
    301, 618, 621, 618, 618, 618, 616, 618, 618, 618, 618, 618, 
    618, 618, 618, 618, 618, 618, 618, 618, 618, 618, 618, 618, 
    618, 618, 618, 618, 618, 618, 618, 618, 618, 618, 618, 618, 
    618, 618, 618, 618, 618, 618, 621, 618, 618, 618, 618, 618, 
    618, 618, 618, 618, 618, 618, 618, 618, 618, 618, 618)), row.names = c(NA, 
-79L), class = c("tbl_df", "tbl", "data.frame"))

r mgcv 纵向随机效应

0赞 Nate 5/23/2023

部分解决方案;我想我两者都需要。单项随机截距捕获来自同一站点和/或年份的数据点之间的相似性，而交互项为每个唯一站点|年份块的组合创建随机截距，并表示重复的度量位。来源：stackoverflow.com/questions/63542402/...

1赞 Gavin Simpson 5/24/2023

如果您要发布数据，如果它们已经采用您的模型所需的格式（或提供代码以以这种方式获取它们）会有所帮助;这些是宽格式的，需要长（整齐）的形式进行建模（即它们甚至不包含您在整个问题中提到的 Site 变量）。Stack Overflow 的一般规则是“帮助我们帮助您”。如果你不提供工作代码，或者不清理数据（或清理数据的代码），许多用户就会不看就走过去，因为这太麻烦了，甚至无法让他们开始看到问题所在。

0赞 Nate 5/24/2023

嗨，Gavin，我以这种格式添加数据，希望它只能作为采样设计问题的视觉辅助工具，而不是建模问题本身。

1赞 Gavin Simpson 5/24/2023

但是，不熟悉{mgcv}但熟悉混合效应的人可能会从您的问题中继续前进，因为要回答这个问题，他们需要查看从您的数据中构建的模型矩阵，如果不整理您的数据并试图弄清楚您使用的调用是什么，他们将无法做到这一点。如果你包含了适当的结构化数据和模型代码，你很可能已经得到了答案（对于涉及 R 的问题添加 r 标签总是有帮助的）gam()gam()

答：

1赞 Gavin Simpson 5/24/2023 #1

这是否仅在每个级别中存在相同的 12 个站点时才有效？fCYR

试试看。TLDR：是的，它有效;不，它不需要每年出现相同的 12 个站点。关键是，独特的网站应该被唯一地编码。但是，使用您在模型中提到的两种随机效应形式可能是有益的，因为它可以导致更精简的拟合。

假设您的数据符合我认为您的模型的格式并正确格式化df

library("mgcv")
library("tidyr")
library("dplyr")
library("janitor")

df2 <- pivot_longer(df, cols = c(-Month, -CYR), names_to = "site",
             values_to = "value", names_prefix = "^Coll\\. ") |>
janitor::clean_names() |>
mutate(fcyr = factor(cyr), site = factor(site))

我们可以拟合三个模型来比较随机效应的作用：

knots <- list(month = c(0.5, 12.5))
m1 <- gam(value ~ s(month, bs = "cc") +
            s(fcyr, bs = "re") + s(site, bs = "re"),
          data = df2, family = nb(), knots = knots)
m2 <- gam(value ~ s(month, bs = "cc") + s(fcyr, site, bs = "re"),
          data = df2, family = nb(), knots = knots)
m3 <- gam(value ~ s(month, bs = "cc") +
            s(fcyr, bs = "re") + s(site, bs = "re") +
            s(fcyr, site, bs = "re"),
          data = df2, family = nb(), knots = knots)

由于有 132 个独特的（站点、年份）对

> with(df2, nlevels(interaction(fcyr, site, drop = TRUE)))                    
[1] 132

m2使用近 131 个 EDF（由于模型截距 IIRC，一个 df 丢失）

> model_edf(m1)                                                               
# A tibble: 1 × 2
  model   edf
  <chr> <dbl>
1 m1     19.8

r$> model_edf(m2)                                                               
# A tibble: 1 × 2
  model   edf
  <chr> <dbl>
1 m2     132.

r$> model_edf(m3)                                                               
# A tibble: 1 × 2
  model   edf
  <chr> <dbl>
1 m3     58.1

而其中包含所有三个 ranef 术语，则在 ~58 EDF 上使用。m3

正如您的评论所推测的那样，如果您包括站点和年份的随机效应以及它们的组合，则与仅按照 .site:fcyrm2

两者都允许您估计所有 132 个（站点，年份）对的随机截距。那么有什么区别呢？m2m3

正如您的评论所推测的那样，两个“主要”随机效应编码了每个效应的平均效应和每个的平均效应。然后，该项估计单个（站点，年份）对与站点和年份平均值的偏差。一旦您估计了平均值和效应，由于年份导致的单个站点周围的偏差（或有关单个站点均值的等效年份偏差）相对较小。因此，该术语仅使用 ~40 EDF;大部分潜在偏差被缩小了，因为随机和效应已经在模型中考虑在内。siteyears(site, fcyr, bs = "re")siteyears(site,fcyr, bs = "re")siteyear

因此，在这种情况下，您可以使用任何一种形式：

~ s(fcyr, bs = "re") + s(site, bs = "re") + s(fcyr, site, bs = "re")或
~ s(fcyr, site, bs = "re")

但选项 1 更为简洁。然而，情况不一定如此。如果站点内和年内变化很大，则分解形式（1.）不会比2没有优势，这并不是不可想象的。

随机效果的主要目的是确保组被唯一地标记。假设您有 3 个字段，并从每个字段中的 6 个位置重复进行测量。您不会在每个字段中将六个位置编码为，， ...，因为模型会认为任何等于 to 的行都是相同的位置。假设字段被标记为，并且，您可以将位置编码为、、 ...、、 ...，以避免任何歧义。在 R 的混合效果软件世界中，这种问题经常出现;在那里，通过公式接口的方式，例如和工作，如果您说嵌套在 .因此，关于交叉或嵌套随机效应的讨论。如果你正确地对分组变量进行编码，整个讨论就会消失。loc1loc2loc6locationloc1ABClocA1locA2locB1locC6lme()lmer()loc1locationfield

使用 mgcv，如果您使用表格，然后拟合公式 1。模型的形式，你会把事情弄错loc1

y ~ s(field, bs="re") + s(location, bs="re") + s(field, location, bs="re")

因为该术语会认为所有标记的观察都是相同的位置（“主题”），而事实并非如此。在这种情况下，2.模型的形式将是正确的s(location, bs="re")loc1

y ~ s(field, location, bs="re")

总而言之，您使用哪种形式实际上取决于您如何编码分组因子，而不是设计是嵌套还是交叉。

上一个：ggplot2：在 aes（）中将 1 个级别的分组变量拆分为多个几何项？

下一个：无法找到签名“BAM”的函数“draw”的继承方法

mgcv 语法中的嵌套示例设计

Nested sample design in mgcv syntax

评论