VIF for Categorical Variable with More Than 2 Categories

Question

I'm trying to detect multicollinearity using VIF in both Python and R. Based on my knowledge, the VIF should be less than 10 if there is no multicollinearity. However, for the categorical variable with more than 2 categories, the VIF of some categories are very high. My data include the variable more than 10 categories. Here is what I did in Python:

y, X = dmatrices('InvoiceUnitPrice~NewWidth+NewLength+NewThickness+InvoiceQuantity+Weight+SUPP_CD', data=ga_for_model, return_type='dataframe')
vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif["features"] = X.columns
vif

Out[198]: 
    VIF Factor            features
0   171.420478           Intercept
1    16.307844         SUPP_CD[W2]
2     7.677684         SUPP_CD[W3]
3     5.200108         SUPP_CD[Y0]
4     1.033676         SUPP_CD[Y4]
5     1.324480         SUPP_CD[Y1]
6     1.030234         SUPP_CD[H0]
7     1.220017         SUPP_CD[L0]
8     1.067945         SUPP_CD[L1]
9     1.163532         SUPP_CD[X1]
...   ...              ...
83    2.692464            NewWidth
84    2.729983           NewLength
85    1.744165        NewThickness
86    1.426814     InvoiceQuantity
87    1.079581              Weight

[88 rows x 2 columns]

The SUPP_CD[W2] has a very high VIF as it showed. Then I use vif() from car package in R to run the result again:

> vif(model)
                             GVIF Df GVIF^(1/(2*Df))
for_R$NewWidth           2.780087  1        1.667359
for_R$NewLength          2.834620  1        1.683633
for_R$SUPP_CD         7419.836402 82        1.055845
for_R$NewThickness       2.367231  1        1.538581
for_R$Type            8406.690333 21        1.240062
for_R$InvoiceQuantity    1.495487  1        1.222901
for_R$Weight             1.142044  1        1.068665

The difference between these two results makes me confused. For the result in R, I've looked up the difference between GVIF Df and GVIF^(1/(2*Df)) from
Which variance inflation factor should I be using: $\text{GVIF}$ or $\text{GVIF}^{1/(2\cdot\text{df})}$?

"Georges Monette and I introduced the GVIF in the paper "Generalized collinearity diagnostics," JASA 87:178-183, 1992 (link). As we explained, the GVIF represents the squared ratio of hypervolumes of the joint-confidence ellipsoid for a subset of coefficients to the "utopian" ellipsoid that would be obtained if the regressors in this subset were uncorrelated with regressors in the complementary subset. In the case of a single coefficient, this specializes to the usual VIF. To make GVIFs comparable across dimensions, we suggested using GVIF^(1/(2*Df)), where Df is the number of coefficients in the subset. In effect, this reduces the GVIF to a linear measure, and for the VIF, where Df = 1, is proportional to the inflation due to collinearity in the confidence interval for the coefficient."

So I think the results from the R points out no multicollinearity by looking at GVIF^(1/(2*Df)) (Please correct me if I'm wrong.)

But for the result in Python, it gives VIF for each category. I don't know how to interpret them and how to deal with them.

Although Paul Allison introduced 3 situations that can ignore high VIF values in When Can You Safely Ignore Multicollinearity?, he mentions dummy variables only. Not suitable for my problem.

3. The variables with high VIFs are indicator (dummy) variables that represent a categorical variable with three or more categories. If the proportion of cases in the reference category is small, the indicator variables will necessarily have high VIFs, even if the categorical variable is not associated with other variables in the regression model.
   Suppose, for example, that a marital status variable has three categories: currently married, never married, and formerly married. You choose formerly married as the reference category, with indicator variables for the other two. What happens is that the correlation between those two indicators gets more negative as the fraction of people in the reference category gets smaller. For example, if 45 percent of people are never married, 45 percent are married, and 10 percent are formerly married, the VIFs for the married and never-married indicators will be at least 3.0.

I know I could convert categorical variables to dummy variables, but the VIF function still works without the conversion. Please help. Thank you!

Answer 1

The "generalized variance inflation factors" (GVIF) implemented in the vif() function of the R car package were designed by Fox and Monette specifically to handle situations like this, where there are groups of predictor variables that should be considered together rather than separately. Such situations include multi-level categorical variables and polynomial terms in a single variable.

The standard VIF calculation described on the Wikipedia page (and evidently as implemented in the Python variance_inflation_factor() function) treats each predictor separately. A $k$-level categorical variable then counts as $k-1$ predictors, and the result of that type of VIF calculation will depend on how that variable is coded, specifically which category is considered the reference level. Allison alluded to that in the post you linked, recommending use of the most frequent category as the reference when performing that type of VIF calculation.

The GVIF approach provides a combined measure of collinearity for each group of predictors that should be considered together, like each of your multi-level categorical variables. It does this in a way that is independent of the details of how those predictors are coded. The GVIF^(1/(2*Df)) calculation then provides comparability among predictor sets having different dimensions.

So in your case the GVIF approach is most applicable, and there seems to be no substantial evidence of multicollinearity.