Making a matrix from lsmeans contrasts return

To create the data frame:

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)

I get this return from:

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

I am looking for a simple way to turn this output (just p values) into:

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

I have about 20 sets of these that I need to turn into tables, so the simpler and more general the better.

Bonus points if the output is tab-deliminated, etc, so that I can easily paste into word/excel.

Thanks!

Solution

Here's a function that works...

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
}

Illustration:

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

Notes

As provided, this does not support by variables. Accordingly, the first line of the function disables them.
Using pairs(..., reverse = TRUE) generates the P values in the correct order needed later for upper.tri()
you can pass arguments to test() via ...

To create a tab-delimited version, use the clipr package:

clipr::write_clip(pm)

What you need is now in the clipboard and ready to paste into a spreadsheet.

Addendum

Answering this question inspired me to add a new function pwpm() to the emmeans package. It will appear in the next CRAN release, and is available now from the github site. It displays means and differences as well as P values; but the user may select which to include.

> 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