A predictor that "becomes" categorical when larger than a cutoff

Question

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)) 

Answer 1

Having two features like those described by you is a perfectly normal and valid solution. It also arises in many different scenarios. For example, if you have two features: age and gender, and consider interaction age * gender, the values for it would be 0 for people whose gender was coded as 0 and non-negative otherwise (see also the link mentioned in the comment by EdM).

I cannot comment on what the statistician said because you didn't give us a full quote and context. Maybe they meant some specific scenario, but in general, such features and their interactions are commonly used.

Answer 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.

Answer 3

Other answers are right but I think we should take in account that by including or not including the binary variable we are adjusting different models with underlying different assumptions.

Dropping the binary variable and only using the continuous variable means assuming that the effects of smoking are continuous at 0. That is, the effect of smoking very little (approaching zero) approaches the effect of not smoking.

You can see that in your example this assumption is false, because the expected response for a non smoker is 20 but the expected response for a smoker of 0 cigarettes is 25. Therefore for this example the model can match better the underlying problem when both predictors are included.

From actual data we might not know how the response behaves near zero. For some problems the knowledge on the problem may give some clue: for example, some contaminants are known not to have any safe dose, so we can suppose that very little exposure is going to produce a different response than no exposure at all and using both predictors may be useful. Other phenomena can behave differently.

In case of no prior knowledge, you can test if the categorical predictor is significant, just as with any other predictor. Here we are working in the opposite way, and finding reasons to discard it as not significant may suggest that the response is continuous at zero smoking and that may be an interesting result in itself.

Additionally, if your actual problem involves more variables, by discarding the categorical predictor you are assuming that the effects of those other variables are the same for smokers and non smokers, but by including the categorical predictor and its interactions with other variables you are actually adjusting two different models for smokers and non smokers. If you have an enough large dataset to adjust such a model without overfitting, it can give interesting results.

Answer 4

I would agree with your statistician. There is a common question, with no perfect answer. Ultimately, it comes down to what your question is. The two most common approaches I see in addiction research are (1) define 'control' as drug-exposed non-users and then define a level of drug-use above-which a user is defined as a 'case' (non-exposed non-users and exposed low-users are removed). Then use these two groups for your analysis; or (2) exclude non-exposed non-users and treat level of use as a continuous variable (drug-exposed non-users are simply '0' on this scale).

If you check the correlation between your variables:

df$NoSmoke2 = (df$NoSmoke == "NoSmoke")*1
cor(df[,-1])

You'll see that NoSmoke and Cig are very highly correlated (with set.seed(1234) r = -0.73). This makes it quite hard to interpret the regression coefficients, or the significance of either variable. For example, let's say you want to know the effect of smoking on the response variable. What's the effect? Compare your effect of Non-smoking and your effect of cigarettes when they are entered together or separately.

The issue of sample restrictions has to do with trying to move beyond pure correlations. If you want to say something about the effect of smoking, then you should only be examining people who had an opportunity to smoke. Otherwise, reverse effects (Response leads to people being more likely to smoke) is a very real possibility.

Answer 5

I would have one variable such as $\tilde x=\max(x,0)$. This type of a variable is used in a linear spline in some packages such as Stata, where it can be created with mkspline function.

Normally, these are used to model varying slopes, e.g. you may have a model where the response to negative predictor is different from response to the positive one. In this case you create two variables: $x_+=\max(0,x),x_-=\min(0,x)$.