# Linear Regression with Only Categorical Features: Evaluating the Model

Big Idea: This might seem a bit rambly, but there is a unified theme: how good is my model, and can I trust the predictions it's giving me?

Background: I am performing a linear regression (not logistic regression), where the response is continuous but all the features are categorical. I asked this question yesterday, but I have a different question today. A priori, I wouldn't expect such a model to perform terribly well: it doesn't make intuitive sense that only categorical features would be able adequately to predict where a continuous variable should be.

$$\mathbf{R}^2\!\!:$$ In my model, I've got a pretty low $$R^2$$ value, $$0.13.$$ Again, that's to be expected. But here's the kicker: the model is actually giving me very reasonable predictions. I'm wondering if the $$R^2$$ value (which I know is not the whole story of model evaluation - residual plot coming up) is just not a good number for evaluating linear regressions where all the features are categorical? Is there a better number for when all the features are categorical?

Residuals vs Fitted: As promised, here is a plot of the residuals against the fitted:

Would it make sense to try to reproduce this plot, while using some additional information to color the plot (by values of various categories)? And if the variance appears to be different for different values of a particular category, would that be evidence of heteroskedasticity? I have my suspicions about the small cluster of points to the right.

Also, is this residual plot telling me I have some systematic phenomenon going on that the model isn't capturing?

Other Diagnostics: The Q-Q plot is long-tailed, but I figure as long as I'm doing prediction (I just want to know the effect of switching from one value to another of the categorical variables) and not inference, I'm OK there. Leverages seem reasonable. There's no reason to suspect correlated errors (nothing is a time series or correlated in time). No doubt there are some outliers, but again, I'm not doing inference.

Summary: Like I said, the predictions the model is giving me seem reasonable, but I need more than that; I need to take this model to leadership and I need to be able to say that the model is reliable enough to mean what it's saying, with authority.

I've seen a related thread here, but it doesn't answer my question. Should I look into quantile regression, as this post mentions?

• I'm not sure why a linear model with only categorical predictors would necessarily produce a bad fit? If the within-category means are sufficiently far apart compared to within-category variances, the fit can be excellent. A more likely explanation seems that that isn't true for your data. Commented Feb 16 at 15:52
• There may be pre-packaged solutions out there, I would start by simply calculating and plotting means +- standard deviations across categories (probably not all 960 combinations at once though). Commented Feb 16 at 16:20
• Modeling a continuous outcome with only categorical predictors is what classical ANOVA does. It has a history of being useful if its assumptions are met. Is your model a simple linear model with two multi-level categorical predictors without interactions, as suggested in your other question? Then it's even simpler than a 2-way ANOVA that incorporates all interactions. Also, are your data necessarily non-negative as in the other question? Then reconsider your model, which is providing negative fitted values.
– EdM
Commented Feb 16 at 18:09
• @AdrianKeister If your outcome can not be negative, you can either use linear regression (making the assumption that the model will be a good approximation anyway. You would need to verify this), change the likelihood (e.g. use gamma regression), or take the log of the outcome model $E[\log(y) \mid X]$. Commented Feb 19 at 19:38
• @AdrianKeister explainability will be effected only in so far as the coefficients are concerned. The estimated effects will multiply rather than add so you're going to interpret each $\beta$ as the factor by which the expected value is multiplied by for a unit change in $x$. Commented Feb 19 at 20:20

Since the response is the Net Operating Income (NOI), if a strictly positive continuous variable, we should consider gamma regression that allows a specific pattern of variance increasing with the mean. If the response have negative values, we can shift all values by a constant before fitting, such as I(noi + 1000), and remove the constant from all predictions.

