拟合`gls`模型后,将`joint_tests`映射到列表 [英] Map `joint_tests` to a list after fitting a `gls` model
问题描述
我正在尝试从具有以下代码的列表中获取具有 emmeans :: joint_tests()
的3型ANOVA表.我不完全理解该错误信息.指导我的代码来自 http://pages.stat.wisc.edu/~yandell/R_for_data_sciences/curate/tidyverse.html
I am trying to get the type 3 ANOVA table with emmeans::joint_tests()
from a list with the following code. I don't fully understand the error message.
The code that tutors me came from
http://pages.stat.wisc.edu/~yandell/R_for_data_sciences/curate/tidyverse.html
library(dplyr)
library(nlme)
library(emmeans)
data("diamonds")
diamonds %>%
split(.$cut) %>%
map(~ gls(price ~ x + y + z,
weights = varIdent(form = ~ 1|color),
data = .))%>%
map(summary)
错误消息似乎表明我以某种方式保存了我的个人模型,然后应用了 joint_tests
.我不明白的是,为什么工作流适用于 summary
,而不适用于 joint_tests
.当我们拟合单个模型时,分别打印摘要表或ANOVA表的是 summary(model)
或 joint_tests(model)
.
The error message seems to suggest that I save my individual models somehow and then apply joint_tests
. What I don't understand is why the workflow works for summary
but not for joint_tests
. When we fit single models, it's summary(model)
or joint_tests(model)
that prints the summary table or the ANOVA table, respectively.
data("diamonds")
diamonds %>%
split(.$cut) %>%
map(~ gls(price ~ x + y + z,
weights = varIdent(form = ~ 1|color),
data = .))%>%
map(joint_tests)
Error in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : Perhaps a 'data' or 'params' argument is needed
使用 map(〜joint_tests)
给出了一个奇怪的列表
Using map(~ joint_tests)
gave a weird list like this
$Fair
function (object, by = NULL, show0df = FALSE, cov.reduce = range,
...)
{
if (!inherits(object, "emmGrid")) {
args = .zap.args(object = object, cov.reduce = cov.reduce,
..., omit = "submodel")
object = do.call(ref_grid, args)
}
facs = setdiff(names(object@levels), by)
do.test = function(these, facs, result, ...) {
if ((k <- length(these)) > 0) {
emm = emmeans(object, these, by = by, ...)
tst = test(contrast(emm, interaction = "consec"),
joint = TRUE, status = TRUE)
tst = cbind(ord = k, `model term` = paste(these,
collapse = ":"), tst)
result = rbind(result, tst)
last = max(match(these, facs))
}
else last = 0
if (last < (n <- length(facs)))
for (i in last + seq_len(n - last)) result = do.test(c(these,
facs[i]), facs, result, ...)
result
}
result = suppressMessages(do.test(character(0), facs, NULL,
...))
result = result[order(result[[1]]), -1, drop = FALSE]
if (!show0df)
result = result[result$df1 > 0, , drop = FALSE]
class(result) = c("summary_emm", "data.frame")
attr(result, "estName") = "F.ratio"
attr(result, "by.vars") = by
if (any(result$note != "")) {
msg = character(0)
if (any(result$note %in% c(" d", " d e")))
msg = .dep.msg
if (any(result$note %in% c(" e", " d e")))
msg = c(msg, .est.msg)
attr(result, "mesg") = msg
}
else result$note = NULL
result
}
<bytecode: 0x7ff68eb4a0a8>
<environment: namespace:emmeans>
$Good
function (object, by = NULL, show0df = FALSE, cov.reduce = range,
...)
{
这是我在没有列表的情况下进行 joint_tests
的方式.
Here is how I did joint_tests
without a list.
diamond.gls <- gls(price ~ x + y + z,
weights = varIdent(form = ~ 1|color),
data = diamonds)
joint_tests(diamond.gls)
model term df1 df2 F.ratio p.value
x 1 14311.72 4898.859 <.0001
y 1 12964.08 46.231 <.0001
z 1 8380.71 24.576 <.0001
请查看如何解决该问题.谢谢.
Please see how I can fix it. Thank you.
推荐答案
出于我将要调查的原因, joint_tests()
是 data
参数时需要> gls 模型,至少是在函数体内调用时.为了克服这个问题,我们需要编写一个适合模型并运行 joint_tests()
的函数.这是一个平行的例子:
For reasons I will investigate, joint_tests()
needs the data
argument when it is a gls
model, at least when called within the body of a function. To overcome this, we need to write a function that fits the model and runs joint_tests()
. Here is a parallel illustration:
mod_jt = function(dat) {
mod = gls(breaks ~ tension, data = dat)
joint_tests(mod, data = dat)
}
warpbreaks %>% split(.$wool) %>% map(mod_jt)
...,我们得到结果:
... and we get the results:
$A
model term df1 df2 F.ratio p.value
tension 2 24 7.288 0.0034
$B
model term df1 df2 F.ratio p.value
tension 2 24 4.059 0.0303
我认为您拥有的代码将与 lm
模型一起使用,至少可以与 emmeans *
I think the code you had will work with an lm
model, at least with the newest version of emmeans*
这篇关于拟合`gls`模型后,将`joint_tests`映射到列表的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!