How can I estimate variable influence of the factors in my model rather than just factor levels?

Question

I created a regression model And i want to estimate an influence analysis for each factor.

Meaning take the square Wald-estimation (z-value) for a specific factor and divide it by sum of squares of their Wald-estimation. and that how I estimate the influence of specific factor.

My problem is that the factors are divided by their levels.

I will give an example:

model<-glm(formula = form,
            family = binomial("logit"),
            data   = Train)
View(summary(model)$coefficients)

In the table we can see that the factor dom_time_Colnames is divided into 4 levels. Same thing happened with first_byte_downdload_Colnames.

I want to take the factors z-values and not their levels z-values.

How I do it? anova() is a good idea but it doesn't stop running for me. I search for creative solution that give me output like the the z-value in glm summary or deviance in anova for the all factor and not for theirs levels.

Here is a reproducible example:

Data<-data.frame(Species=iris$Species)
for(i in 1:ncol(iris)){
  if(is.numeric(iris[,i])){
    result=quantile(x = iris[,i],probs = seq(0,1,0.1)) 
    out<-cut(iris[,i], breaks = unique(result),include.lowest = TRUE)
    Data<-data.frame(Data,out)
    colnames(Data)[length(Data)]<-colnames(iris)[i]
  } else {
    next()
  }
}
Data$y<-rbinom(n = nrow(Data),size = 1,prob = 0.1)
form<-formula(y~.)
model<-glm(formula = form,
            family = binomial("logit"),
            data = Data)
View(summary(model)$coefficients)

we can see that the factor sepal or Petal is divided to its levels.

Answer 1

As I already mentioned is R converting factor variables into dummy variables when using glm. This makes I think intuitive sense because otherwise it would use arbitrary numeric values which are used for the coding of the factor levels.

You can use an F-test to look if multiple variables together are significant. Here we want to test whether all dummies together are significant (in other words the whole factor variable). This can be done by first fitting one model with the factor included and one where you skipped the factor variable. In a second step you can pass both models to the test function (here we use lmtest::waldtest).

Here is my suggestion:

One remark: I changed the probability of the y variable from 0.1 to 0.3. Otherwise the P-value almost always becomes 1.

library(lmtest)

## The full model containing all variables
full_formula      <- formula(y~ Species + Sepal.Length + Sepal.Width + Petal.Length + Petal.Width)
full_model        <- glm(formula = full_formula,
                         family = binomial("logit"),
                         data = Data)

## Reduced models (leaving out each variable once)
reduced_models <- list()
for(i in 1:(ncol(Data)-1)) {
  reduced_formulas  <- as.formula(paste0("y ~ ", paste(names(Data)[c(-i ,-ncol(Data))], collapse= " + " )))
  reduced_models[[i]] <- glm(formula = reduced_formulas,
                         family = binomial("logit"),
                         data = Data)
}
## Test the full model against all reduced models
result        <- lapply(reduced_models, waldtest, full_model)

##Add the name of the tested variable
names(result) <- names(Data[ ,-ncol(Data)])

The first three results:

> result[1:3]
$Species
Wald test

Model 1: y ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
Model 2: y ~ Species + Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
  Res.Df Df      F Pr(>F)
1    115                 
2    113  2 0.1116 0.8945

$Sepal.Length
Wald test

Model 1: y ~ Species + Sepal.Width + Petal.Length + Petal.Width
Model 2: y ~ Species + Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
  Res.Df Df      F Pr(>F)
1    122                 
2    113  9 0.7051 0.7031

$Sepal.Width
Wald test

Model 1: y ~ Species + Sepal.Length + Petal.Length + Petal.Width
Model 2: y ~ Species + Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
  Res.Df Df      F Pr(>F)
1    121                 
2    113  8 0.1743 0.9939

The first list entry tested whether the factor variable species is significant, the second whether Sepal.Length is significant and so on.

It is no surprise that all tests have high P-values since the y variable was completely random.

*edit

You can also use anova():

result_anova        <- lapply(reduced_models, anova, full_model)
names(result_anova) <- names(Data[ ,-ncol(Data)])

Only the result for the first variable:

> result_anova[[1]]
Analysis of Deviance Table

Model 1: y ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
Model 2: y ~ Species + Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
  Resid. Df Resid. Dev Df Deviance
1       115     149.16            
2       113     146.87  2   2.2914

Answer 2

Given your clarification in the comments that what you mean by "variable importance" is actually statistical significance, I suggest you use a likelihood ratio test. For simple-to-simple hypotheses, these are uniformly most powerful (UMP) tests and for composite hypotheses are still ideal.

The Wald test option, which other @Alex suggests above, is another option. The LRT and Wald test are asymptotically equivalent. Although, for simple v simple hypotheses, the LRT is the UMP and thus preferred.

# using your Data argument, here are nested models:
model1<-glm(y ~ ., family = binomial("logit"), data = Data)
model2 <- glm(y ~ Sepal.Length + Sepal.Width + Petal.Length, 
              family= binomial("logit"), data= Data)

# LRT model1 vs the null hypothesis
pchisq(model1$null.deviance - model1$deviance, 
       model1$df.null - model1$df.residual, lower= FALSE)
[1] 0.1721306

# nested model LRT
pchisq(model2$deviance - model1$deviance, 
       model2$df.residual - model1$df.residual, lower= FALSE)
[1] 0.05313724

In both cases, since pchisq > 0.05 (a standard default significance level) we do not reject the null hypothesis that the additional covariates are needed. For this dataset, we would choose the null model. More relevant to your question: for the nested model LRT, the additional terms (eg. model1) are not needed.

Again, this is not a metric of "variable importance" (eg. Breiman 2001). It's a metric assessing goodness of fit.