How to handle a relationship for only non-zero values of dependent variables in linear regression

Question

I have a dependent variable in my model which is non-zero only for part of all observations. The example dataset could look like this:

library(dplyr)
library(ggplot2)
df <- data.frame(x = c(rep(0, 10), 
                       runif(12, 0, 2000), 
                       rep(0, 8), 
                       runif(14, 0, 2000), 
                       rep(0, 23), 
                       runif(33, 0, 2000))) %>% 
  mutate(y = exp(x / 1000) + rnorm(100, 0, 1)) %>% 
  mutate(y = ifelse(x == 0, runif(100, 0, 10), y))

df %>% 
  ggplot() + 
  geom_point(aes(x, y))

So as you can see, there is a visible relationship for non-zero values of x and y, but if x == 0 the values of response variable are still varied. It's quite obvious that the problem occurs due to omitted variable bias and although I have other dependent variables to feed into my model, they can't really capture all the points of interest. I know that one way would be to limit my dataset to only those observations for which x > 0, however it's unfortunetely not an option for me. So I wanted to ask how can I handle such a relationship in multivariable linear regression so I could fit a non-linear function to x (probably exponens in the case shown)? I was thinking about something like creating a dummy variable that would indicate if x > 0 and potentially capture a little bit of variance from x == 0 by predicting the mean of all those points. However, I'm not really sure if that's a valid approach and was thinking about possible alternatives.

Answer 1

You can probably fit this with a model that looks like this:

$$E[y \vert x]= \alpha+\beta \cdot I(x>0) + \gamma \cdot x \cdot I(x>0)+ \delta \cdot x^2 \cdot I(x>0), $$

where $I()=1$ if $x>0$ and $0$ otherwise.

There is probably a more elegant way to do this in R, but this should do it:

library(dplyr)
library(ggplot2)
df <- data.frame(x = c(rep(0, 10), 
                       runif(12, 0, 2000), 
                       rep(0, 8), 
                       runif(14, 0, 2000), 
                       rep(0, 23), 
                       runif(33, 0, 2000))) %>% 
  mutate(y = exp(x / 1000) + rnorm(100, 0, 1)) %>% 
  mutate(y = ifelse(x == 0, runif(100, 0, 10), y))

df$pos_x <- ifelse(df$x >0 , 1, 0)
ols<-lm(y~pos_x + pos_x:I(x) + pos_x:I(x^2),data=df)

df %>% 
  ggplot() + 
  geom_point(aes(x, y)) +    
  geom_point(aes(x,predict(ols), color='red'))

The predictions from the model are shown in red:

Answer 2

This is a common problem in many data-sets. Look into zero-inflated models and hurdle models. Depending on the specifics of your variables, a version of one of those should be what you are looking for.