A predictor that "becomes" categorical when larger than a cutoff I have a dataset where a predictor can have 2 states - if it's in state 1, then the value is always the same. If it's in state 2, then the value is continuous and changes. An example to this can be number of cigarettes smoked per day. If the person is a non-smoker, the value is always 0. If the person is a smoker, the value can be 0 or higher. In my model, I have 2 predictor variables - NoSmoke and Cig. Until recently, I thought this was the correct approach (based on another SO post that I can no longer find). However, I was recently told by a statistician that this approach is problematic - partly because the 2 predictors are independent, and partly because they're correlated.
When I run the model and plot the predictions, things seem to be working out. Is there an issue with this modeling approach? And if so - what's the correct way to analyze this, while preserving the continuous predictor values?
toy example with simulated data:
library(plyr)
library(dplyr)
library(tidyr)
library(ggplot2)

df1 <- data.frame(NoSmoke = "NoSmoke", Cig = rep(0, 100), Response = rnorm(100, 20, 5))
df2 <- data.frame(NoSmoke = "Smoke", Cig = rpois(100, 2)) %>%
        mutate(Response = 5 + 2 * Cig + rnorm(100, 20, 5))
df <- rbind(df1, df2)

mod <- lm(Response ~ NoSmoke + Cig, data = df)

new <- data.frame(NoSmoke = c("NoSmoke", rep("Smoke", 7)), Cig = c(0, 0:6))
preds <- as.data.frame(predict(mod, newdata = new, interval = "confidence"))

new$Pred <- preds$fit
new$Lower <- preds$lwr
new$Upper <- preds$upr

ggplot(df) +
    geom_point(aes(x = Cig, y = Response, colour = NoSmoke), position = position_dodge(width = 0.2)) +
    geom_point(data = new, aes(x = Cig, y = Pred, fill = NoSmoke), size = 2, shape = 21, colour = "black", 
        position = position_dodge(width = 0.2)) +
    geom_errorbar(data = new, aes(x = Cig, ymin = Lower, ymax = Upper, colour = NoSmoke), 
        position = position_dodge(width = 0.2)) 

EDIT
I left this out of the question for simplicity, but following some of the answers/comments - my actual continuous variable is actually continuous, as opposed to the number of cigarettes in the example (apologies, I was trying for simplicity). In terms of correlation between the NoSmoke and the Cig variables - in my analysis, I standardize Cig, so that the standardized value for the Smoke=="Smoke" cases is standardized as normal, and the standardized value for the Smoke=="NoSmoke" is zero. This allows the model to be Response = Intercept + beta_cig*Cig for the smokers and Response = Intercept + beta_nosmoke for the non-smokers. The correlation between the NoSmoke and the CigS is low - ~0.4 in the example, and <0.1 in my real data.
Here's the toy example with the standardization:
df <- df %>%
        group_by(NoSmoke) %>%
        mutate(CigS = (Cig - mean(Cig)) / sd(Cig)) %>%
        ungroup() %>%
        mutate(CigS = ifelse(NoSmoke == "NoSmoke", 0, CigS)) 
                NoSmoke = factor(NoSmoke, levels = c("Smoke", "NoSmoke")))
cor(as.numeric(df$NoSmoke), df$CigS)

mod <- lm(Response ~ NoSmoke + CigS, data = df)

new <- data.frame(NoSmoke = c("NoSmoke", rep("Smoke", 7)), Cig = c(0, 0:6)) %>%
        left_join(unique(select(df, NoSmoke, Cig, CigS)))
preds <- as.data.frame(predict(mod, newdata = new, interval = "confidence"))

new$Pred <- preds$fit
new$Lower <- preds$lwr
new$Upper <- preds$upr

ggplot(df) +
    geom_point(aes(x = Cig, y = Response, colour = NoSmoke), position = position_dodge(width = 0.2)) +
    geom_point(data = new, aes(x = Cig, y = Pred, fill = NoSmoke), size = 2, shape = 21, colour = "black", 
        position = position_dodge(width = 0.2)) +
    geom_errorbar(data = new, aes(x = Cig, ymin = Lower, ymax = Upper, colour = NoSmoke), 
        position = position_dodge(width = 0.2)) 


@Tim, as usual, summarizes this well (+1): there is no problem with your "perfectly normal and valid solution." The following illustrates in a bit more detail.
On further review, I found your initial model much easier to interpret, with nonsmokers the reference group for the indicator variable and raw, non-negative and non-standardized values for Cig. That's the way that @whuber suggested in response to a very similar question. I used set.seed(20230212) before running your code. Then:
coef(mod)
#  (Intercept) NoSmokeSmoke          Cig 
#    19.920394     6.956709     1.511727 

is easy to interpret: (Intercept) is the outcome estimate for non-smokers, NoSmokeSmoke is the extra outcome estimate for smokers if they had Cig=0, and Cig is the extra outcome beyond that per Cig. No need to deal with un-de-meaning, interpreting coefficients in a standard-deviation-of-the-predictor scale, or similar complications.
The correlation between the two predictors isn't a problem here. It seems large:
cor(df$NoSmoke=="Smoke",df$Cig)
# [1] 0.6980322

as another answer notes. Yes, it inflates the variances of the individual coefficient estimates, but not by much:
car::vif(mod)
#  NoSmoke      Cig 
# 1.950264 1.950264 

A variance-inflation factor of that size isn't typically considered a problem. When you use the model for predictions there is a counterbalancing negative correlation between the coefficient estimates:
print(cov2cor(vcov(mod)),digits=3)
#              (Intercept) NoSmokeSmoke       Cig
# (Intercept)     1.00e+00       -0.506  9.41e-16
# NoSmokeSmoke   -5.06e-01        1.000 -6.98e-01
# Cig             9.41e-16       -0.698  1.00e+00

that leads to perfectly reasonable (and precise) predictions when the coefficient covariances are properly taken into account. You'll note that the regression restricted to the smokers gives the same result as the combined model:
coef(lm(Response~Cig,data=df,subset=NoSmoke=="Smoke"))
# (Intercept)         Cig 
#   26.877102    1.511727  

with an (Intercept) that's the sum of the original model's (Intercept) plus its NoSmokeSmoke coefficient:
sum(coef(mod)[1:2])
# [1] 26.8771

So your solution for this type of data does not have the problems that one might have feared.
