问题描述
我一直在使用 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,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[,"clr")])
nn$Locality=0 # turns random effect off
# after https://stats.stackexchange.com/q/131106/214127
# predict:
pred <- predict(object = gam_G1,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,k = 20) # interaction
+ s(wl,by = Locality,bs="tp",m=1)
# "tp" is short for "thin plate [regression spline]"
+ s(Locality,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)
尝试根据 Taxon 为固定效应（wl 和 gam_GI1）绘制预测：
# dataset of predictor values to estimate response values for:
nn <- unique(ss[,"clr")])
nn$Locality=0 # turns random effect off
# after https://stats.stackexchange.com/q/131106/214127
# predict:
pred <- predict(object = gam_GI1,# 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,bs="fs",m=1),# "fs" is short for "factor-smoother [interaction]"
family="quasipoisson",method = 'REML'
)
为什么 gam_GI1 和 gam_GS1 不对其固定效应产生预测，我如何获得它们？
模型可能需要几分钟才能运行。为了节省时间，他们的输出可以从 here 下载为 RData 文件。我的 R 脚本（包括绘制图形的代码）可用 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() 显示。
您还需要在 Taxon 中包含 terms 术语，我认为应该是：
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。由于这是一个 数字 变量，您将平滑乘以 dummy，因此为什么将其设置为 0 排除了该术语。您的代码不起作用的原因是因为您确实希望为每个 wl 分别对 Locality 进行平滑处理； Locality 是您的数据/模型中感兴趣的实际变量，而不是我们为实现从模型中排除项的目的而创建的虚拟变量。
希望现在您能明白为什么 exclude 和 terms 是比 dummy 技巧更好的解决方案。
仅供参考，在 bs = "tp" 中，"tp" 并不意味着张量积平滑。这意味着薄板回归样条（TPRS）。您只能通过 te()、t2() 或 ti() 项获得张量积平滑。