Transformations and link functions in emmeans

emmeans package, Version 1.10.0

Contents

This vignette covers the intricacies of transformations and link functions in emmeans.

  1. Overview
  2. Re-gridding
  3. Link functions
  4. Graphing transformations and links
  5. Both a response transformation and a link
  6. Special transformations
  7. Specifying a transformation after the fact
  8. Auto-detected transformations
  9. Standardized response
  10. Faking a log transformation
    1. Faking other transformations
    2. Alternative scale
  11. Bias adjustment

Index of all vignette topics

Overview

Consider an example with the pigs dataset that is used in many of these vignettes; but this time, we use a log transformation of the response:

pigs.lm <- lm(log(conc) ~ source + factor(percent), data = pigs)

This model has two factors, source and percent (coerced to a factor), as predictors; and log-transformed conc as the response. Here we obtain the EMMs for source, examine its structure, and finally produce a summary, including a test against a null value of log(35):

emm.src <- emmeans(pigs.lm, "source")
str(emm.src)
## 'emmGrid' object with variables:
##     source = fish, soy, skim
## Transformation: "log"
summary(emm.src, infer = TRUE, null = log(35))
##  source emmean     SE df lower.CL upper.CL null t.ratio p.value
##  fish     3.39 0.0367 23     3.32     3.47 3.56  -4.385  0.0002
##  soy      3.67 0.0374 23     3.59     3.74 3.56   2.988  0.0066
##  skim     3.80 0.0394 23     3.72     3.88 3.56   6.130  <.0001
## 
## Results are averaged over the levels of: percent 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Now suppose that we want the EMMs expressed on the same scale as conc. This can be done by adding type = "response" to the summary() call:

summary(emm.src, infer = TRUE, null = log(35), type = "response")
##  source response   SE df lower.CL upper.CL null t.ratio p.value
##  fish       29.8 1.09 23     27.6     32.1   35  -4.385  0.0002
##  soy        39.1 1.47 23     36.2     42.3   35   2.988  0.0066
##  skim       44.6 1.75 23     41.1     48.3   35   6.130  <.0001
## 
## Results are averaged over the levels of: percent 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## Tests are performed on the log scale

Note: Looking ahead, this output is compared later in this vignette with a bias-adjusted version.

Timing is everything

Dealing with transformations in emmeans is somewhat complex, due to the large number of possibilities. But the key is understanding what happens, when. These results come from a sequence of steps. Here is what happens (and doesn’t happen) at each step:

  1. The reference grid is constructed for the log(conc) model. The fact that a log transformation is used is recorded, but nothing else is done with that information.
  2. The predictions on the reference grid are averaged over the four percent levels, for each source, to obtain the EMMs for sourcestill on the log(conc) scale.
  3. The standard errors and confidence intervals for these EMMs are computed – still on the log(conc) scale.
  4. Only now do we do back-transformation…
    1. The EMMs are back-transformed to the conc scale.
    2. The endpoints of the confidence intervals are back-transformed.
    3. The t tests and P values are left as-is.
    4. The standard errors are converted to the conc scale using the delta method. These SEs were not used in constructing the tests and confidence intervals.

The model is our best guide

This choice of timing is based on the idea that the model is right. In particular, the fact that the response is transformed suggests that the transformed scale is the best scale to be working with. In addition, the model specifies that the effects of source and percent are linear on the transformed scale; inasmuch as marginal averaging to obtain EMMs is a linear operation, that averaging is best done on the transformed scale. For those two good reasons, back-transforming to the response scale is delayed until the very end by default.

Back to Contents

Re-gridding

As well-advised as it is, some users may not want the default timing of things. The tool for changing when back-transformation is performed is the regrid() function – which, with default settings of its arguments, back-transforms an emmGrid object and adjusts everything in it appropriately. For example:

str(regrid(emm.src))
## 'emmGrid' object with variables:
##     source = fish, soy, skim
summary(regrid(emm.src), infer = TRUE, null = 35)
##  source response   SE df lower.CL upper.CL null t.ratio p.value
##  fish       29.8 1.09 23     27.5     32.1   35  -4.758  0.0001
##  soy        39.1 1.47 23     36.1     42.2   35   2.827  0.0096
##  skim       44.6 1.75 23     40.9     48.2   35   5.446  <.0001
## 
## Results are averaged over the levels of: percent 
## Confidence level used: 0.95

