【发布时间】：2021-07-19 09:18:31
【问题描述】：
我一直在 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_G1Locality 级别的预测。右图：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_GI1Locality 级别的预测。右面板（空白）：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 下载为 RData 文件。我的 R 脚本（包括绘制图形的代码）可在here 获得。
【问题讨论】：
标签： r hierarchical-data mixed-models gam mgcv