Comparing weights between significant predictors in GLM with binomial outcome

Question

I need to compare the explanatory power of several predictor variables $\{x_1,x_2,x_3\} \in \mathbb{R}$ for a response variable $y$. For this problem, we can use logistic regression as a generalized linear model:

$$logit(p_i) = \beta_0 + \beta_1 x_1,i + \beta_2 x_2,i + \beta_3 x_3,i + \beta_4 x_4,i$$

$\forall i \in [1,N]$, where $N$ is the number of observations.

Any statistical software can fit the model and provide weights and p-values for each predictor in the model. What is the simplest method to compare the resulting weights that are significant? From my reading, it appears one method of making the weights comparable in standard GLM approaches is to scale $\beta_j$ by the standard deviation of the predictor variable:

$$\beta_j' = \beta_j \sigma_j$$

However, is this approach still valid for a GLM with binomial response variable? If not, is there another straightforward way to do this?

Answer 1

In my opinion, the $\beta'_j$ measure you propose doesn't have a useful role in understanding logistic regression and you would be better of using traditional generalized linear model summaries.

Here is a simple logistic regression example to show that the $\beta_j'$ measures are misleading:

> x1 <- seq(from=-1, to=1, length=1000)
> beta <- 0.5
> p <- exp(x1*beta) / (1+exp(x1*beta))
> y <- rbinom(1000, prob=p, size=1)
> x2 <- rep(0,1000)
> x2[500] <- 1
> y[500] <- 1
> fit <- glm(y~x1+x2, family=binomial, maxit=100, epsilon=1e-18)
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred

In this example x1 is highly significant whereas x2 is not significant at all. x2 has little or no predictive power: We can see this either from the t-statistics and Wald p-values:

> summary(fit)

Call:
glm(formula = y ~ x1 + x2, family = binomial, maxit = 100, epsilon = 1e-18)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4218  -1.1601   0.9551   1.1401   1.3907  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 3.481e-02  6.400e-02   0.544    0.587    
x1          5.234e-01  1.117e-01   4.686 2.78e-06 ***
x2          3.153e+01  6.711e+07   0.000    1.000    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1386.0  on 999  degrees of freedom
Residual deviance: 1362.2  on 997  degrees of freedom
AIC: 1368.2

Number of Fisher Scoring iterations: 30

or from analysis of deviance and likelihood ratio tests:

> anova(fit, test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: y

Terms added sequentially (first to last)


     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                   999     1386.0              
x1    1  22.4111       998     1363.6 2.201e-06 ***
x2    1   1.3513       997     1362.2     0.245    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Nevertheless x2 has a much higher $\beta_j'$ measure than x1:

> beta1 <- coef(fit)[2]
> beta2 <- coef(fit)[3]
> beta1*sd(x1)
[1] 0.3026638
> beta2*sd(x2)
[1] 0.9971314

Here the value of $\beta_2'$ would actually be infinite if we were able to fit the glm in exact arithmetic rather than on a floating point computer.

Measuring predictive power in a meaningful way requires comparing the predicted probabilities to the actual outcome values. One useful summary is the area under the receiver operating curve (auROC):

The fitted model with x1 has auROC of 59%, regardless of whether x2 is in the model as well:

> library(limma)
> auROC(y, fitted(fit))
[1] 0.5870262
> fit1 <- glm(y~x1, family=binomial)
> auROC(y,fitted(fit1))
[1] 0.5861539

x2 alone is hardly better than random:

> fit2 <- glm(y~x2, family=binomial)
> auROC(y, fitted(fit2))
[1] 0.5009823