Notice that the structure no longer includes the transformation. That’s because it is no longer relevant; the reference grid is on the conc scale, and how we got there is now forgotten. Compare this summary() result with the preceding one, and note the following:

Understood, right? But think carefully about how these EMMs were obtained. They are back-transformed from emm.src, in which the marginal averaging was done on the log scale. If we want to back-transform before doing the averaging, we need to call regrid() after the reference grid is constructed but before the averaging takes place:

pigs.rg <- ref_grid(pigs.lm)
remm.src <- emmeans(regrid(pigs.rg), "source")
summary(remm.src, infer = TRUE, null = 35)
##  source response   SE df lower.CL upper.CL null t.ratio p.value
##  fish       30.0 1.10 23     27.7     32.2   35  -4.585  0.0001
##  soy        39.4 1.49 23     36.3     42.5   35   2.927  0.0076
##  skim       44.8 1.79 23     41.1     48.5   35   5.486  <.0001
## 
## Results are averaged over the levels of: percent 
## Confidence level used: 0.95

These results all differ from either of the previous two summaries – again, because the averaging is done on the conc scale rather than the log(conc) scale.

Note: For those who want to routinely back-transform before averaging, the regrid argument in ref_grid() simplifies this. The first two steps above could have been done more easily as follows:

remm.src <- emmeans(pigs.lm, "source", regrid = "response")

But don’t get regrid and type confused. The regrid argument is passed to regrid() (as the transform argument) after the reference grid is constructed, whereas the type argument is simply remembered and used by summary(). So a similar-looking call:

emmeans(pigs.lm, "source", type = "response")

will compute the results we have seen for emm.src – back-transformed after averaging on the log scale.

Remember again: When it comes to transformations, timing is everything.

Back to Contents

Graphing transformations and links

There are a few options for displaying transformed results graphically. First, the type argument works just as it does in displaying a tabular summary. Following through with the neuralgia example, let us display the marginal Treatment EMMs on both the link scale and the response scale (we are opting to do the averaging on the link scale):

neur.Trt.emm <- suppressMessages(emmeans(neuralgia.glm, "Treatment"))
plot(neur.Trt.emm)   # Link scale by default

Plot A: Display of the results of confint(neur.Trt.emm)

plot(neur.Trt.emm, type = "response")

Plot B: Display of the results of confint(neur.Trt.emm, type = "response"). These intervals are markedly skewed right|left for low|high estimated probabilities

Besides whether or not we see response values, there is a dramatic difference in the symmetry of the intervals.

For emmip() and plot() only (and currently only with the “ggplot” engine), there is also the option of specifying type = "scale", which causes the response values to be calculated but plotted on a nonlinear scale corresponding to the transformation or link:

plot(neur.Trt.emm, type = "scale")

Plot C: On the inside, this plot looks exactly like Plot A above, but the scale is transformed to show the values in Plot B. However, there are not enough tick marks.

Notice that the interior part of this plot is identical to the plot on the link scale. Only the horizontal axis is different. That is because the response values are transformed using the link function to determine the plotting positions of the graphical elements – putting them back where they started.

As is the case here, nonlinear scales can be confusing to read, and it is very often true that you will want to display more scale divisions, and even add minor ones. This is done via adding arguments for the function ggplot2::scale_x_continuous() (see its documentation):

plot(neur.Trt.emm, type = "scale", breaks = seq(0.10, 0.90, by = 0.10),
     minor_breaks = seq(0.05, 0.95, by = 0.05))

Plot D: Same as Plot C except there are more tick marks so we can discern the values better

When using the "ggplot" engine, you always have the option of using ggplot2 to incorporate a transformed scale – and it doesn’t even have to be the same as the transformation used in the model. For example, here we display the same results on an arcsin-square-root scale.

plot(neur.Trt.emm, type = "response") +
  ggplot2::scale_x_continuous(trans = scales::asn_trans(),
                              breaks = seq(0.10, 0.90, by = 0.10))

