How to calculate percentage difference of geometric means with emmeans?

Question

I want to estimate the estimated marginal means (EMMs) of a log-normal variable across different groups. Let’s take the dataset of pigs from emmeans package for example. I first build a linear model taking the log(conc) as response variable and source and percent as predicator variables. Then, passing the model to the function of emmeans(). The EMMs can be back-transformed to the original scale with the argument of type=”response”. In the result, except for the EMMs, I also get the contrasts between groups, with a ratio column representing the exponential(mean difference of log value). But I want to get the value of 100* exponential(mean difference of log value) -100, which is more frequently presented in articles. I know I can get this value by extracting the ratio column and doing a little transformation, but how can I get the confidence interval?

library(emmeans)
data(pigs)
fit1 <- lm(log(conc)~source+factor(percent),data=pigs)
emmeans(fit1,spec= pairwise~source,type="response",infer=T)

Answer 1

Here's a way to do it using the emmeans package. It involves creating an "identity" contrast, and using the scale and offset arguments. An example using the pigs dataset follows:

library(emmeans)

pigs.lm = lm(log(conc) ~ source * factor(percent), data = pigs)
EMM = emmeans(pigs.lm, ~ percent | source, type = "response")
PRS = pairs(EMM)

# custom contrast function:
iden.emmc = function(levs, ...) {
    M = as.data.frame(diag(length(levs)))
    names(M) = levs
    M
}

confint(contrast(regrid(PRS), "iden", scale = 100, offset = -100, 
                 estName = "pct.diff"))
## source = fish:
##  contrast pct.diff    SE df lower.CL upper.CL
##  (9/12)    -16.877  8.82 17    -35.5     1.73
##  (9/15)    -17.251  9.62 17    -37.5     3.05
##  (9/18)    -20.549  8.43 17    -38.3    -2.76
##  (12/15)    -0.449 10.57 17    -22.7    21.84
##  (12/18)    -4.418  9.07 17    -23.6    14.72
##  (15/18)    -3.986 10.19 17    -25.5    17.51
## 
## source = soy:
##  contrast pct.diff    SE df lower.CL upper.CL
##  (9/12)    -13.166  8.24 17    -30.6     4.22
##  (9/15)    -12.274  8.33 17    -29.8     5.29
##  (9/18)    -19.782 10.77 17    -42.5     2.94
##  (12/15)     1.027  9.59 17    -19.2    21.26
##  (12/18)    -7.619 12.40 17    -33.8    18.55
##  (15/18)    -8.558 12.28 17    -34.5    17.34
## 
## source = skim:
##  contrast pct.diff    SE df lower.CL upper.CL
##  (9/12)    -18.480  7.74 17    -34.8    -2.15
##  (9/15)    -29.114  7.52 17    -45.0   -13.24
##  (9/18)    -41.167  7.90 17    -57.8   -24.50
##  (12/15)   -13.045  9.23 17    -32.5     6.43
##  (12/18)   -27.830  9.69 17    -48.3    -7.39
##  (15/18)   -17.004 11.82 17    -41.9     7.93
## 
## Confidence level used: 0.95

^{Created on 2021-11-11 by the reprex package (v2.0.0)}

Note that we had to also use regrid(PRS) because the ratios are not actually stored as such in the PRS object -- just the instructions for how to summarize it.

Another approach

Another approach is to replace the log transformation with a custom one that back-transforms to the desired values:

mytran = list(
    linkfun = function(mu) log(mu/100 + 1),
    linkinv = function(eta) 100 * (exp(eta) - 1),
    mu.eta = function(eta) 100 * exp(eta),
    name = "pct.diff tran"
)

update(PRS, tran = mytran, inv.lbl = "pct.diff", infer = TRUE)
## source = fish:
##  contrast pct.diff    SE df lower.CL upper.CL t.ratio p.value
##  9 / 12    -16.877  8.82 17    -38.5    12.39  -1.742  0.3339
##  9 / 15    -17.251  9.62 17    -40.5    15.16  -1.629  0.3897
##  9 / 18    -20.549  8.43 17    -41.2     7.43  -2.168  0.1723
##  12 / 15    -0.449 10.57 17    -26.4    34.60  -0.042  1.0000
##  12 / 18    -4.418  9.07 17    -27.0    25.19  -0.476  0.9634
##  15 / 18    -3.986 10.19 17    -29.0    29.82  -0.383  0.9802
## 
## source = soy:
##  contrast pct.diff    SE df lower.CL upper.CL t.ratio p.value
##  9 / 12    -13.166  8.24 17    -33.7    13.73  -1.487  0.4662
##  9 / 15    -12.274  8.33 17    -33.0    14.90  -1.380  0.5281
##  9 / 18    -19.782 10.77 17    -45.2    17.49  -1.642  0.3829
##  12 / 15     1.027  9.59 17    -22.9    32.32   0.108  0.9995
##  12 / 18    -7.619 12.40 17    -36.9    35.30  -0.590  0.9336
##  15 / 18    -8.558 12.28 17    -37.6    33.93  -0.666  0.9082
## 
## source = skim:
##  contrast pct.diff    SE df lower.CL upper.CL t.ratio p.value
##  9 / 12    -18.480  7.74 17    -37.8     6.77  -2.152  0.1767
##  9 / 15    -29.114  7.52 17    -47.6    -4.15  -3.242  0.0225
##  9 / 18    -41.167  7.90 17    -59.8   -13.83  -3.952  0.0051
##  12 / 15   -13.045  9.23 17    -35.7    17.57  -1.317  0.5651
##  12 / 18   -27.830  9.69 17    -50.7     5.70  -2.430  0.1090
##  15 / 18   -17.004 11.82 17    -44.6    24.40  -1.309  0.5699
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 4 estimates 
## Intervals are back-transformed from the pct.diff tran scale 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the pct.diff tran scale

This produces the same estimates, but different confidence intervals. The reason is that the intervals (and tests) are still performed on the scale of the fitted model, rather than using the SEs of the ratios that we obtain after applying regrid(PRS) in the previous solution. For that reason, this latter method produces the same confidence limits as we'd get by converting the original limits to the percent-difference scale:

with(confint(PRS), 100*(lower.CL - 1))
##  [1] -38.52415 -40.53746 -41.24007 -26.37452 -27.02186 -28.99035 ...

Hence I think that this latter approach is "better."