无法从 mgcv 中的分层 GAM 预测固定效果
Troubles predicting fixed effects from a hierarchical GAM in mgcv
我一直在使用 R 中的 mgcv 拟合不同的分层 GAM（以下简称：HGAM）。我可以毫无问题地提取和绘制它们对随机效应的预测。相反，提取和绘制他们对固定效应的预测仅适用于某些模型，我不知道为什么。
这是一个实际的例子，它指的是在不同地点采样的两个物种（Taxon）的花朵的色谱（也讨论了here）：
rm(list=ls()) # wipe R's memory clean
library(pacman) # load packages, installing them from CRAN if needed
p_load(RCurl) # allows accessing data from URL
ss <- read.delim(text=getURL("https://raw.githubusercontent.com/marcoplebani85/datasets/master/flower_color_spectra.txt"))
head(ss)
ss$density <- ifelse(ss$density<0, 0, ss$density) # set spurious negative reflectance values to zero
ss$clr <- ifelse(ss$Taxon=="SpeciesB", "red", "black")
ss <- with(ss, ss[order(Locality, wl), ])
这些是两个物种在种群水平上的平均色谱（使用滚动方法）：
每种颜色代表不同的物种。每行代表不同的地区。
根据Pedersen et al.'s classification (2019)，以下模型是 G 类型的 HGAM，它没有给出任何问题：
gam_G1 <- bam(density ~ Taxon # main effect
+ s(wl, by = Taxon, k = 20) # interaction
+ s(Locality, bs="re"), # "re" is short for "random effect"
data = ss, method = 'REML',
family="quasipoisson"
)
# gam.check(gam_G1)
# k.check(gam_G1)
# MuMIn::AICc(gam_G1)
# gratia::draw(gam_G1)
# plot(gam_G1, pages=1)
# use gam_G1 to predict wl by Locality
# dataset of predictor values to estimate response values for:
nn <- unique(ss[, c("wl", "Taxon", "Locality", "clr")])
# predict:
pred <- predict(object= gam_G1, newdata=nn, type="response", se.fit=T)
nn$fit <- pred$fit
nn$se <- pred$se.fit
# use gam_G1 to predict wl by Taxon
# dataset of predictor values to estimate response values for:
nn <- unique(ss[, c("wl",
"Taxon",
"Locality",
"clr")])
nn$Locality=0 # turns random effect off
# after https://stats.stackexchange.com/q/131106/214127
# predict:
pred <- predict(object = gam_G1,
type="response",
newdata=nn,
se.fit=T)
nn$fit <- pred$fit
nn$se <- pred$se.fit
R 警告我 factor levels 0 not in original fit，但它执行任务没有问题：
左图：gam_G1 Locality 级别的预测。右图：gam_G1 对固定效应的预测。
麻烦的模型
以下模型是“GI”类型的 HGAM sensu Pedersen et al. (2019)。它在 Locality 水平产生更准确的预测，但我只能得到 NA 作为固定效应水平的预测：
# GI: models with a global smoother for all observations,
# plus group-level smoothers, the wiggliness of which is estimated individually
start_time <- Sys.time()
gam_GI1 <- bam(density ~ Taxon # main effect
+ s(wl, by = Taxon, k = 20) # interaction
+ s(wl, by = Locality, bs="tp", m=1)
# "tp" is short for "thin plate [regression spline]"
+ s(Locality, bs="re"),
family="quasipoisson",
data = ss, method = 'REML'
)
end_time <- Sys.time()
end_time - start_time # it took ~2.2 minutes on my computer
# gam.check(gam_GI1)
# k.check(gam_GI1)
# MuMIn::AICc(gam_GI1)
尝试根据 gam_GI1:
绘制固定效应（Taxon 和 wl）的预测
# dataset of predictor values to estimate response values for:
nn <- unique(ss[, c("wl",
"Taxon",
"Locality",
"clr")])
nn$Locality=0 # turns random effect off
# after https://stats.stackexchange.com/q/131106/214127
# predict:
pred <- predict(object = gam_GI1,
type="response",
# exclude="c(Locality)",
# # this should turn random effect off
# # (doesn't work for me)
newdata=nn,
se.fit=T)
nn$fit <- pred$fit
nn$se <- pred$se.fit
head(nn)
# wl Taxon Locality clr fit se
# 1 298.34 SpeciesB 0 red NA NA
# 2 305.82 SpeciesB 0 red NA NA
# 3 313.27 SpeciesB 0 red NA NA
# 4 320.72 SpeciesB 0 red NA NA
# 5 328.15 SpeciesB 0 red NA NA
# 6 335.57 SpeciesB 0 red NA NA
左图：gam_GI1 Locality 级别的预测。右图（空白）：gam_GI1 对固定效应的预测。
以下模型，其中包括所有观察的全局平滑器，以及组级平滑器，都具有相同的“波动性”，也不提供固定效应预测：
gam_GS1 <- bam(density ~ Taxon # main effect
+ s(wl, by = Taxon, k = 20) # interaction
+ s(wl, by = Locality, bs="fs", m=1),
# "fs" is short for "factor-smoother [interaction]"
family="quasipoisson",
data = ss, method = 'REML'
)
为什么 gam_GI1 和 gam_GS1 不能对它们的固定效应进行预测，我怎样才能得到它们？
模型可能需要几分钟才能 运行。为了节省时间，可以从 here as an RData file. My R scripts (which include the code for plotting the figures) are available here.
下载他们的输出
我觉得你把几件事混为一谈了； by 关闭随机效果的技巧仅适用于 bs = "re" 平滑。 Locality 是一个因素（否则你的随机效应不是随机拦截）并将其设置为 0 正在创建一个新级别（尽管它可能会创建一个 NA，因为 0 不是' t 在原始级别中。
如果你想做的是关闭任何与 Locality 有关的东西，你应该使用 exclude；但是你有错误的调用。它不起作用的原因是因为您正在创建具有单个元素 "c(Locality)" 的字符向量。一旦您意识到 c(Locality) 与您的模型中的任何内容都不相关，这就会因显而易见的原因而失败。您需要在此处提供的是由 summary() 打印的光滑名称向量 。例如，要排除平滑的 s(Locality, bs = "re")，{mgcv} 将其识别为 s(Locality)，因此您可以使用 exclude = "s(Locality)".
在您的情况下，为每个平滑输入所有 "s(wl):LocalityLevelX" 标签是乏味的。由于您只有两个分类群，因此使用免费参数 terms 会更容易，您可以在其中列出要在模型中 包含 的平滑标签。所以你可以为这些平滑做 terms = c("s(wl):TaxonSpeciesB", "s(wl):TaxonSpeciesC") 或任何 summary() 显示。
您还需要在 terms 中包含 Taxon 项，我认为需要是：
terms = c("TaxonSpeciesB", TaxonSpeciesC",
"s(wl):TaxonSpeciesB", "s(wl):TaxonSpeciesC")
如果您安装并加载我的 {gratia} 包，您可以使用 smooths(gam_GI1) 列出 {mgcv} 知道的所有平滑标签。
by 技巧是这样工作的：
gam(y ~ x + s(z) + s(id, bs = "re", by = dummy)
其中 dummy 在拟合时设置为 数值 值 1，在预测时设置为 0。因为这是一个 numeric 变量，所以你将平滑乘以 dummy，因此为什么将其设置为 0 排除了该术语。您的代码不起作用的原因是因为您确实希望每个 Locality 的 wl 单独平滑； Locality 是您 data/model 中感兴趣的实际变量，而不是我们为实现从模型中排除某个项的目的而创建的虚拟变量。
希望现在您能明白为什么 exclude 和 terms 是比 dummy 技巧更好的解决方案。
仅供参考，在 bs = "tp" 中，"tp" 并不意味着张量积平滑。这意味着薄板回归样条（TPRS）。您只能通过 te()、t2() 或 ti() 项获得张量积平滑。
我一直在使用 R 中的 mgcv 拟合不同的分层 GAM（以下简称：HGAM）。我可以毫无问题地提取和绘制它们对随机效应的预测。相反，提取和绘制他们对固定效应的预测仅适用于某些模型，我不知道为什么。
这是一个实际的例子，它指的是在不同地点采样的两个物种（Taxon）的花朵的色谱（也讨论了here）：
rm(list=ls()) # wipe R's memory clean
library(pacman) # load packages, installing them from CRAN if needed
p_load(RCurl) # allows accessing data from URL
ss <- read.delim(text=getURL("https://raw.githubusercontent.com/marcoplebani85/datasets/master/flower_color_spectra.txt"))
head(ss)
ss$density <- ifelse(ss$density<0, 0, ss$density) # set spurious negative reflectance values to zero
ss$clr <- ifelse(ss$Taxon=="SpeciesB", "red", "black")
ss <- with(ss, ss[order(Locality, wl), ])
这些是两个物种在种群水平上的平均色谱（使用滚动方法）：
每种颜色代表不同的物种。每行代表不同的地区。
根据Pedersen et al.'s classification (2019)，以下模型是 G 类型的 HGAM，它没有给出任何问题：
gam_G1 <- bam(density ~ Taxon # main effect
+ s(wl, by = Taxon, k = 20) # interaction
+ s(Locality, bs="re"), # "re" is short for "random effect"
data = ss, method = 'REML',
family="quasipoisson"
)
# gam.check(gam_G1)
# k.check(gam_G1)
# MuMIn::AICc(gam_G1)
# gratia::draw(gam_G1)
# plot(gam_G1, pages=1)
# use gam_G1 to predict wl by Locality
# dataset of predictor values to estimate response values for:
nn <- unique(ss[, c("wl", "Taxon", "Locality", "clr")])
# predict:
pred <- predict(object= gam_G1, newdata=nn, type="response", se.fit=T)
nn$fit <- pred$fit
nn$se <- pred$se.fit
# use gam_G1 to predict wl by Taxon
# dataset of predictor values to estimate response values for:
nn <- unique(ss[, c("wl",
"Taxon",
"Locality",
"clr")])
nn$Locality=0 # turns random effect off
# after https://stats.stackexchange.com/q/131106/214127
# predict:
pred <- predict(object = gam_G1,
type="response",
newdata=nn,
se.fit=T)
nn$fit <- pred$fit
nn$se <- pred$se.fit
R 警告我 factor levels 0 not in original fit，但它执行任务没有问题：
左图：gam_G1 Locality 级别的预测。右图：gam_G1 对固定效应的预测。
麻烦的模型
以下模型是“GI”类型的 HGAM sensu Pedersen et al. (2019)。它在 Locality 水平产生更准确的预测，但我只能得到 NA 作为固定效应水平的预测：
# GI: models with a global smoother for all observations,
# plus group-level smoothers, the wiggliness of which is estimated individually
start_time <- Sys.time()
gam_GI1 <- bam(density ~ Taxon # main effect
+ s(wl, by = Taxon, k = 20) # interaction
+ s(wl, by = Locality, bs="tp", m=1)
# "tp" is short for "thin plate [regression spline]"
+ s(Locality, bs="re"),
family="quasipoisson",
data = ss, method = 'REML'
)
end_time <- Sys.time()
end_time - start_time # it took ~2.2 minutes on my computer
# gam.check(gam_GI1)
# k.check(gam_GI1)
# MuMIn::AICc(gam_GI1)
尝试根据 gam_GI1:
Taxon 和 wl）的预测
# dataset of predictor values to estimate response values for:
nn <- unique(ss[, c("wl",
"Taxon",
"Locality",
"clr")])
nn$Locality=0 # turns random effect off
# after https://stats.stackexchange.com/q/131106/214127
# predict:
pred <- predict(object = gam_GI1,
type="response",
# exclude="c(Locality)",
# # this should turn random effect off
# # (doesn't work for me)
newdata=nn,
se.fit=T)
nn$fit <- pred$fit
nn$se <- pred$se.fit
head(nn)
# wl Taxon Locality clr fit se
# 1 298.34 SpeciesB 0 red NA NA
# 2 305.82 SpeciesB 0 red NA NA
# 3 313.27 SpeciesB 0 red NA NA
# 4 320.72 SpeciesB 0 red NA NA
# 5 328.15 SpeciesB 0 red NA NA
# 6 335.57 SpeciesB 0 red NA NA
左图：gam_GI1 Locality 级别的预测。右图（空白）：gam_GI1 对固定效应的预测。
以下模型，其中包括所有观察的全局平滑器，以及组级平滑器，都具有相同的“波动性”，也不提供固定效应预测：
gam_GS1 <- bam(density ~ Taxon # main effect
+ s(wl, by = Taxon, k = 20) # interaction
+ s(wl, by = Locality, bs="fs", m=1),
# "fs" is short for "factor-smoother [interaction]"
family="quasipoisson",
data = ss, method = 'REML'
)
为什么 gam_GI1 和 gam_GS1 不能对它们的固定效应进行预测，我怎样才能得到它们？
模型可能需要几分钟才能 运行。为了节省时间，可以从 here as an RData file. My R scripts (which include the code for plotting the figures) are available here.
下载他们的输出我觉得你把几件事混为一谈了； by 关闭随机效果的技巧仅适用于 bs = "re" 平滑。 Locality 是一个因素（否则你的随机效应不是随机拦截）并将其设置为 0 正在创建一个新级别（尽管它可能会创建一个 NA，因为 0 不是' t 在原始级别中。
如果你想做的是关闭任何与 Locality 有关的东西，你应该使用 exclude；但是你有错误的调用。它不起作用的原因是因为您正在创建具有单个元素 "c(Locality)" 的字符向量。一旦您意识到 c(Locality) 与您的模型中的任何内容都不相关，这就会因显而易见的原因而失败。您需要在此处提供的是由 summary() 打印的光滑名称向量 。例如，要排除平滑的 s(Locality, bs = "re")，{mgcv} 将其识别为 s(Locality)，因此您可以使用 exclude = "s(Locality)".
在您的情况下，为每个平滑输入所有 "s(wl):LocalityLevelX" 标签是乏味的。由于您只有两个分类群，因此使用免费参数 terms 会更容易，您可以在其中列出要在模型中 包含 的平滑标签。所以你可以为这些平滑做 terms = c("s(wl):TaxonSpeciesB", "s(wl):TaxonSpeciesC") 或任何 summary() 显示。
您还需要在 terms 中包含 Taxon 项，我认为需要是：
terms = c("TaxonSpeciesB", TaxonSpeciesC",
"s(wl):TaxonSpeciesB", "s(wl):TaxonSpeciesC")
如果您安装并加载我的 {gratia} 包，您可以使用 smooths(gam_GI1) 列出 {mgcv} 知道的所有平滑标签。
by 技巧是这样工作的：
gam(y ~ x + s(z) + s(id, bs = "re", by = dummy)
其中 dummy 在拟合时设置为 数值 值 1，在预测时设置为 0。因为这是一个 numeric 变量，所以你将平滑乘以 dummy，因此为什么将其设置为 0 排除了该术语。您的代码不起作用的原因是因为您确实希望每个 Locality 的 wl 单独平滑； Locality 是您 data/model 中感兴趣的实际变量，而不是我们为实现从模型中排除某个项的目的而创建的虚拟变量。
希望现在您能明白为什么 exclude 和 terms 是比 dummy 技巧更好的解决方案。
仅供参考，在 bs = "tp" 中，"tp" 并不意味着张量积平滑。这意味着薄板回归样条（TPRS）。您只能通过 te()、t2() 或 ti() 项获得张量积平滑。