3
$\begingroup$

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)
$\endgroup$
5
  • $\begingroup$ It's a linear transformation, so just do the same thing to the confidence limits. $\endgroup$
    – Russ Lenth
    Nov 10, 2021 at 22:58
  • $\begingroup$ You should be compute the ratio of geometric means. $\endgroup$ Nov 11, 2021 at 14:49
  • 1
    $\begingroup$ @frankharrell that is what Kyle was getting, and the question is how to get CLs for the percent differences of the GMs. $\endgroup$
    – Russ Lenth
    Nov 11, 2021 at 14:58
  • 1
    $\begingroup$ The the title of the post needs editing. $\endgroup$ Nov 11, 2021 at 15:08
  • $\begingroup$ True. I suggest "how to compute percentage difference ..." $\endgroup$
    – Russ Lenth
    Nov 12, 2021 at 14:44

1 Answer 1

2
$\begingroup$

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."

$\endgroup$
5
  • $\begingroup$ BTW, you can add parens = NULL to remove the parentheses that contrast adds to the levels. $\endgroup$
    – Russ Lenth
    Nov 11, 2021 at 15:20
  • $\begingroup$ @kyle Please see the addendum to my answer. It is a better method. $\endgroup$
    – Russ Lenth
    Nov 13, 2021 at 16:36
  • $\begingroup$ I really appreciate your detailed answer. Thanks a lot! $\endgroup$
    – Kyle
    Nov 22, 2021 at 11:02
  • $\begingroup$ @RussLenth would it be correct to include a bias adjustment in the update() line of this code? $\endgroup$
    – myfatson
    Sep 4, 2022 at 1:22
  • $\begingroup$ Yes, I think that would be OK $\endgroup$
    – Russ Lenth
    Sep 5, 2022 at 15:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.