Generalized Linear Models in R, Part 3: Plotting Predicted Probabilities

by David Lillis, Ph.D.

In our last article, we learned about model fit in Generalized Linear Models on binary data using the glm() command. We continue with the same glm on the mtcars data set (regressing the vs variable on the weight and engine displacement).

Now we want to plot our model, along with the observed data.

Although we ran a model with multiple predictors, it can help interpretation to plot the predicted probability that vs=1 against each predictor separately.  So first we fit a glm for only one of our predictors, wt.

model_weight

summary(model_weight)

Call:
glm(formula = vs ~ wt, family = binomial, data = mtcars)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9003  -0.7641  -0.1559   0.7223   1.5736

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   5.7147     2.3014   2.483  0.01302 * 
wt           -1.9105     0.7279  -2.625  0.00867 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 43.860  on 31  degrees of freedom
Residual deviance: 31.367  on 30  degrees of freedom
AIC: 35.367

Number of Fisher Scoring iterations: 5

To plot our model we need a range of values of weight for which to produce fitted values. This range of values we can establish from the actual range of values of wt.

range(mtcars$wt)
 [1] 1.513 5.424

A range of wt values between 0 and 6 would be ideal. So we create a sequence of values between 0 and 6 in increments of 0.01. Joining such a large number of closely spaced points will give a smooth appearance to our model.

xweight <- seq(0, 6, 0.01)

Now we use the predict() function to create the model for all of the values of xweight.

yweight <- predict(model_weight, list(wt = xweight),type="response")

Now we plot.

plot(mtcars$wt, mtcars$vs, pch = 16, xlab = "WEIGHT (g)", ylab = "VS")

lines(xweight, yweight)

We can do the same for displacement.

model_disp

summary(model_disp)

Call:
glm(formula = vs ~ disp, family = binomial, data = mtcars)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8621  -0.4143  -0.0872   0.5694   1.8157

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept)  4.137827   1.389354   2.978  0.00290 **
disp        -0.021600   0.007131  -3.029  0.00245 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 43.860  on 31  degrees of freedom
Residual deviance: 22.696  on 30  degrees of freedom
AIC: 26.696

Number of Fisher Scoring iterations: 6

range(mtcars$disp)
[1]  71.1 472.0

xdisp

ydisp

plot(mtcars$disp, mtcars$vs, pch = 16, xlab = "DISPLACEMENT (cubic inches)", ylab = "VS")

lines(xdisp, ydisp)

We can see that for both predictors, there is a negative  relationship between the probability that vs=1 and the predictor variable.  As the predictor increases, the probability decreases.

That wasn’t so hard! In our next article, we will look at other applications of the glm() function.

About the Author: David Lillis has taught R to many researchers and statisticians. His company, Sigma Statistics and Research Limited, provides both on-line instruction and face-to-face workshops on R, and coding services in R. David holds a doctorate in applied statistics.

Tagged as: Generalized Linear Model, GLM, Logistic Regression, R, sigmoidal curve