Plot E: An alternative to Plot D using an arcsin scaling. This scale is less nonlinear than in plot D so the intervals are somewhat skewed, but less so than plot B

This comes across as a compromise: not as severe as the logit scaling, and not as distorted as the linear scaling of response values.

Again, the same techniques can be used with emmip(), except it is the vertical scale that is affected.

Back to Contents

Special transformations

The make.tran() function provides several special transformations and sets things up so they can be handled in emmeans with relative ease. (See help("make.tran", "emmeans") for descriptions of what is available.) make.tran() works much like stats::make.link() in that it returns a list of functions linkfun(), linkinv(), etc. that serve in managing results on a transformed scale. The difference is that most transformations with make.tran() require additional arguments.

To use this capability in emmeans(), it is fortuitous to first obtain the make.tran() result, and then to use it as the enclosing environment for fitting the model, with linkfun as the transformation. For example, suppose the response variable is a percentage and we want to use the response transformation \(\sin^{-1}\sqrt{y/100}\). Then proceed like this:

tran <- make.tran("asin.sqrt", 100)
my.model <- with(tran, 
    lmer(linkfun(percent) ~ treatment + (1|Block), data = mydata))

Subsequent calls to ref_grid(), emmeans(), regrid(), etc. will then be able to access the transformation information correctly.

The help page for make.tran() has an example like this using a Box-Cox transformation.

Back to Contents

Specifying a transformation after the fact

It is not at all uncommon to fit a model using statements like the following:

mydata <- transform(mydata, logy.5 = log(yield + 0.5))
my.model <- lmer(logy.5 ~ treatment + (1|Block), data = mydata)

In this case, there is no way for ref_grid() to figure out that a response transformation was used. What can be done is to update the reference grid with the required information:

my.rg <- update(ref_grid(my.model), tran = make.tran("genlog", .5))

Subsequently, use my.rg in place of my.model in any emmeans() analyses, and the transformation information will be there.

For standard transformations (those in stats::make.link()), just give the name of the transformation; e.g.,

model.rg <- update(ref_grid(model), tran = "sqrt")

Auto-detected response transformations

As can be seen in the initial pigs.lm example in this vignette, certain straightforward response transformations such as log, sqrt, etc. are automatically detected when emmeans() (really, ref_grid()) is called on the model object. In fact, scaling and shifting is supported too; so the preceding example with my.model could have been done more easily by specifying the transformation directly in the model formula:

my.better.model <- lmer(log(yield + 0.5) ~ treatment + (1|Block), data = mydata)

The transformation would be auto-detected, saving you the trouble of adding it later. Similarly, a response transformation of 2 * sqrt(y + 1) would be correctly auto-detected. A model with a linearly transformed response, e.g. 4*(y - 1), would not be auto-detected, but 4*I(y + -1) would be interpreted as 4*identity(y + -1). Parsing is such that the response expression must be of the form mult * fcn(resp + const); operators of - and / are not recognized.

Back to Contents

Faking a log transformation

The regrid() function makes it possible to fake a log transformation of the response. Why would you want to do this? So that you can make comparisons using ratios instead of differences.

Consider the pigs example once again, but suppose we had fitted a model with a square-root transformation instead of a log:

pigroot.lm <- lm(sqrt(conc) ~ source + factor(percent), data = pigs)
logemm.src <- regrid(emmeans(pigroot.lm, "source"), transform = "log")
confint(logemm.src, type = "response")
##  source response   SE df lower.CL upper.CL
##  fish       29.8 1.32 23     27.2     32.7
##  soy        39.2 1.54 23     36.2     42.6
##  skim       45.0 1.74 23     41.5     48.7
## 
## Results are averaged over the levels of: percent 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
pairs(logemm.src, type = "response")
##  contrast    ratio     SE df null t.ratio p.value
##  fish / soy  0.760 0.0454 23    1  -4.591  0.0004
##  fish / skim 0.663 0.0391 23    1  -6.965  <.0001
##  soy / skim  0.872 0.0469 23    1  -2.548  0.0457
## 
## Results are averaged over the levels of: percent 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

