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

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"? • 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). Oct 19, 2020 at 1:28 ## 1 Answer 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.875689 exp(coef(model)) #> (Intercept) #> 1.875689 rmse = Metrics::rmse(log(y), predict(model)) median(y) #>  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. • @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? Oct 18, 2020 at 17:45
• 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)? Oct 18, 2020 at 18:04
• @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? Oct 18, 2020 at 19:07
• Posterior is in lognormal scale. How should I transform them into its original scale? Just exp(b_Intercept) or should I integrate sigma also? Oct 18, 2020 at 19:36
• @st4co4 see my edit Oct 18, 2020 at 20:11