4
$\begingroup$

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

$\endgroup$
1
  • 1
    $\begingroup$ You're not estimating the mean of the logs with your glm call, but instead you're estimating the log of the mean. When you exponentiate it back you're back to estimating the mean. Instead, try glm(log(y) ~ 1, data = df, family = gaussian) ... or the simpler lm(log(y)~1,data=df). $\endgroup$
    – Glen_b
    Commented Oct 19, 2020 at 1:28

1 Answer 1

6
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ @EdM is may be possible OP is just computing the median/means from the model incorrectly. Perhaps OP can use the approach I've presented here and include some normal qqplots of $\log(y)$ so we can better diagnose the problem? $\endgroup$ Commented Oct 18, 2020 at 17:45
  • $\begingroup$ Thank you! Could you clarify two more things. (1) The estimate you calculated from your model is median, geometric mean or how should I call it on plots? (2) Also, I used your generated data and run similar bayesian model and its posterior samples (data frame) has two columns (b_Intercept and sigma) (ibb.co/n60fZWk). To calculate your shown estimate, how should I do the same calculation that you showed, using these two columns? exp(b_Intercept/sigma)? $\endgroup$
    – st4co4
    Commented Oct 18, 2020 at 18:04
  • $\begingroup$ @st4co4 1) It is the median. You can confirm this by looking at the lognormal wikipedia entry. 2) Why not just generate data from the posterior predictive distribution and compute the median? $\endgroup$ Commented Oct 18, 2020 at 19:07
  • $\begingroup$ Posterior is in lognormal scale. How should I transform them into its original scale? Just exp(b_Intercept) or should I integrate sigma also? $\endgroup$
    – st4co4
    Commented Oct 18, 2020 at 19:36
  • $\begingroup$ @st4co4 see my edit $\endgroup$ Commented Oct 18, 2020 at 20:11

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.