我有一个GAM模型(用mgcv
拟合),它在by
参数中有两个因子的交互作用
gam(y ~ s(percent, by = source:treatment), data = dat)
如果我试图从那个模型中提取emmeans,它会给我一个错误
Error in `[.data.frame`(tbl, , vars, drop = FALSE) :
undefined columns selected
这里有一个完全可复制的例子
library(tidyverse)
library(mgcv)
library(emmeans)
# Simulate Data -----------------------------------------------------------
# similar to ?emmeans::pigs
mycurve <- function(x, h = 4, s = 2){
30-h + (2*h)/(1+exp((9-x)/s))
}
d_source <- tibble(source = factor(c("fish", "soy", "skim")),
hi = 4:6)
d_treatment <- tibble(treatment = factor(c("control", "treatment")),
si = c(2, 0.4))
dat <- expand_grid(d_source, d_treatment,
percent = seq(9, 18, by = 0.5))
set.seed(1)
dat <- dat %>%
rowwise() %>%
mutate(mu = mycurve(x = percent, h = hi, s = si),
y = mu + rnorm(n = 1, mean = 0, sd = 0.1))
dat %>%
ggplot(aes(x = percent, y = y, color = source, shape = treatment,
group = interaction(treatment, source))) +
geom_point() +
geom_line()
# Fit model ---------------------------------------------------------------
fit <- gam(y ~ s(percent, by = source:treatment), data = dat)
plot(fit, pages = 1)
# Extract EMM -------------------------------------------------------------
emmeans(fit, specs = "percent", by = c("treatment", "source"),
at = list(percent = 9:18))
你知道为什么这不起作用吗?或者是emmeans::recover_data.gam()
函数中的一个bug?
一个简单的解决方法是定义一个新的因素,但这不便于以后进行具体的对比。
dat <- dat %>%
mutate(st = interaction(source, treatment))
fit2 <- gam(y ~ s(percent, by = st), data = dat)
emmeans(fit2, specs = "percent", by = "st",
at = list(percent = 9:18))
1条答案
按热度按时间c9qzyr3d1#
我不确定,但如果你使用
emmeans::ref_grid(fit, at = list(percent = 9:18))
,它会显示从你拟合的模型中获得的参考网格的摘要,包括你可能合法使用的变量的名称。也许gam
已经为source:treatment
创建了一些变量,你可以将其用作by
变量。顺便说一句,所有的tibles和管道的东西都让你很难看到你在做什么。你发现的问题是关于
emmeans()
的,如果你保留它的其余部分,我会很感激。