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?

  • $\begingroup$ The simplest method to assess the importance of the different predictors is to look at the p-values. Small means important. Can you explain in your question why or why not the p-values already serve your purpose? What else are you after? $\endgroup$ – Gordon Smyth May 8 at 0:43
  • $\begingroup$ Thanks for the comment - my question was meant to focus on the set of significant predictors and comparing between these. I'll edit accordingly now. $\endgroup$ – user59071 May 8 at 2:38
  • $\begingroup$ The edits don't help. You've missed my point, which was to examine the size of the p-values, not just apply preset significant cutoff. I will write an answer to clarify. $\endgroup$ – Gordon Smyth May 9 at 8:21

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)

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  

             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
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.