时间戳

首页 > 解决方案 > 从 lsmeans 生成矩阵对比返回

r - 从 lsmeans 生成矩阵对比返回

问题描述

创建数据框：

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)

我得到这个回报：

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 中。

谢谢！

标签： rlsmeans

解决方案

这是一个有效的功能......

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)以稍后需要的正确顺序生成 P 值upper.tri()
  • 您可以将参数传递给test()via...

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

clipr::write_clip(pm)

您现在需要的是剪贴板，可以粘贴到电子表格中。

附录

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

> 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

推荐阅读