• See a tutorial at https://data.library.virginia.edu/getting-started-with-gamma-regression/ where glm(family = Gamma) estimates mean equation and a constant dispersion, despite that the dispersion estimate is inaccurate according to MASS::gamma.shape(glm()). See a comparison between gamma regression and linear regression https://civil.colorado.edu/%7Ebalajir/CVEN6833/lectures/GammaGLM-01.pdf.
• I recommend glmmTMB(noi ~ 1 + (x1 + x2 + x3)^3, dispformula = ~ (x1 + x2 + x3)^3, family = Gamma(link = "log")) which reports an accurate dispersion estimate if constant and allows the dispersion parameter to vary with some predictors. I found that glmmTMB::glmmTMB() actually fits log(shape) = log(1/dispersion) = b * x in the dispersion model and that predict(type = "disp") gives sigma instead of dispersion = sigma^2, neither of which is well documented. I reported the issue at https://github.com/glmmTMB/glmmTMB/issues/990 on Feb 20, 2024.
• Some useful relationship in gamma regression: Mean = shape * scale. Variance = mean * scale = shape * scale^2 = dispersion * mean^2 = mean^2/shape. Shape = 1/dispersion = mean^2 / var. Scale = 1/rate = var/mean = mean/shape = mean * dispersion. Dispersion = 1/shape = sigma()^2.

Nevertheless, OLS and residual plots provides useful information about data patterns. $$R^2$$ has nothing to do with how reasonable or reliable predictions are. To evaluate models, we need to compare different models and find better ones.

• $$R^2$$ is useful in selecting linear models. In a linear model, $$R^2$$ measures the proportion of variation in the response explained by all the included predictors, no matter whether the predictors are categorical or continuous. In many empirical econometric studies, $$R^2 = .13$$ is usually considered pretty high. For generalized linear models, the popular choice is AIC for model or variable selection. BIC prefers more parsimonious models. Comparatively, I have not found many versions of pseudo-R2 for generalized linear models any more useful.
• To remove incidental choice in model selection and evaluate predicting power of different models, we can use (repeated) cross validation within a training sample (e.g. 80% of cases). With a good sample size, we can consider a 10-fold cross validation repeated 5 times, so we estimate a model specification 50 times over different subsamples. We track the performance each time on the validation set, using metrics such as RMSE for a continuous response, which functions like sigma or standard error in a fitted lm(). There is also multimodel inference (see glmulti).
• To assess validity of selected models on future occasions, use performance metrics from the reserved test set that is never used in model fitting or cross validation.

Concerning in the residual plots are skewness (the red curve is downward), heteroscedasticity (spread of residual vary over fitted response), and outliers (a small cluster with huge response values). Even if the residual plots show undesired patterns, the linear model has the property that the point estimates of coefficients are unbiased. It states that the predicted mean income equals the expected value conditional on the predictors, assuming the model formula is correctly specified so that linearity is satisfied. Making predictions is conducting inference, however, so the residual diagnosis is relevant once p values, standard errors, and confidence intervals come into play.

• To test skewness or rather linearity, I recommend Ramsey's RESET test for functional form using the fitted-value approach, lmtest::resettest(lm(), power = 2:3, type = "fitted"). Other tests for linearity, such as Harvey-Collier and rainbow test, require ordering the observations properly. We should also use MacKinnon-White-Davidson PE test for comparing linear vs. log-linear specifications: lmtest::petest(lm(noi ~ .), lm(log(noi) ~ .), data = list(), vcov. = NULL, ...). To correct skewness when predictors are all categorical, we should add interaction terms, such as lm(noi ~ x1 * x2 * x3) or lm(noi ~ (x1 + x2 + x3)^2) for two-way interaction and lm(noi ~ (x1 + x2 + x3)^3) for three-way interaction. With six categorical variables of 2, 2, 3, 4, 4, and 5 levels each, the most complex model is 2 × 2 × 3 × 4 × 4 × 5 six-way interaction. With many categorical predictors, the number of possible interactional terms can be gigantic and forward stepwise selection step(lm(), direction = "forward") is useful.
• To test heteroscedasticity, I recommend White’s general test, implemented as the fitted-value approach of a Breusch-Pagan test, lmtest::bptest(lm(), varformula = ~ fitted + I(fitted^2) + I(fitted^3)) after saving fitted values into the data frame as fitted. The original Breusch-Pagan test detects linear forms of heteroskedasticity by predicting squared residuals with all predictors. White's test expands it by including squares and interactions of all predictors https://en.wikipedia.org/wiki/White_test that detects heteroskedasticity, skewness and kurtosis of the error distribution but can still be implemented as bptest() (e.g. https://www.statology.org/white-test-in-r/ and White's test for heteroskedasticity in R). Woodbridge suggests the fitted-value approach for White's test https://cran.r-project.org/web/packages/whitestrap/whitestrap.pdf. To correct heteroscedasticity, popular methods include either robust standard errors, a standard practice adopted in Stata, or generalized least squares that requires predictors of error variances. In R, the former is done in sandwich::vcovHC() and lmtest::coeftest(, vcov = ...), which only corrects standard errors without changing coefficient estimates; the latter is done by nlme::gls(, weights = varIdent(form = ~ 1 | x1 * x2 * x3)) after which we need residual diagnosis based on standardized residuals, resid(gls.fit, type = "pearson"). With many predictors, however, the pattern of heteroskedasticity can be complex, for which we may consider auxiliary (non)linear variance models to estimate error variances in a heteroskedastic linear regression model https://cran.r-project.org/web/packages/skedastic/skedastic.pdf. With six categorical predictors, defining a few heteroscedastic groups in gls seems feasible.
• To test outliers, use car::outlierTest(lm()). We can find influential cases by influence.measures(lm()) and plotting car::dfbetasPlots(lm()), where extreme values in X (not Y), or unusual combination level of categorical variables, affect coefficients the most. The small high-income cluster may not be a problem once appropriate interaction terms are included. If the sample has natural grouping, such as geographical units, hierarchical companies, repeated measurements, or accounting suppliers, we should consider clustered standard errors through sandwich::vcovCL() and mixed-effects modeling such as lme4::lmer() and nlme::lme() for linear models and lme4::glmer() and glmmTMB::glmmTMB() for generalized linear models.
• Still working my way through your very detailed answer - thank you! One note: to get correct RMSE values from a glm(... family=Gamma()) model you have to use residuals(mod, type="response") and the same type="response" in the predict function if you go that way. Commented Feb 22 at 22:55

