Errores estándar en R, paquete emmeans

Question

Estoy ajustando un modelo logit multinomial en R utilizando el multinom() del paquete nnet. Me gustaría obtener las proporciones de cada clase para los dos grupos. Si utilizo el paquete emmeans para hacerlo, obtengo los resultados que se indican a continuación. Si utilizo el método delta del paquete car, obtengo las mismas proporciones transformadas a la inversa, pero errores estándar diferentes. No puedo entender la razón de tal diferencia.

##########################################################################
#load data
rm(list=ls())
freq <- c(21, 23, 17, 24, 19, 13, 34, 23, 21, 26, 29, 40, 37, 31, 23,
           32, 32, 22, 29, 28, 28, 33, 20, 34, 42, 45, 40, 42, 16, 7, 16,
           13, 28, 27, 25, 19, 33, 40, 37, 36, 47, 18, 20, 13, 12, 13, 13, 
           12, 40, 22, 33, 30, 15, 11, 8, 10, 26, 27, 21, 22, 30, 24, 31,
           30, 28, 25, 27, 21, 16, 24, 13, 8, 29, 12, 18, 8, 8, 4, 5, 4, 8,
           14, 16, 7, 27, 20, 13, 5, 20, 22, 14, 17, 13, 12, 20, 19, 14, 10,
           8, 3, 3, 10, 9, 1, 8, 8, 7, 10, 19, 15, 19, 20, 20, 13, 18, 11,
           21, 26, 22, 16)
Group <- rep(c("P1", "P2"), each = 12, length.out = length(freq))
Class <- rep( paste("c", 1:5, sep = ""), each=24)

datasetR <- data.frame(
  Photo = rep(rep(1:24, 5), times = freq), 
  Group = rep(Group, times = freq),
  Class = rep(paste("c", 1:5, sep = ""), times = c(639, 614, 542, 335, 311)) )
rm(freq, Group)

#Fit multinomial model
library(nnet)
mmod <- multinom(Class ~ Group, data=datasetR)

#Calculate proportions
emmeans::emmeans(mmod, ~Class|Group, mode="prob")

#Use delta method to derive standard errors
library(car)
coefs <- as.vector(coef(mmod))
names(coefs) <- paste("b", 1:8, sep="")
vcovMmod <- vcov(mmod)
p1 <- deltaMethod(coefs, g="1/(1+(exp(b1)+exp(b2)+exp(b3)+exp(b4)))",  vcov.=vcovMmod)
p2 <- deltaMethod(coefs, g="exp(b1)/(1+(exp(b1)+exp(b2)+exp(b3)+exp(b4)))", vcov.=vcovMmod)
p3 <- deltaMethod(coefs, g="exp(b2)/(1+(exp(b1)+exp(b2)+exp(b3)+exp(b4)))", vcov.=vcovMmod)
p4 <- deltaMethod(coefs, g="exp(b3)/(1+(exp(b1)+exp(b2)+exp(b3)+exp(b4)))", vcov.=vcovMmod)
p5 <- deltaMethod(coefs, g="exp(b4)/(1+(exp(b1)+exp(b2)+exp(b3)+exp(b4)))", vcov.=vcovMmod)
p6 <- deltaMethod(coefs, g="exp(0)/(1+(exp(b1 +b5)+exp(b2+b6)+exp(b3+b7)+exp(b4+b8)))",  vcov.=vcovMmod)
p7 <- deltaMethod(coefs, g="exp(b1+b5)/(1+(exp(b1+b5)+exp(b2+b6)+exp(b3+b7)+exp(b4+b8)))", vcov.=vcovMmod)
p8 <- deltaMethod(coefs, g="exp(b2+b6)/(1+(exp(b1+b5)+exp(b2+b6)+exp(b3+b7)+exp(b4+b8)))", vcov.=vcovMmod)
p9 <- deltaMethod(coefs, g="exp(b3+b7)/(1+(exp(b1+b5)+exp(b2+b6)+exp(b3+b7)+exp(b4+b8)))", vcov.=vcovMmod)

p10 <- deltaMethod(coefs, g="exp(b4+b8)/(1+(exp(b1+b5)+exp(b2+b6)+exp(b3+b7)+exp(b4+b8)))", vcov.=vcovMmod)
result <- rbind(p1=p1,p2=p2,p3=p3,p4=p4,p5=p5,p6=p6,p7=p7,p8=p8,p9=p9,p10=p10)
result

Answer 1

No estoy seguro, pero observo que los SEs de emmeans son proporcionales a sqrt(prob*(1-prob)) dentro de cada grupo:

> x = summary(emmeans(mmod, ~Class|Group, mode="prob"))
> with(x, sqrt(prob*(1-prob)) / SE)
 [1] 33.15117 33.15117 33.15117 33.15117 33.15117 
 [6] 36.63332 36.63332 36.63332 36.63332 36.63332

Esto parece, intuitivamente, que podría ser correcto.

Apéndice

Voy a decir además que estoy seguro de que los SE de emmeans son los correctos. Se trata de un conjunto de datos simple que tiene 2 grupos, y proporciones de cada categoría en cada grupo. Observe estos resultados:

> sum(datasetR$Group == "P1")
[1] 1099
> sqrt(.Last.value)
[1] 33.15117

> sum(datasetR$Group == "P2")
[1] 1342
> sqrt(.Last.value)
[1] 36.63332

Compáralos con los resultados que he mostrado antes. Aquí están los cálculos manuales de la primera estimación y del SE:

> sum((datasetR$Group == "P1") & (datasetR$Class == "c1"))
[1] 290

> ( p11 = 290 / 1099 )
[1] 0.2638763

> sqrt(p11 * (1 - p11) / 1099)
[1] 0.01329464

Compara con:

> emmeans::emmeans(mmod, ~Class|Group, mode="prob")
Group = P1:
 Class       prob          SE df   lower.CL  upper.CL
 c1    0.26387376 0.013294604  8 0.23321635 0.2945312