r语言 - 用emmeans::ref_grid中的两个变量的组合重建一个参考网格

Question

我的真实数据与emmeans:MOats想要传达的想法具有相似的复杂性。我以MOats为例。

library(emmeans)

MOats.lm = lm(yield ~ Block + Variety, data = MOats)
ref_grid(MOats.lm)
'emmGrid' object with variables:
Block = VI, V, III, IV, II, I
Variety = Golden Rain, Marvellous, Victory
rep.meas = multivariate response levels: 0, 0.2, 0.4, 0.6
# Silly illustration of how to use 'mult.levs' to make comb's of two factors
ref_grid(MOats.lm, mult.levs = list(T=LETTERS[1:2], U=letters[1:2]))

假设MOats.lm中的Block因子不是实验设计中常见的阻断因子，而是Oat的一个特性。

主要问题:我想从Variety和Block的组合中创建一个新的变量，称为eater，具有add_grouping语法，这样如果Variety = Golden Rain x Block = I则eater = fox，如果Variety = Golden Rain x Block = II则eater = fox，如果Variety = Marvellous x Block = II则eater = cat，以此类推，使12个组合(12只是任意的，有些动物吃更多品种，有些只吃一种)。我认为我需要使Block x Variety的虚拟变量，然后分配所需的eater。最后，我想对每个品种的食客进行对比。

eater <- factor(c("fox", "cat","mouse","frog"), levels = c("fox", "cat","frog", "mouse"))

我该如何继续?add_grouping的例子只有单因素重构。如果Block的水平不能被Variety的水平整除怎么办?例如，Block有9个等级，Variety有4个等级。https://rdrr.io/cran/emmeans/man/add_grouping.html

fiber.lm <- lm(strength ~ diameter + machine, data = fiber)
( frg <- ref_grid(fiber.lm) )
# Suppose the machines are two different brands
brands <- factor(c("FiberPro", "FiberPro", "Acme"), levels = c("FiberPro", "Acme"))
( gfrg <- add_grouping(frg, "brand", "machine", brands) )

附带问题:rep.meas = multivariate response levels: 0, 0.2, 0.4, 0.6从何而来?View(MOats)中没有这样的列。

我还没有弄清楚如何从这里https://rdrr.io/github/rvlenth/emmeans/src/R/ref-grid.R的源代码以Factor1 = Factor2*Factor3的形式构造一个新变量。如有任何线索，不胜感激。

更新:以下几行添加了新的分组变量，但删除了原来的分组变量Variety和Block。

eater <- rep(LETTERS[1:3],6)
RG_add2 <- add_grouping(RG, "eater", "BV", eater)
RG_add2
'emmGrid' object with variables:
BV = 6 G, 5 G, 3 G, 4 G, 2 G, 1 G, 6 M, 5 M, 3 M, 4 M, 2 M, 1 M, 6 V, 5 V, 3 V, 4 V, 2 V, 1 V
rep.meas = multivariate response levels: 0.0, 0.2, 0.4, 0.6
eater = A, B, C
Nesting structure:  BV %in% eater

RG_add <- add_grouping(RG, "eater", "BVlev", eater)  
Error in add_grouping(RG, "eater", "BVlev", eater) : 
Length of 'newlevs' doesn't match # levels of 'BVlev'

我不明白这个错误，因为

length(BV)
[1] 18
length(eater)
[1] 18
BV
[1] "6 G" "5 G" "3 G" "4 G" "2 G" "1 G" "6 M" "5 M" "3 M" "4 M" "2 M" "1 M"
[13] "6 V" "5 V" "3 V" "4 V" "2 V" "1 V"
BVlev
[1] "6 G" "5 G" "3 G" "4 G" "2 G" "1 G" "6 M" "5 M" "3 M" "4 M" "2 M" "1 M"
[13] "6 V" "5 V" "3 V" "4 V" "2 V" "1 V"

最后，我想做emmeans(RG_add, ~ Variety|eater)

Answer 1

add_grouping()函数目前需要一个嵌套的因子。所以你需要创造一个因子。这可以使用levels<-方法完成:

library(emmeans)
MOats.lm = lm(yield ~ Block + Variety, data = MOats)
RG = ref_grid(MOats.lm)
RG
## 'emmGrid' object with variables:
##     Block = VI, V, III, IV, II, I
##     Variety = Golden Rain, Marvellous, Victory
##     rep.meas = multivariate response levels: 0, 0.2, 0.4, 0.6
BVlev = do.call(paste, expand.grid(c(6, 5, 3, 4, 2, 1), c("G", "M", "V")))
levels(RG) = list(BV = BVlev, rep.meas = c(0, 0.2, 0.4, 0.6))
RG
## 'emmGrid' object with variables:
##     BV = 6 G, 5 G, 3 G, 4 G, 2 G, 1 G, 6 M, 5 M, 3 M, 4 M, 2 M, 1 M, 6 V, 5 V, 3 V, 4 V, 2 V, 1 V
##     rep.meas = multivariate response levels: 0.0, 0.2, 0.4, 0.6

^{由reprex包(v2.0.0)在2021-08-17创建}

现在可以继续add_grouping(RG, "eater", "VB", eaters)。eaters的长度必须为18，这样每个元素都指定了与BV的每个级别相关联的eater。

在替换级别时，需要注意保持级别列表中各因素的相对顺序。组合的因素必须是连续的。