log-transformation bias correction "in reverse" when creating simulated dataset from regression predictions

Question

I am developing 1000 simulation datasets where I add random noise to predictions from a regression model where the y-variable has been log-transformed to account for non-normality. I have notice that my simulation datasets exhibit a bias that I think is related to adding values in log-units then back-transforming. I would like to know if my method for correcting the bias is appropriate. Below is reproducible code.

First, create a dummy dataset:

set.seed(8675309)
Y_var <- rlnorm(50, meanlog = 5, sdlog = 0.5)
X_var <- Y_var * exp(rnorm(50, 0, 0.55)) 

Next, I fit the model with log-tranformed variables and perform bias-correction to back-transform the predictions from log-units to original units:

lfit <- lm(log(Y_var) ~ log(X_var))
summ_lfit <- summary(lfit)
RSE <- summ_lfit$sigma
bc_pred <- exp(predict(lfit)) * exp(0.5 * RSE^2)
mean(bc_pred)
> [1] 171.4359

I then create a matrix of my bias-corrected predictions (each column is a copy of the predictions which will form the basis for the simulated datasets) that is log-transformed so I can add random noise. The noise is created from a normal distribution with mean equal to zero and SD equal to the residual standard error (RSE) from the regression model. Since the y-variable was log-transformed, the residual standard error is in log-units (hence the need to log-transform the bias-corrected predictions):

log_bc_pred <- replicate(1000, log(bc_pred)) ### 1000 simulated datasets
noise <- replicate(1000, rnorm(length(bc_pred), 0, (RSE)))

I now will add the noise to the my predictions, and back-transform to get 1000 simulated dataset with noise in the original units. However, note that the mean of means of the individual simulated datasets does not match the mean of the original predictions (bc_pred). The histogram is also not centered in the mean of bc_pred (red dashed line):

bc_pred_with_noise <- exp(log_bc_pred) * exp(noise)
mean(apply(bc_pred_with_noise, 2, mean))
> [1] 181.8478
hist(apply(bc_pred_with_noise, 2, mean))
abline(v = mean(bc_pred), col="red", lwd=3, lty=2)

However, if I apply the bias-correction "in reverse" the mean of the means of the simulated datasets match the mean of bc_pred and the histogram is centered in the mean of bc_pred:

bc_pred_with_noise2 <- exp(log_bc_pred) * exp(noise) / exp(0.5 * RSE^2)
mean(apply(bc_pred_with_noise2, 2, mean))
> [1] 171.0853
hist(apply(bc_pred_with_noise2, 2, mean))
abline(v = mean(bc_pred), col="red", lwd=3, lty=2)

It appears the applying the bias-correction "in reverse" now has my simulated datasets centered on the original dataset (as I would like it to be). However, I've not seen this application before. Can any tell me whether this is an appropriate bias-correction in my situation?

Answer 1

The phenomenon is related entirely to the method adding noise, and the setup involving a log-scale regression problem doesn't look directly relevant. Let $X$ be a generic distribution, corresponding to bc_pred in the original question. If 1000 versions of $X$ are created with additive gaussian noise, their means estimate the mean of $X$ itself, due to the linearity of expectation -- $E[X+N] = E[X]+E[N] = E[X]$. A small demonstration of this baseline model is

set.seed(42)
L <- 50
X <- rnorm(L, mean = 5, sd = 0.5)
RSE <- 1

X_resampled_noisily <- replicate(1000, X + rnorm(L,0,RSE))
mean(apply(X_resampled_noisily, 2, mean))
> [1] 4.979469

hist(apply(X_resampled_noisily, 2, mean))
abline(v = mean(X), col="red", lwd=3, lty=2)

In the question, the noise added to bc_pred is multiplicative. When $X$ and $N$ are uncorrelated, it holds in expectation that $E[XN] = E[X] E[N]$. But moreover, the original noise distribution applied in the question is a log-normal distribution with meanlog 0 and meansd RSE, whose actual mean is not one but instead $\exp(\mathrm{RSE}/2)$. So of course, it follows that applying the scale that is greater than 1 on average increases the mean of the toy sample:

set.seed(42)
L <- 50
X <- rnorm(L, mean = 5, sd = 0.5)
RSE <- 1

X_resampled_noisily <- replicate(1000, X * rlnorm(L,0,RSE))
mean(apply(X_resampled_noisily, 2, mean))

> [1] 8.216646

hist(apply(X_resampled_noisily, 2, mean))
abline(v = mean(X), col="red", lwd=3, lty=2)

The final "bias correction" example in the question replaces the log-normal multiplicative noise with a different log-normal distribution whose actual mean is 1, by multiplying the noise by $\exp(-\mathrm{RSE}/2)$, equivalent to shifting the meanlog by $-\mathrm{RSE}/2$. And indeed this works as it should:

set.seed(42)
L <- 50
X <- rnorm(L, mean = 5, sd = 0.5)
RSE <- 1

X_resampled_noisily <- replicate(1000, X * rlnorm(L,0,RSE) * exp(-RSE/2))
mean(apply(X_resampled_noisily, 2, mean))
> [1] 4.983648

hist(apply(X_resampled_noisily, 2, mean))
abline(v = mean(X), col="red", lwd=3, lty=2)

The question seems to be asking if this correction is related to procedures taken in the original regression problem. By the given definition of the regression problem, $Y$ is a variable that is proportional to the covariate $X$ by some log-normal noise whose distribution has a mean greater than 1. So, it stands to reason that the set of points $\{Y_i\}$ should have mean greater than the set of points $\{X_i\}$, and any toy dataset created from $\{X_i\}$ with noise properties that simulate the relationship $Y ~ X$ should also have a mean greater than the mean of $X$. If that is not what was expected, then either the regression model has different behavior than intended, or that expectation -- the motivation behind the "bias correction" procedure -- was wrong, and the greater mean of the resampled data should be accepted as an accurate simulation.