从 lsmeans contrasts return 制作矩阵

Question

创建数据框：

num <- sample(1:25, 20)
x <- data.frame("Day_eclosion" = num, "Developmental" = c("AP", "MA", 
"JU", "L"), "Replicate" = 1:5)

model <- glmer(Day_eclosion ~ Developmental +  (1 | Replicate), family = 
"poisson", data= x)

我从 return 那里得到这个 return：

a <- lsmeans(model, pairwise~Developmental, adjust = "tukey")
a$contrasts

contrast estimate     SE  df z.ratio p.value
 AP - JU    0.2051 0.0168 Inf 12.172  <.0001 
 AP - L     0.3009 0.0212 Inf 14.164  <.0001 
 AP - MA    0.3889 0.0209 Inf 18.631  <.0001 
 JU - L     0.0958 0.0182 Inf  5.265  <.0001 
 JU - MA    0.1839 0.0177 Inf 10.387  <.0001 
 L - MA     0.0881 0.0222 Inf  3.964  0.0004

我正在寻找一种简单的方法将此输出（仅 p 值）转换为：

     AP     MA     JU    L
AP   -   <.0001  <.0001 <.0001 
MA   -       -   <.0001 0.0004  
JU   -       -      -   <.0001
L    -       -            -

我有大约 20 套这些需要转换成表格，所以越简单、越通用越好。

如果输出以制表符分隔等方式加分，这样我就可以轻松粘贴到 word/excel。

谢谢！

Answer 1

这是一个有效的函数...

pvmat = function(emm, ...) {
    emm = update(emm, by = NULL)    # need to work harder otherwise
    pv = test(pairs(emm, reverse = TRUE, ...)) $ p.value
    fmtpv = sprintf("%6.4f", pv) 
    fmtpv[pv < 0.0001] = "<.0001"
    lbls = do.call(paste, emm@grid[emm@misc$pri.vars])
    n = length(lbls)
    mat = matrix("", nrow = n, ncol = n, dimnames = list(lbls, lbls))
    mat[upper.tri(mat)] = fmtpv
    idx = seq_len(n - 1)
    mat[idx, 1 + idx]   # trim off last row and 1st col
}

插图：

require(emmeans)

> warp.lm = lm(breaks ~ wool * tension, data = warpbreaks)
> warp.emm = emmeans(warp.lm, ~ wool * tension)

> warp.emm
 wool tension emmean   SE df lower.CL upper.CL
 A    L         44.6 3.65 48     37.2     51.9
 B    L         28.2 3.65 48     20.9     35.6
 A    M         24.0 3.65 48     16.7     31.3
 B    M         28.8 3.65 48     21.4     36.1
 A    H         24.6 3.65 48     17.2     31.9
 B    H         18.8 3.65 48     11.4     26.1

Confidence level used: 0.95 

> pm = pvmat(warp.emm, adjust = "none")
> print(pm, quote=FALSE)
    B L    A M    B M    A H    B H   
A L 0.0027 0.0002 0.0036 0.0003 <.0001
B L        0.4170 0.9147 0.4805 0.0733
A M               0.3589 0.9147 0.3163
B M                      0.4170 0.0584
A H                             0.2682

备注

根据规定，这不支持 by 变量。因此，函数的第一行禁用它们。
使用 pairs(..., reverse = TRUE) 以正确的顺序生成稍后 upper.tri()
您可以通过 ...

test()

要创建制表符分隔版本，请使用 clipr 包：

clipr::write_clip(pm)

您需要的内容现在已在剪贴板中，随时可以粘贴到电子表格中。

附录

回答这个问题启发我在 emmeans 包中添加了一个新函数 pwpm()。它将出现在下一个 CRAN 版本中，现在可以从 the github site 获得。它显示平均值和差异以及 P 值；但用户可能 select 要包含的内容。

> pwpm(warp.emm)

wool = A
       L      M      H
L [44.6] 0.0007 0.0009
M 20.556 [24.0] 0.9936
H 20.000 -0.556 [24.6]

wool = B
       L      M      H
L [28.2] 0.9936 0.1704
M -0.556 [28.8] 0.1389
H  9.444 10.000 [18.8]

Row and column labels: tension
Upper triangle: P values   adjust = “tukey”
Diagonal: [Estimates] (emmean) 
Upper triangle: Comparisons (estimate)   earlier vs. later

从 lsmeans contrasts return 制作矩阵

Making a matrix from lsmeans contrasts return

r

lsmeans

附录