我希望将我的全局效果s(CYR.std)和我的站点偏差s(CYR.std, fSite, bs = "fs")项一起绘制在ggploy中.我可以一次画一个,但想不出如何把这两个都包括在同一个情节中.我不确定我是否在一路上正确地指定了一些东西,或者如何将s(CYR.std)项添加到因数平滑项的顶部,以进行比较.

set.seed(12345)

library(mgcv)
library(gratia)

# Hypothetical fish counts from negative binomial distribution
df <- as.data.frame(rnbinom(1000, mu = 0.6971, size = 1))

df$year <- rep(2011:2020, each=100)
df$CYR.std <- df$year - min(df$year)
df$fCYR <- as.factor(df$year)

df$site <- seq(1, 50, 1)
df$fSite <- as.factor(df$site)

df$season <- rep(c("DRY", "WET"), each=50)
df$fSeason <- as.factor(df$season)

# Depth (continuous covariate)
df$sal <- sample(0.5:40, 1000, replace = TRUE)

names(df)[1] <- "count"

m <- bam(count ~ s(sal) + 
           s(CYR.std) + 
           fSeason + 
           s(CYR.std, by = fSeason) + 
           s(CYR.std, fSite, bs = "fs") + 
           s(fCYR, bs = "re"), 
         method = "fREML",
         discrete = TRUE,
         select = TRUE, 
         family = nb(link = "log"),
         data = df)


# Site deviations from global term
ds <- data_slice(m, CYR.std = evenly(CYR.std),
                 fSite = evenly(fSite))

fv <- fitted_values(m, data = ds, scale = "response") # density

ggplot(fv, aes(x = CYR.std, y = fitted, color = fSite)) +
  # geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  geom_line()


# Global term
ds2 <- data_slice(m, CYR.std = evenly(CYR.std))

fv2 <- fitted_values(m, data = ds2, scale = "response") # density

ggplot(fv2, aes(x = CYR.std, y = fitted)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  geom_line()

Is this the correct way to specify plotting the 'fs' term and how do I combine both plots?

推荐答案

虽然Allan的答案从绘图的Angular 来看是正确的,但OP错误地指定了如何生成预测值,如果目标是在现场特定效应之上绘制全局效应.

以下是更正后的版本

library(mgcv)
library(gratia)
library(dplyr)
library(ggplot2)

您的模型具有以下平滑效果:

gratia::smooths(m)
> gratia::smooths(m)
[1] "s(sal)"                "s(CYR.std)"            "s(CYR.std):fSeasonDRY"
[4] "s(CYR.std):fSeasonWET" "s(CYR.std,fSite)"      "s(fCYR)"

还有参数fSeason效应.

当您从GAM生成预测时,您需要提供模型中使用的所有变量的值.但由于这个模型是累加的,我们可以省略一个或多个项的影响.我们通过 Select 术语或排除呼叫predict.gam()gratia::fitted_values()中的术语来做到这一点.

如果我们看看你的前data_slice()

# Site deviations from global term
ds <- data_slice(m, CYR.std = evenly(CYR.std),
                 fSite = evenly(fSite))

我们看到已经为模型中的所有未声明变量创建了值,包括salfSeason

> ds
# A tibble: 5,000 × 5
   CYR.std fSite fSeason   sal fCYR 
     <dbl> <fct> <fct>   <dbl> <fct>
 1       0 1     DRY      19.5 2011 
 2       0 2     DRY      19.5 2011 
 3       0 3     DRY      19.5 2011 
 4       0 4     DRY      19.5 2011 
 5       0 5     DRY      19.5 2011 
 6       0 6     DRY      19.5 2011 
 7       0 7     DRY      19.5 2011 
 8       0 8     DRY      19.5 2011 
 9       0 9     DRY      19.5 2011 
10       0 10    DRY      19.5 2011 
# ℹ 4,990 more rows
# ℹ Use `print(n = ...)` to see more rows

而预测是以all个这些值为条件的.

除了s(CYR.std)s(CYR.std,fSite)之外,我们可以go 除所有因素的影响,但要让这些预测与赛季无关更困难,因为处理方式不同;截取的术语将是参考赛季DRY.要么你需要接受这一点,并指出这些都是根据赛季而定的,要么就是对两个赛季都做出预测.如果你这样做了,那么你应该包括CYR.std)的平滑效果的因子吗?

让我们假设你不在乎这些预测是关于DRY赛季的.然后,排除除两个必需的影响之外的所有内容:

fv <- fitted_values(m, data = ds, scale = "response",
  terms = c("(Intercept)", "s(CYR.std)", "s(CYR.std,fSite)"))

这将只给出两个要求的平滑和模型常数项的影响.

同样,仅就全球而言,我们需要排除某些影响:

# Global term
ds2 <- data_slice(m, CYR.std = evenly(CYR.std))
fv2 <- fitted_values(m, data = ds2, scale = "response",
  terms = c("(Intercept)", "s(CYR.std)"))

我在这里使用terms而不是exclude,因为在这种情况下, Select 要包括的词条比列出要排除的词条更容易.

现在您可以遵循Allan的示例,因为我已经修改了返回的变量的名称,所以我已经修改了它以使用开发版本的{grata}:

ggplot(fv2, aes(x = CYR.std, y = .fitted)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci), alpha = 0.2) +
  geom_line(data = fv, aes(x = CYR.std, y = .fitted, color = fSite)) +
  geom_line(linewidth = 1, linetype = 2)

这给了我们

enter image description here

现在注意全局平滑上的可信区间实际上是如何用于此平滑的(加上截距)

R相关问答推荐

使用scale_x_continuous复制ggplot 2中的离散x轴

创建重复删除的唯一数据集组合列表

任意列的欧几里得距离

如何在modelsummary中重命名统计数据?

使用data.table::fcase()而不是dplyr::case_When()时保持值

如何在R中描绘#符号?

根据现有列的名称和字符串的存在进行变异以创建多个新列

派生程序包| ;无法检索';return()';的正文

查找所有站点的最小值

根据纬度和距离连接两个数据集

手动指定从相同数据创建的叠加图的 colored颜色

创建列并对大型数据集中的特定条件进行成对比较的更高效程序

为什么我对圆周率图的蒙特卡罗估计是空的?

以不同于绘图中元素的方式对GG图图例进行排序

SHILINY中DT列的条件着色

位置_道奇在geom_point图中不躲避

以R表示的NaN值的IS.NA状态

将每晚的平均值与每晚的值进行比较,统计是否有效?

合并多个数据帧,同时将它们的名称保留为列名?

打印的.txt文件,将值显示为&Quot;Num&Quot;而不是值