Lognormally distributed outcome variable: regression results are way higher than crude data medians

Question

Sorry if this is trivial for you, but it's a "problem" that I am facing.

I have a lognormally distributed, extremely skewed, outcome variable. Thus, I report its value using the median instead of the mean. This is really important in this study since mean is much higher, depending more on extreme values. Therefore, reporting mean would not describe the real situation (readers somewhat get a wrong opinion on Y variable).

median(df$y)

7.5

However, when modelling this using log-link function (I need some adjusted analyses also)

model = glm(y ~ 1, data = df, family = gaussian(link = "log"))

Intercept = 2.513

Exponentiated Intercept = exp(2.5) = 12.3 (similar to mean of Y, not median of Y). Or in other words, I should report a value which is almost two times higher!

Basically, reporting modelling results means that I am not describing the real situation (y variable values are dependent on extremes). When reporting modelling results I reporting somewhat a different world from the reality? I can not throw out the extreme values as they can not be considered as outliers.

How to overcome such "problem"?

Answer 1

This is a commonly misunderstood property of the lognormal.

If $$ y \sim \operatorname{lognormal}(\mu, \sigma^2)$$

Then $E(y) = \exp(\mu + \sigma^2/2)$. This is the expectation of the lognormal random variable. If you want the median, you want $\exp(\mu)$. Remember, $\mu, \sigma^2$ are the parameters of $\log(y)$, not $y$. So, if you want to report the median of the random variable using glm you need to account for the extra factor of $\exp(\sigma^2/2)$.

Using glm,

    # Generate
    set.seed(0)
    N = 10000
    y = exp(rnorm(N, 0.5, 0.5))


    model = glm(y~1, family = gaussian(link = 'log'))


    mean(y)
    #> [1] 1.875689
    exp(coef(model))
    #> (Intercept) 
    #>    1.875689

    rmse = Metrics::rmse(log(y), predict(model))
    median(y)
    #> [1] 1.656802
    exp(coef(model))/exp(rmse^2/2)
    #> (Intercept) 
    #>    1.644235

Since you have no covariates, you could also just do...

mu = mean(log(y))
exp(mu)

EDIT: The bayesian approach is a little different.

library(tidyverse)
library(rstanarm)
library(tidybayes)


# Generate
set.seed(0)
N = 10000
y = exp(rnorm(N, 0.5, 0.5))
d = tibble(y)


model = stan_glm(log(y)~1, 
                 data = d, 
                 family = gaussian(), 
                 adapt_delta = 0.8,
                 prior_intercept = normal(0, 10))


model %>% 
  spread_draws(`(Intercept)`, sigma) %>% 
  rename(b0 = `(Intercept)` ) %>% 
  mutate(med = exp(b0)) %>% 
  pull(med) %>% 
  hist

There is a lot to consider about this problem. I wrote a little blog post discussing some nuance.