These results are not identical, but very similar to the back-transformed confidence intervals above for the EMMs and the pairwise ratios in the “comparisons” vignette, where the fitted model actually used a log response.

Faking other transformations

It is possible to fake transformations other than the log. Just use the same method, e.g. (results not displayed)

regrid(emm, transform = "probit")

would re-grid the existing emm to the probit scale. Note that any estimates in emm outside of the interval \((0,1)\) will be flagged as non-estimable.

The section on standardized responses gives an example of reverse-engineering a standardized response transformation in this way.

Revisiting the neuralgia example

As was mentioned before in the neuralgia example, some users prefer ratios of probabilities (risk ratios) rather than odds ratios. The additional machinery of regrid() makes this possible. First, do

log.emm <- regrid(neuralgia.emm, "log")

While each node of neuralgia.emm is an estimate of the logit of a probability, the corresponding nodes of log.emm are estimates of the log of the same probability. Thus, pairwise comparisons are differences of logs, which are logs of ratios. The risk ratios are thus obtainable by

pairs(log.emm, reverse = TRUE, type = "response")
##  contrast ratio    SE  df null z.ratio p.value
##  B / A    0.572 0.484 Inf    1  -0.661  0.7863
##  P / A    4.107 2.221 Inf    1   2.613  0.0244
##  P / B    7.181 5.185 Inf    1   2.730  0.0174
## 
## Results are averaged over the levels of: Sex 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

The test statistics and P values differ somewhat from those for the odds ratios because they are computed on the log scale rather than the original logit scale.

We were able to obtain both odds ratios and risk ratios for neuralgia.glm. But what if we had not used the logit link? Then the odds ratios would not just fall out naturally. However, we can regrid() to the "logit" scale if we want odds ratios, or to "log" scale if we want risk ratios. For example,

neuralgia.prb <- glm(Pain ~ Treatment * Sex + Age, family = binomial(link = "probit"), 
                     data = neuralgia)
prb.emm <- suppressMessages(emmeans(neuralgia.prb, "Treatment"))
pairs(regrid(prb.emm, "logit"), type = "response", reverse = TRUE)
##  contrast odds.ratio     SE  df null z.ratio p.value
##  B / A         0.517  0.509 Inf    1  -0.670  0.7809
##  P / A        15.561 14.358 Inf    1   2.975  0.0082
##  P / B        30.077 33.225 Inf    1   3.081  0.0058
## 
## Results are averaged over the levels of: Sex 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log odds ratio scale

These are vaguely comparable to the odds ratios we obtained with neuralgia.glm, Similar re-gridding for the log will give us risk ratios.

Alternative scale

It is possible to create a report on an alternative scale by updating the tran component. For example, suppose we want percent differences instead of ratios in the preceding example with the pigs dataset. This is possible by modifying the reverse transformation: since the uusual reverse transformation is a ratio of the form \(r = a/b\), we have that the percentage difference between \(a\) and \(b\) is \(100(a-b)/b = 100(r-1)\). Thus,

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

update(pairs(logemm.src, type = "response"), 
       tran = pct.diff.tran, inv.lbl = "pct.diff", adjust = "none",
       infer = c(TRUE, TRUE))
##  contrast    pct.diff   SE df lower.CL upper.CL t.ratio p.value
##  fish / soy     -24.0 4.54 23    -32.8   -13.97  -4.591  0.0001
##  fish / skim    -33.7 3.91 23    -41.3   -25.08  -6.965  <.0001
##  soy / skim     -12.8 4.69 23    -22.0    -2.54  -2.548  0.0180
## 
## Results are averaged over the levels of: percent 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log(pct.diff) scale 
## Tests are performed on the log(pct.diff) scale

Another way to obtain the same estimates is to directly transform the estimated ratios to \(100r - 100\):

contrast(regrid(pairs(logemm.src)), "identity", scale = 100, offset = -100,
         infer = c(TRUE, TRUE))
