# Ran a logistic regression and calculated the probabilities but they appear too high

I ran a logistic regression in R using driving data from about 10,000 people. The model included age, years of driving experience, as well as 4 driving test results. The dependent variable was whether or not they had been involved in a crash recently (yes or no, a categorical variable).

The coefficients of the model are given below:

                    Estimate     Std.Err       z value     Pr(>|z|)
(Intercept)        -1.450041     0.207144      -7.000      2.56e-12 ***
riding experience  -0.014115     0.003697      -3.818      0.000134 ***
age                -0.034544     0.003608      -9.575       < 2e-16 ***
test 1              0.261485     0.088645       2.950      0.003180 **
test 2              0.090102     0.051328       1.755      0.079184 .
test 3              0.228918     0.073666       3.108      0.001887 **
test 4              0.070106     0.063652       1.101      0.270729


Firstly, with 10,000 people am I right in thinking that p-values aren't going to be that useful?

I calculated the probabilities of being involved in a crash with a 1 unit increase in each variable by doing exp(variable) to get the odds and then, probability = odds/(1+odds). It gave me:

(Intercept)        0.1899952
ridingexp          0.4964712
age                0.4913648
test 1             0.5650012
test 2             0.5225104
test 3             0.5569810
test 4             0.5175193


These seem awfully high! It is like saying that an increase in age of 1 year makes you 49% less likely to be involved in a crash? Surely that can't be right.

I believe the mistake you are making is the connection between the coefficients and the probability.

Defining the linear predictor (i.e. the estimated log odds) as

$\eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ...$

You are right that $\exp(\beta_1)$ represents the fold difference in odds associated with a one unit increase in $x_1$. So for age, a one year increase in age is associated with a exp(-0.0345) = 0.966 fold change in the odds of being in an accident. As you expected, not a huge difference.

It appears that you just took 0.966 as the odds and converted this to a probability. This was the fundamental mistake: 0.996 is not the odds, but the multiplicative difference in odds.

I assume you want to go back to the probability scale because that's the scale that's most natural to you. However, it's not so simple to convert a change in the odds to a change in the probability; how much the change a multiplicative change in the odds changes the probability (on the multiplicative scale) is not constant, but rather depends on the current value of the odds.

For example, suppose we have that the base odds are 1, and multiply this by 2. The corresponding probabilities are 1/(1 + 1) = 1/2 and 2/(2+1) = 2/3. Now suppose our base odds are 2 and we multiply this by 2. The corresponding probabilities are 2/(2+1) = 2/3 and 4/(4+1) = 4/5.

So it's not so simple to summarize a change in the odds on the probability scale.

As @CliffAB mentions, interpreting changes in the odds on the predicted probabilities scale can be tough. But this is where graphing comes in handy. For instance, let's generate some data following a simple version of the model you provided, then estimate the model:

set.seed(2)
n <- 1000
experience <- rnorm(n)
age <- sample(16:99, n, replace=TRUE)
z <- -1.45 + (-0.01 * experience) + (-0.03 * age)
y <- rbinom(n, 1, plogis(z))
mod <- glm(y ~ experience + age, family="binomial")
summary(mod)

Call:
glm(formula = y ~ experience + age, family = "binomial")

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5963  -0.3857  -0.2556  -0.1709   3.0434

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.085033   0.332718  -3.261  0.00111 **
experience  -0.016846   0.140535  -0.120  0.90459
age         -0.037357   0.007272  -5.137 2.79e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 402.9  on 999  degrees of freedom
Residual deviance: 370.5  on 997  degrees of freedom
AIC: 376.5

Number of Fisher Scoring iterations: 6


We can then extract the predicted probabilities for each data point in our data set with the predict function (here just looking at the first 6 observations):

> head(predict(mod, type="response"))
1          2          3          4          5          6
0.02194610 0.02319195 0.01355816 0.04226472 0.09935722 0.06117067


Using type="response" is the same as taking the linear predictor with predict and running it through plogis:

> head(cbind(predict(mod, type="response"),
+            plogis(predict(mod))))
[,1]       [,2]
1 0.02194610 0.02194610
2 0.02319195 0.02319195
3 0.01355816 0.01355816
4 0.04226472 0.04226472
5 0.09935722 0.09935722
6 0.06117067 0.06117067


Now, we can plot the predicted probability of having an accident at different ages, according to our data set:

> plot(age, predict(mod, type="response"), ylab="Predicted Probability")


However, this doesn't keep experience constant, so we can create a new data set of "covariate profiles" we may be interested in, where the only thing that changes is age, holding experience at its mean:

new <- data.frame("(Intercept)"=1,
"experience"=mean(experience),
"age"=16:99)

plot(predict(mod, newdata=new, type="response"), type="l", lwd=2, col="red",
ylab="Predicted Probability")