# Regression models for log transformed data without multiplicative error

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.

## 1 Answer

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)), b = coef(fit1)))


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)), b = coef(fit1)),
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?

• The log-log relationship that I am using, is standard practice for my data. Therefore, I am looking for options within the log-log framework. Firstly, thank you for your explanation - it is very thorough! However, a model which allows for larger variance for the larger values can surely work both ways, as I have seen that it under-predicts large values using my data. This means that the sum of the predicted values is a lot lower than the actual figures. Would the other models have this same issue? – sym246 Jan 9 '17 at 10:45
• The weighted nls model actually models the variance. If the variance of your residuals is not related to the y values it would fit a coefficient close to zero for the variance model. Similarly, the gamma distribution is quite flexible and can approach a gaussian distribution. If you want to put more weight on larger values, you can of course also do that and use the additive error models. However, for your artificial example data it wouldn't appear justified. – Roland Jan 9 '17 at 11:24
• Note that the log transformation is often used for the purpose of modelling data with increasing variance in increasing values. It's a variance-stabilizing transformation. That might be a reason why the log-log model is common in this specific domain. – Roland Jan 9 '17 at 11:28