# How to check for multicollinearity in Poisson regression model with robust error variance

I have created a Poisson regression model with robust error variance (https://academic.oup.com/aje/article/159/7/702/71883) to calculate relative risks.

This is the Poisson regression model:

glm.poisson <- glm(new_anydiagnosis ~ gender + age + socioeco_status+dsdm_category+family_history+mental_before_hiv+relations+sexual_life+stigma_discrimination,
data=finaldata)


To calculate the robust standard errors, I have used the package "sandwich" and "lmtest":


library("sandwich")
library("lmtest")

glm.robust <- coeftest(glm.poisson, vcov = sandwich)


And I get the following estimates:

z test of coefficients:

Estimate  Std. Error  z value  Pr(>|z|)
(Intercept)             -1.6661719   0.6595846  -2.5261  0.011534 *
gender2                  0.2785553   0.2568445   1.0845  0.278130
age                     -0.0039597   0.0087759  -0.4512  0.651846
socioeco_status2        -0.2993100   0.2941404  -1.0176  0.308880
socioeco_status3        -0.2157970   0.3404303  -0.6339  0.526149
socioeco_status5       -14.2642777   0.5838599 -24.4310 < 2.2e-16 ***
dsdm_category2           0.4010045   0.2631345   1.5240  0.127521
dsdm_category3           0.3511796   0.5260824   0.6675  0.504429
dsdm_category4           0.6580096   0.3078882   2.1372  0.032584 *
family_history2         -0.5735175   0.2567557  -2.2337  0.025502 *
mental_before_hiv2       0.8308926   0.5205853   1.5961  0.110472
mental_before_hiv3       1.4173136   0.5843978   2.4253  0.015298 *
relations2               0.0348879   0.2586927   0.1349  0.892721
relations3              27.5722112   1.6086563  17.1399 < 2.2e-16 ***
sexual_life2             0.5890264   0.2753578   2.1391  0.032425 *
sexual_life3           -13.7536832   0.7330256 -18.7629 < 2.2e-16 ***
stigma_discrimination2   0.7728856   0.2483147   3.1125  0.001855 **
stigma_discrimination3 -13.9276696   1.0337112 -13.4735 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


However, now I want to check if multicollinearity between the variables of the model exists. For this, I was considering using the variance inflation factor ("vif").

But when trying this, I do not get estimates of the "vif" for all of the parameters in the model. I only get one estimate instead of estimates of the vif for all variables of the model:

> vif(glm.robust)
Estimate Std. Error    z value   Pr(>|z|)
11.381583   2.220626  10.048533   1.309115


The problem (I think) is that my glm.robust model is not defined as a model, but regarded as a vector (?), since it is defined based on the "coeftest".

Do you know how I can test for multicollinearity in the Poisson regression model with robust error variance?

I was also considering finding the "vif" for the parameters of the original Poisson model (glm.poisson), but as this model does not use the robust error variance, I am not sure if it would be fine to do.

I am not at all experienced in using Poisson regression, so maybe I am missing something.

The most useful test for multicollinearity in a generalized linear model, particularly with multi-level categorical predictors like yours, is based on the variance-covariance matrix of the coefficient estimates. If the vif() function you invoked is from the R car package, then you are on the right track.
As I understand the coeftest() function, however, it does not return the robust variance-covariance matrix. It generates and uses that matrix for significance testing of the individual coefficients, but then discards it.
With the car vif() function, you need to pass to it an object that has the variance-covariance matrix available via vcov(object). This might already be implemented somewhere, but you could generate the robust matrix directly from the sandwich package, define a new class for that matrix, and write a vcov() function for that class that just returns the matrix. Then submit the newly-classed matrix to vif(). Or you could adapt the code of vif() as shown in the page linked above to work directly on the matrix.