##  contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
##  (fish/soy)     -24.0 4.54 23    -33.3    -14.6  -5.281  <.0001
##  (fish/skim)    -33.7 3.91 23    -41.8    -25.6  -8.613  <.0001
##  (soy/skim)     -12.8 4.69 23    -22.5     -3.1  -2.730  0.0119
## 
## Results are averaged over the levels of: percent 
## Confidence level used: 0.95

While the estimates are the same, the tests and confidence intervals are different because they are computed on the re-gridded scale using the standard errors shown, rather than on the link scale as in the first results.

Standardized response

In some disciplines, it is common to fit a model to a standardized response variable. R’s base function scale() makes this easy to do; but it is important to notice that scale(y) is more complicated than, say, sqrt(y), because scale(y) requires all the values of y in order to determine the centering and scaling parameters. The ref_grid() function (called by `emmeans() and others) tries to detect the scaling parameters. To illustrate:

fiber.lm <- lm(scale(strength) ~ machine * scale(diameter), data = fiber)
emmeans(fiber.lm, "machine")   # on the standardized scale
##  machine   emmean    SE df lower.CL upper.CL
##  A        0.00444 0.156  9   -0.349    0.358
##  B        0.28145 0.172  9   -0.109    0.672
##  C       -0.33473 0.194  9   -0.774    0.105
## 
## Results are given on the scale(40.2, 4.97) (not the response) scale. 
## Confidence level used: 0.95
emmeans(fiber.lm, "machine", type = "response")   # strength scale
##  machine response    SE df lower.CL upper.CL
##  A           40.2 0.777  9     38.5     42.0
##  B           41.6 0.858  9     39.7     43.5
##  C           38.5 0.966  9     36.3     40.7
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the scale(40.2, 4.97) scale

More interesting (and complex) is what happens with emtrends(). Without anything fancy added, we have

emtrends(fiber.lm, "machine", var = "diameter")
##  machine diameter.trend     SE df lower.CL upper.CL
##  A                0.222 0.0389  9   0.1339    0.310
##  B                0.172 0.0450  9   0.0705    0.274
##  C                0.174 0.0418  9   0.0791    0.268
## 
## Confidence level used: 0.95

These slopes are (change in scale(strength)) / (change in diameter); that is, we didn’t do anything to undo the response transformation, but the trend is based on exactly the variable specified, diameter. To get (change in strength) / (change in diameter), we need to undo the response transformation, and that is done via regrid (which invokes regrid() after the reference grid is constructed):

emtrends(fiber.lm, "machine", var = "diameter", regrid = "response")
##  machine diameter.trend    SE df lower.CL upper.CL
##  A                1.104 0.194  9    0.666     1.54
##  B                0.857 0.224  9    0.351     1.36
##  C                0.864 0.208  9    0.394     1.33
## 
## Confidence level used: 0.95

What if we want slopes for (change in scale(strength)) / (change in scale(diameter))? This can be done, but it is necessary to manually specify the scaling parameters for diameter.

with(fiber, c(mean = mean(diameter), sd = sd(diameter)))
##      mean        sd 
## 24.133333  4.323799
emtrends(fiber.lm, "machine", var = "scale(diameter, 24.133, 4.324)")
##  machine scale(diameter, 24.133, 4.324).trend    SE df lower.CL upper.CL
##  A                                      0.960 0.168  9    0.579     1.34
##  B                                      0.745 0.195  9    0.305     1.19
##  C                                      0.751 0.181  9    0.342     1.16
## 
## Confidence level used: 0.95

This result is the one most directly related to the regression coefficients:

coef(fiber.lm)[4:6]
##          scale(diameter) machineB:scale(diameter) machineC:scale(diameter) 
##                0.9598846               -0.2148202               -0.2086880

There is a fourth possibility, (change in strength) / (change in scale(diameter)), that I leave to the reader.

What to do if auto-detection fails

Auto-detection of standardized responses is a bit tricky, and doesn’t always succeed. If it fails, a message is displayed and the transformation is ignored. In cases where it doesn’t work, we need to explicitly specify the transformation using make.tran(). The methods are exactly as shown earlier in this vignette, so we show the code but not the results for a hypothetical example.

One method is to fit the model and then add the transformation information later. In this example, some.fcn is a model-fitting function which for some reason doesn’t allow the scaling information to be detected.

mod <- some.fcn(scale(RT) ~ group + (1|subject), data = mydata)
emmeans(mod, "group", type = "response",
        tran = make.tran("scale", y = mydata$RT))

The other, equivalent, method is to create the transformation object first and use it in fitting the model:

mod <- with(make.tran("scale", y = mydata$RT),
            some.fcn(linkfun(RT) ~ group + (1|subject), data = mydata))
emmeans(mod, "group", type = "response")

Reverse-engineering a standardized response

An interesting twist on all this is the reverse situation: Suppose we fitted the model without the standardized response, but we want to know what the results would be if we had standardized. Here we reverse-engineer the fiber.lm example above:

fib.lm <- lm(strength ~ machine * diameter, data = fiber)

# On raw scale:
emmeans(fib.lm, "machine")
##  machine emmean    SE df lower.CL upper.CL
##  A         40.2 0.777  9     38.5     42.0
##  B         41.6 0.858  9     39.7     43.5
##  C         38.5 0.966  9     36.3     40.7
## 
## Confidence level used: 0.95
# On standardized scale:
tran <- make.tran("scale", y = fiber$strength)
emmeans(fib.lm, "machine", regrid = tran)
##  machine   emmean    SE df lower.CL upper.CL
##  A        0.00444 0.156  9   -0.349    0.358
##  B        0.28145 0.172  9   -0.109    0.672
##  C       -0.33473 0.194  9   -0.774    0.105
## 
## Results are given on the scale(40.2, 4.97) (not the response) scale. 
## Confidence level used: 0.95

In the latter call, the regrid argument causes regrid() to be called after the reference grid is constructed.

Back to Contents

Bias adjustment

So far, we have discussed ideas related to back-transforming results as a simple way of expressing results on the same scale as the response. In particular, means obtained in this way are known as generalized means; for example, a log transformation of the response is associated with geometric means. When the goal is simply to make inferences about which means are less than which other means, and a response transformation is used, it is often acceptable to present estimates and comparisons of these generalized means. However, sometimes it is important to report results that actually do reflect expected values of the untransformed response. An example is a financial study, where the response is in some monetary unit. It may be convenient to use a response transformation for modeling purposes, but ultimately we may want to make financial projections in those same units.

In such settings, we need to make a bias adjustment when we back-transform, because any nonlinear transformation biases the expected values of statistical quantities. More specifically, suppose that we have a response \(Y\) and the transformed response is \(U\). To back-transform, we use \(Y = h(U)\); and using a Taylor approximation, \(Y \approx h(\eta) + h'(\eta)(U-\eta) + \frac12h''(\eta)(U-\eta)^2\), so that \(E(Y) \approx h(\eta) + \frac12h''(\eta)Var(U)\). This shows that the amount of needed bias adjustment is approximately \(\frac12h''(\eta)\sigma^2\) where \(\sigma\) is the error SD in the model for \(U\). It depends on \(\sigma\), and the larger this is, the greater the bias adjustment is needed. This second-order bias adjustment is what is currently used in the emmeans package when bias-adjustment is requested. There are better or exact adjustments for certain cases, and future updates may incorporate some of those.

Pigs example revisited

Let us compare the estimates in the overview after we apply a bias adjustment. First, note that an estimate of the residual SD is available via the sigma() function:

sigma(pigs.lm)
## [1] 0.115128

This estimate is used by default. The bias-adjusted EMMs for the sources are:

summary(emm.src, type = "response", bias.adj = TRUE)
##  source response   SE df lower.CL upper.CL
##  fish       30.0 1.10 23     27.8     32.4
##  soy        39.4 1.48 23     36.5     42.6
##  skim       44.9 1.77 23     41.3     48.7
## 
## Results are averaged over the levels of: percent 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## Bias adjustment applied based on sigma = 0.11513

These estimates (and also their SEs) are slightly larger than we had without bias adjustment. They are estimates of the arithmetic mean responses, rather than the geometric means shown in the overview. Had the value of sigma been larger, the adjustment would have been greater. You can experiment with this by adding a sigma = argument to the above call.