如何在 mgcv 中使用 s(y1, by=y2:y3)) 指定具有两个分类和连续预测变量的分层 GAM (HGAM) 模型？
How to specify a hierarchical GAM (HGAM) model with two categorical & a continuous predictor using s(y1, by=y2:y3)) in mgcv?
从问题的虚拟数据开始：
zone <- c(rep("Z1",1000),rep("Z2",1000),rep("Z3",1000),rep("Z4",1000))
scenario <- rep(c(rep("S1",250),rep("S2",250),rep("S3",250),rep("S4",250)),4)
model <- rep(rep(c(rep("m1",50),rep("m2",50),rep("m3",50),rep("m4",50),rep("m5",50)),4),4)
time <- rep(seq(2021,2070), 80)
response <- rep(c(rep(rnbinom(50, mu =exp(4), size = 50), 5),
rep(rnbinom(50, mu =exp(4.5), size = 50), 5),
rep(rnbinom(50, mu =exp(5), size = 50), 5),
rep(rnbinom(50, mu =exp(6), size = 50), 5)),4)
df <- cbind(zone, scenario, model, time, response)
df <- as.data.frame(df)
df$zone <- as.factor(df$zone)
df$scenario <- as.factor(df$scenario)
df$model <- as.factor(df$model)
df$time <- as.numeric(df$time)
df$response <- as.numeric(df$response)
我正在尝试在 R 中的 mgcv 中拟合一个相当简单的 GAM 模型。我对响应变量如何随时间变化（2021 年至 2070 年之间的预测）感兴趣。我预计响应将是非线性的，所以我首先适应随着时间的推移更平滑：
library(mgcv)
m1 <- response ~ s(time)
在这种情况下，有四个模型（m1、m2、m3、m4），每个模型都预测在强度增加的四种不同“场景”（S1、S2、S3、S4）下响应在未来可能会如何变化).我对模型如何变化并不明确感兴趣（我希望它们如此），但我感兴趣的是场景（来自所有四个模型预测）在未来可能会如何变化。为了解释这种数据结构（即每个场景包含四个模型预测），我将“模型”作为随机效应包括在内：
m2 <- gam(response~s(time, k=-1) + scenario + s(model, bs='re'), data=df)
鉴于“响应”和“时间”之间的非线性响应可能因场景而异，我使用“by”函数为每个级别的场景安装了一个具有单独平滑器的模型：
m3 <- gam(response~s(time, by=scenario, k=-1) + scenario + s(model, bs='re'), data=df)
为了增加复杂性，预测的建模区域被分成四个独立的区域，因此模型中包含一个附加项“区域”：
m4 <- gam(response~s(time, by=scenario, k=-1) + scenario + zone + s(model, bs='re'), data=df)
我的问题是：我已经通过包含 s(time, by=scenario, k=-1) 来解释 m3 中随着时间响应的“摆动”可能会因场景而异的预期，但是。 .. 我如何扩展该模型以考虑响应不仅因场景而且因区域而异的可能性？我设想的是这样的：
m5 <- gam(response~s(time, by=scenario:zone, k=-1) + scenario + zone + s(model, bs='re'), data=df)
但这不是正确的语法。在多次阅读 PeerJ 中的 excellent 论文“Hierarchical generalized additive models in ecology: an introduction with mgcv” 之后，我正在努力理解简单 x~s( y) gam 模型以及函数的形状如何在两个不同的分组级别之间变化。
据我所知，我需要“单个通用平滑器加上具有不同摆动度的组级平滑器（模型 GI）”方法，以允许每个特定于组的平滑器（此处为“场景”和“ zone") 有自己的平滑参数，因此有自己的摆动程度，但我正在努力了解这个模型如何适应 HGAM 方法。如何调整 M5 以适应 HGAM 框架？
因子 by 变量实际上只是对输入数据中的行进行编码，这些行应该保持在一起并使其平滑。
所以你的直觉是对的，你需要 scenario 和 zone 的交互，但实际上你不能通过 s() 或 te()等
您需要做的是首先在您的数据中创建交互，然后将其传递给 by 参数：
df <- transform(df, zone_in_scenario = interaction(zone, scenario, drop = TRUE))
您可以在平滑中使用 by = zone_in_scenario。
您还需要将因子固定效应调整为 zone * scenario，因为您需要这两个因子的每个组合的组均值，这将为您提供每个平滑的组均值，现在适合。
从问题的虚拟数据开始：
zone <- c(rep("Z1",1000),rep("Z2",1000),rep("Z3",1000),rep("Z4",1000))
scenario <- rep(c(rep("S1",250),rep("S2",250),rep("S3",250),rep("S4",250)),4)
model <- rep(rep(c(rep("m1",50),rep("m2",50),rep("m3",50),rep("m4",50),rep("m5",50)),4),4)
time <- rep(seq(2021,2070), 80)
response <- rep(c(rep(rnbinom(50, mu =exp(4), size = 50), 5),
rep(rnbinom(50, mu =exp(4.5), size = 50), 5),
rep(rnbinom(50, mu =exp(5), size = 50), 5),
rep(rnbinom(50, mu =exp(6), size = 50), 5)),4)
df <- cbind(zone, scenario, model, time, response)
df <- as.data.frame(df)
df$zone <- as.factor(df$zone)
df$scenario <- as.factor(df$scenario)
df$model <- as.factor(df$model)
df$time <- as.numeric(df$time)
df$response <- as.numeric(df$response)
我正在尝试在 R 中的 mgcv 中拟合一个相当简单的 GAM 模型。我对响应变量如何随时间变化（2021 年至 2070 年之间的预测）感兴趣。我预计响应将是非线性的，所以我首先适应随着时间的推移更平滑：
library(mgcv)
m1 <- response ~ s(time)
在这种情况下，有四个模型（m1、m2、m3、m4），每个模型都预测在强度增加的四种不同“场景”（S1、S2、S3、S4）下响应在未来可能会如何变化).我对模型如何变化并不明确感兴趣（我希望它们如此），但我感兴趣的是场景（来自所有四个模型预测）在未来可能会如何变化。为了解释这种数据结构（即每个场景包含四个模型预测），我将“模型”作为随机效应包括在内：
m2 <- gam(response~s(time, k=-1) + scenario + s(model, bs='re'), data=df)
鉴于“响应”和“时间”之间的非线性响应可能因场景而异，我使用“by”函数为每个级别的场景安装了一个具有单独平滑器的模型：
m3 <- gam(response~s(time, by=scenario, k=-1) + scenario + s(model, bs='re'), data=df)
为了增加复杂性，预测的建模区域被分成四个独立的区域，因此模型中包含一个附加项“区域”：
m4 <- gam(response~s(time, by=scenario, k=-1) + scenario + zone + s(model, bs='re'), data=df)
我的问题是：我已经通过包含 s(time, by=scenario, k=-1) 来解释 m3 中随着时间响应的“摆动”可能会因场景而异的预期，但是。 .. 我如何扩展该模型以考虑响应不仅因场景而且因区域而异的可能性？我设想的是这样的：
m5 <- gam(response~s(time, by=scenario:zone, k=-1) + scenario + zone + s(model, bs='re'), data=df)
但这不是正确的语法。在多次阅读 PeerJ 中的 excellent 论文“Hierarchical generalized additive models in ecology: an introduction with mgcv” 之后，我正在努力理解简单 x~s( y) gam 模型以及函数的形状如何在两个不同的分组级别之间变化。
据我所知，我需要“单个通用平滑器加上具有不同摆动度的组级平滑器（模型 GI）”方法，以允许每个特定于组的平滑器（此处为“场景”和“ zone") 有自己的平滑参数，因此有自己的摆动程度，但我正在努力了解这个模型如何适应 HGAM 方法。如何调整 M5 以适应 HGAM 框架？
因子 by 变量实际上只是对输入数据中的行进行编码，这些行应该保持在一起并使其平滑。
所以你的直觉是对的，你需要 scenario 和 zone 的交互，但实际上你不能通过 s() 或 te()等
您需要做的是首先在您的数据中创建交互，然后将其传递给 by 参数：
df <- transform(df, zone_in_scenario = interaction(zone, scenario, drop = TRUE))
您可以在平滑中使用 by = zone_in_scenario。
您还需要将因子固定效应调整为 zone * scenario，因为您需要这两个因子的每个组合的组均值，这将为您提供每个平滑的组均值，现在适合。