# Missing Coefficients in Linear Regression with Multiple Categorical Variables in R

I have an odd scenario where I am trying to regress a numerical variable on several categorical variables, with no other numerical variables. I have roughly 23k rows of data in my real-world example. The R function lm handles categoricals just fine - at least mostly. However, I find that I don't get coefficients for every value of the categorical variables. Here's a minimum working example:

charges = c(1000, 2000, 3000, 4000, 5000)
geo = c("Local", "Regional", "National", "Local", "Regional")
region = c("Texas", "Oklahoma", "Louisiana", "Texas", "Louisiana")

rev = data.frame(charges, geo, region)

rev$$geo = factor(rev$$geo)
rev$$region = factor(rev$$region)

mod = lm(charges~., data=rev)
summary(mod)


The result is:


Call:
lm(formula = charges ~ ., data = rev)

Residuals:
1          2          3          4          5
-1.500e+03  1.484e-13  9.154e-14  1.500e+03  1.484e-13

Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept)        2500       1500   1.667    0.344
geoNational         500       2598   0.192    0.879
geoRegional        2500       2598   0.962    0.512
regionOklahoma    -3000       3000  -1.000    0.500
regionTexas          NA         NA      NA       NA

Residual standard error: 2121 on 1 degrees of freedom
Multiple R-squared:   0.55, Adjusted R-squared:   -0.8
F-statistic: 0.4074 on 3 and 1 DF,  p-value: 0.7848


Now I understand there's a linear dependence problem in the model matrix, which is why regionTexas didn't get coefficients. Supposedly, the Intercept is supposed to be the value for the missing value of the categorical variable; problem is, I'm finding it difficult to interpret that reference value since it's having to do "double duty" as the reference value for more than one categorical variable. Ultimately, I have 4 categorical variables in my real-world application.

How can I fix that problem? I'm willing to pad the data with extra rows, including extra values for the categorical variables, if I can somehow control which level in the factor gets to be the reference level. Ultimately, I have 4 factors for which I need to extend this treatment so as to get coefficients for all the factors.

This post seemed relevant, but didn't quite answer my question. There were several other suggestions on multiple regression with categorical variables, but they don't answer my question.

• I posted an answer to your implicit question of "what is happening here?" As to your explicit question "how can I fix that problem?", it would be good if you could explain what specific problem you are seeing here. Given that everything works as designed, there is no problem that I see, so it would be helpful to know what your expectations were and how they were disappointed. What question are you trying to address in your analysis? Do you want to predict, or test something for significance (what?)? Commented Feb 15 at 14:30
• Looking through the categorical-encoding tag might be helpful. This answer goes into details. Commented Feb 15 at 14:34
• The "problem" with "missing" coefficients can be fixed by using post-modeling tools like those provided by the emmeans package instead of just depending on coefficients reported in the model summary. Set up the model with any proper coding of the categorical predictors, as suggested in the answers, and then the tools in that package can provide predictions based on any combinations of predictor values that you like. You don't need to worry about choice of reference level; the software takes that into account.
– EdM
Commented Feb 16 at 18:01

You need to define a reference level for each separate categorical variable, which will be absorbed into the intercept. (Specifically, R does this automatically, by using the alphabetically first factor level as the reference.)

The intercept is the fitted value for an instance of Local and Louisiana.

An instance of Local and Oklahoma would add (i.e., subtract) the regionOklahoma estimate of -3000.

An instance of National and Louisiana would add the geoNational estimate of 500.

And finally, an instance of National and Oklahoma would add both the regionOklahoma estimate of -3000 and the geoNational estimate of 500.

Everything is as it should be.

In your example, region==Texas exactly where geo==Local, so these two factor levels carry the exact same information, and we can't disentangle the two. R addresses this issue by assigning the estimate to the column that comes first in the data frame, or in the formula. Compare the result of summary(lm(charges~region+geo,rev)), which just switches the order of the two variables - you get an NA elsewhere, because your data is still redundant.

