Regression models for log transformed data without multiplicative error

Question

I have a set of data with one explanatory and one response variable. They are both extremely positively skewed, and so have been transformed using a log to make them 'more normal'.

When I created a linear regression between the two variables, the fit was very good (R squared of 0.85), but the larger values were very highly under-predicted once back-transformed, due to the multiplicative nature of errors using a log transformation.

The following example shows what i mean:

set.seed(10)

x=rlnorm(100,5,1)
y=rlnorm(100,2,2)

x=sort(x, decreasing = FALSE)
y=sort(y, decreasing = FALSE)

DF=data.frame(x=x,y=y)

## Plot relationship between variables
plot(log(y)~log(x))

## Create regression using logged data
fit=lm(log(y)~log(x), data=DF)
summary(fit)

## Plot regression line
plot(log(y)~log(x))
abline(-7.936712,1.990450, col="red")

## Compute predicted y values by back-transforming
DF$Predicted=(exp(-7.936712)*(DF$x^1.990450))

## Calculate sum of actual vs. predicted.
sum(DF$y)
# 4632.657
sum(DF$Predicted)
# 3792.603

## Create model between actual and predicted.
pred_fit=lm(Predicted~y-1, data=DF)
summary(pred_fit)
plot(Predicted~y-1,data=DF)
abline(0,1, col="red")

I have been advised to try other models (such as GLMs), but can't seem to work out exactly how these are applicable. My reason for this is:

• The relationship between the variables seems to be linear once the log transformation of both the response and explanatory variables have been applied. Therefore, a GLM would be subjected to the Gaussian family (correct me if i'm wrong), and so there is no difference to what i have already.

If I apply a GLM to the un-transformed data, using a log-link function, then will this apply the log transformation to my response or explanatory variable (or both), and would I need to back-transform afterward, like i do with the linear model?

Additionally, I don't see if this would solve the multiplicative error problem, which is my motivation for exploring this. Finally, I would like to view the results of this GLM on a plot using the log-log scale, so i can see how well the model fits the data. Not sure if this would be possible, but it would probably help me to understand.

Answer 1

I illustrate five options to fit a model here. The assumption for all of them is that the relationship is actually $y = a \cdot x^b$ and we only need to decide on the appropriate error structure.

1.) First the OLS model $\ln{y} = a + b\cdot\ln{x}+\varepsilon$, i.e., a multiplicative error after back-transformation.

fit1 <- lm(log(y) ~ log(x), data = DF)

I would argue that this is actually an appropriate error model as you clearly have increasing scatter with increasing values.

2.) A non-linear model $y = \alpha\cdot x^b+\varepsilon$, i.e., an additive error.

fit2 <- nls(y ~ a * x^b, data = DF, start = list(a = exp(coef(fit1)[1]), b = coef(fit1)[2]))

3.) A Generalized Linear Model with Gaussian distribution and a log link function. We will see that this is actually the same model as 2 when we plot the result.

fit3 <- glm(y ~ log(x), data = DF, family = gaussian(link = "log"))

4.) A non-linear model as 2, but with a variance function $s^2(y) = \exp(2\cdot t \cdot y)$, which adds an additional parameter.

library(nlme)
fit4 <- gnls(y ~ a * x^b, params = list(a ~ 1, b ~ 1),
             data = DF, start = list(a = exp(coef(fit1)[1]), b = coef(fit1)[2]), 
             weights = varExp(form = ~ y))

5.) A GLM with a gamma distribution and a log link.

fit5 <- glm(y ~ log(x), data = DF, family = Gamma(link = "log"))

Now let's plot them:

plot(y ~ x, data = DF)
curve(exp(predict(fit1, newdata = data.frame(x = x))), col = "green", add = TRUE)
curve(predict(fit2, newdata = data.frame(x = x)), col = "black", add = TRUE)
curve(predict(fit3, newdata = data.frame(x = x), type = "response"), col = "red", add = TRUE, lty = 2)
curve(predict(fit4, newdata = data.frame(x = x)), col = "brown", add = TRUE)
curve(predict(fit5, newdata = data.frame(x = x), type = "response"), col = "cyan", add = TRUE)

legend("topleft", legend = c("OLS", "nls", "Gauss GLM", "weighted nls", "Gamma GLM"),
       col = c("green", "black", "red", "brown", "cyan"),
       lty = c(1, 1, 2, 1, 1))

I hope these fits persuade you that you actually should use a model that allows larger variance for larger values. Even the model where I fit a variance model agrees on that. If you use the non-linear model or Gaussian GLM you place undue weight on the larger values.

Finally, you should consider carefully if the assumed relationship is the correct one. Is it supported by scientific theory?