I'm not 100% certain what you're asking here. It seems like you're lamenting that the $$R^2$$ is not impressive even though the predictions sound reasonable. I'll just give my two cents and can edit this answer when given more concrete details.

• Categorical predictors have no bearing on if $$R^2$$ is a good or bad metric for your use case. That determination will depend on what you're trying to accomplish, but in general I don't see $$R^2$$ as particularly useful, and I certainly don't think it speaks to reliability. Case in point: I can have a model with low $$R^2$$ that correctly identifies a causal estimand of interest form a very well designed experiment. The reliability comes from the design rather than from any single metric.

• What is concerning to me about your model is the switch in skewness of the residuals. Some groups seem to have a very long left tail, while others a very long right tail. This is something all together different than heterogeneity of variance. Can you comment on what you're modelling and tell us if you've cleaned the data in any way?

• I'm wanting to know if there's an equivalent to $R^2$ for when the features are all categorical. I'm wanting to know if the residual plot is concerning. You mention skewness: I'm modeling Net Operating Income as a function of several features we know in advance, all of which are categorical. The skewness switching is something I think I could have predicted in advance. For some kinds of customers, our NOI can be right-tailed, and for other kinds of customers, it can be left-tailed. As to data cleaning, that's already done by the time the data gets to me. Commented Feb 19 at 20:02
• @AdrianKeister to answer your questions: There is no need to search for an "equivalent" $R^2$, the regular will do fine. However, I don't think it gives you want you want (namely, a measure of reliability). Yes, the residual plot is concerning. Have you tried modelling the log of NOI instead? That could stabilize the variance and reign in the tails. Why do these customers have left/right skew when it comes to NOI? In general, can you edit your post to include more context on the problem you're attempting to solve or the question you're trying to answer with your data? Commented Feb 19 at 20:16
• Well, the problem is that this question is for work, and there's only so much I can share. At a basic level, of course, I'm simply trying to predict NOI on the basis of the six categorical variables to which I have access. Or, to put it another way, if I switch one categorical variable from value A to value B, what is the effect on the NOI? I will try the log transform and see if that changes anything. Commented Feb 19 at 21:38
• @AdrianKeister OK, granted you can not share all details. Unless you're trying to predict NOI for a new customer, I would say your goal here is "inference" as opposed to prediction. That being said, I don't think I can be more helpful than I already have been without additional details. Commented Feb 19 at 21:40
• @AdrianKeister. If you are interested in prediction as opposed to inference, why not measure the error the model makes as opposed to the in sample R squared? E.g. why not report (R)MSE or some other loss function? Cab you share any out of sample prediction errors? Commented Feb 20 at 1:01