If you get NA estimates in your actual data, then you almost certainly have exactly this issue: collinearity, AKA this factor level can be exactly determined based on the other variables. Hard to be more specific without digging into your actual data. Especially if you have more predictors, chasing the specific constellation that is redundant can be hard, and usually depends on the exact order of terms in the formula or in the data frame, as above.

In a later comment, you explain that you are looking to get "the effect of region==Louisiana". Now, "the effect of" something implies a comparison. If you want to compare the effect of being in Louisiana vs. the effect of being in Texas, then predict with your models, once for Louisiana and once for Texas, and take the difference. If you want the effect of Louisiana vs. a grand average, then fit a model for Louisiana and another one for all possible regions, taking the grand average, and compare the two fits.

In simple cases, yes, this is just the parameter estimate. In more complex cases, you can derive this difference by using contrasts and the parameter estimates. Contrasts will also help you assess the statistical significance of any such difference, but if all you are interested in is the estimated value of the difference, I would say that simply fitting models and subtracting predicted values is easier to understand.

In any case, you will need to decide on what to do with the other predictors. The "overall typical setting" may be very atypical in Louisiana, or in some other state. (I usually recommend Miller & Chapman, 2001.) So it very much depends on getting the question right that you are looking to get an answer for.

• If you get NA estimates in your actual data, then you almost certainly have exactly this issue: collinearity, AKA this factor level can be exactly determined based on the other variables. Hard to be more specific without digging into your actual data... Commented Feb 15 at 15:06
• It might be useful to start with your actual data and cut it down successively by removing entire columns or rows until you find a smallest specimen that still exhibits the issue. Commented Feb 15 at 15:33
• OK, we are back to my original question: why do you need these rows? You can predict just fine with the models as they are. What is the purpose of your analysis? There is no usual way to get these rows, because these rows do not exist in standard contrast coding. Commented Feb 15 at 15:36
• "The effect of" something implies a comparison. If you want to compare the effect of being in Louisiana vs. the effect of being in Texas, then predict with your models and take the difference. If you want the effect of Louisiana vs. a grand average, then fit and compare the two fits. In simple cases, yes, this is just the parameter estimate. In any case, you will need to decide on what to do with the other predictors. The "overall typical setting" may be very atypical in Louisiana. So it very much depends on getting the question right that you are looking to get an answer for. Commented Feb 15 at 16:28
• Ok, this is working. I used the tidyverse::crossing function to create a results dataframe containing all possible combinations of the factors (960 of them), and then simply ran the predict function on all those rows, and now I can compare any row to any other row, to say what the effect of switching from this region to that is, etc. Many thanks! Commented Feb 15 at 19:26

The intercept is the predicted level of the dependent variable when all the independent variables are 0 (however that is coded in your data).

There are a number of ways to parametrize categorical variables. Most will have the same "problem" as you cite (it's not really a problem). In R you can use lm(y ~ x, contrasts=list(x,"contr.treatment") to list them. This site from UCLA has a list of a lot of possible paremeterizations. You might want deviation coding (each level compared to grand mean).

• Well, I beg to differ: it is a problem, because I need coefficients on those rows in the summary output. In the example in my post, I need information about Texas, which the model isn't giving me (unless I'm missing something). I have my doubts about using contrasts; that's not going to fix the linear dependency problem underlying the NAs, is it? Commented Feb 15 at 14:56
• You can change the reference category (R uses alphabetical list by default) or use deviation coding. Helmert coding will also give you a result for Texas, but probably not the one you want. Commented Feb 15 at 15:06
• @AdrianKeister saying "I need information about Texas" in this context doesn't make sense. You could say "I need information about Texas compared to..." which is what the different contrasts would give you. If you just want estimated means for each state from the model you can get them with (eg) emmeans. Commented Feb 15 at 15:41