Создание матрицы из контраста lsmeans 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)

Я получаю это возвращение от:

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 наборов из них, которые мне нужно превратить в таблицы, поэтому, чем проще и общеизвестнее, тем лучше. et c, так что я могу легко вставить в слово / 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) создает значения P в правильном порядке, необходимом позже для upper.tri()
• , вы можете передавать аргументы в test() через ...

Чтобы создать версию с разделителями табуляции, используйте пакет 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