Interpreting regression coefficients of log(y+1) transformed responses I have measurements $y_1$,...,$y_i$,...,$y_n$ taken from a set of replicates in a factorial designed experiment.
In order to use a linear regression I define my response $z_i = \log(y_i + 1)$. The log is used in order to make the normality assumption hold and the model fit, and the + 1 is used since some $y_i$'s are 0.
If the model is ${z = a + \beta X}$, then the interpretation of $\beta$ is $${\beta = {\sum_{i=1}^nz_i}/n = {\sum_{i=1}^n\log(y_i + 1)}/n}$$ so $${{\log(\prod_{i=1}^n(y_i + 1))} = \beta n},$$ which gives that the geometric mean $$\sqrt[n]{\prod_{i=1}^n(y_i + 1)} = e ^\beta.$$
My question is whether there is some back transformation or some other way to get an interpretation of the geometric mean of $y_i$'s rather than of $(y_i + 1)$'s as a function of $\beta$. The reason I'm asking is that for cases where $y_i$'s are small (close to 0) the value of $\beta$ underestimates the magnitude of the effect.
 A: An alternative to thinking in terms of a GM is semi-elasticity.
Your model for the expected value is something like $$E[\ln(y+1) \vert x]= \alpha + \beta \cdot x + \gamma \cdot z$$
Taking the derivative of that with respect to $x$, you get $$\frac{\partial E[\ln(y+1) \vert x]}{\partial x}= \frac{1}{y+1}\cdot \frac{\partial y}{\partial x} = \beta $$
This can be rewritten as $$\frac{100 \cdot \frac{\Delta y}{y+1}}{\Delta x}=100 \cdot \beta,$$
which is almost the same as the equation for semi-elasticity,
$$\epsilon = \frac{100 \cdot \frac{\Delta y}{y}}{\Delta x}.$$
Common practice is to say that when $y$ is large, the two equations are basically the same, and to interpret $\beta$ as the percentage change in $y$ associated with one additional unit of $x$.
If your X is a set of dummies rather than continuous, you can do this instead.
Having said all that, it is good to think carefully about the origin of the zeros in your model. If there are lots of them, adding .01 can often drastically change your results compared to adding 1 or .001. Since the choice of constant is usually arbitrary, doing some robustness checks with different constants or even tuning the constant can work well. Another option is GLM, like a Poisson with a heteroskedastic variance that relaxes the mean-variance equality assumption, which obviates the need for any transformation entirely. The latter is my favorite option.
If zeros come from a different process than the positives, the interpretation gets more fraught. One example is expenditures on vacation and the number of children, where a big share of HHs spend nothing. More children make it less likely that a family goes on vacation, but conditional on going on one, more children lead to higher expenditures. Here the zeros come from a different choice process than the positive values, and so a more complicated model is needed.
A: Disclaimer: I would highly advise against fitting log(y+1) in order to prevent the logarithm of negative values. This might be commonly used, but that alone does not yet make it a good practice. In many cases there are better techniques. In the regression case one can use GLM or non-linear regression.
We can estimate it with a linear approximation by looking at the derivative of
$$\frac{\partial}{\partial x} f(x) = \frac{\partial}{\partial x}\sqrt[n]{\prod_{i=1}^n(y_i + 1 + x)} $$
We have
$$\frac{\partial}{\partial x} \log f(x) =  \frac{1}{n} \sum_{i=1}^{n}\frac{1}{y_i +1 +x} $$
And
$$\frac{\partial}{\partial x} f(x) = f(x)\frac{\partial}{\partial x} \left[\log f(x)\right] = f(x)   \frac{1}{n} \sum_{i=1}^{n} \frac{1}{y_i +1 +x}$$
The linear approximation is then
$$ \sqrt[n]{\prod_{i=1}^n(y_i)} \approx \sqrt[n]{\prod_{i=1}^n(y_i+1)} \cdot \left( 1-  \overline{1/y}\right)  $$
So the difference between the geometric means of $y$ and $y+1$ relates approximately to the mean of $1/y$ and you could compute bounds of the estimate by the means of $1/y$ and $1/(y+1)$
Computational example:
set.seed(1)
n = 10
y = runif(n,2,10)

# true value 5.896058
prod(y)^(1/n)

# lowerbound 5.673379  
exp(mean(log(y+1))) * (1-mean(1/y))

# upperbound 5.911496
exp(mean(log(y+1))) * (1-mean(1/(y+1)))

