Using the General Social Survey data on those who report having no religious affiliation, we can fit a logit model in which \(\texttt{male}\) and \(\texttt{cohort}\) are our explanatory variables.
# dependencies -- these packages need to be already installed
library(tidyverse)
library(haven)
library(tulaverse)
options(scipen = 9) # reduce amount of scientific notation in results
# read dataset
gss <- read_dta("../dta/gss_norelig_thru2018.dta") %>%
mutate(relig_none = norelig,
raised_as_none = norelig16)model <- glm(relig_none ~ male + cohort, family="binomial", data=gss)
tula(model)Family: binomial / Link: logit
AIC = 43919.582 Number of obs = 64319
BIC = 43946.797 McFadden R-sq = 0.07391
Log likelihood = -21956.791 Nagelkerke R-sq = 0.1017
─────────────────────────────────────────────────────────────────────────────
relig_none │ Coef Std. Err. z P>|z| [95% Conf Interval]
─────────────────────────────────────────────────────────────────────────────
male │ .5712 .02503 22.82 <.0001 .5221 .6203
cohort │ .0341 .0006714 50.8 <.0001 .03279 .03542
(Intercept) │ -68.92 1.316 -52.38 <.0001 -71.49 -66.34
─────────────────────────────────────────────────────────────────────────────. logit norelig male cohort, nolog
Logistic regression Number of obs = 64,319
LR chi2(2) = 3504.53
Prob > chi2 = 0.0000
Log likelihood = -21956.791 Pseudo R2 = 0.0739
------------------------------------------------------------------------------
norelig | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
male | .5711961 .0250283 22.82 0.000 .5221416 .6202507
cohort | .034104 .0006714 50.80 0.000 .0327882 .0354199
_cons | -68.91602 1.3157 -52.38 0.000 -71.49474 -66.33729
------------------------------------------------------------------------------
Examples of how these logit coefficients can be interpreted are:
The main problem with these interpretations is that a change in “log odds” is not a very meaningful metric to anyone not used to thinking in terms of logit coefficients.
One strategy for a more accessible interpretation of logit model results is to present results using odds ratios.
The 2x2 cross-tab below will help us understand what an odds ratio is. Our binary outcome provides the rows of the table and our binary explanatory variable provides the columns.
tulatab(relig_none, male, data=gss) │ Male? │
Religious none? │ female male │ Total
────────────────┼───────────────────────┼──────────
Some religion │ 32,670 24,058 │ 56,728
│ 90.62 84.49 │ 87.92
│ │
No religion │ 3,380 4,417 │ 7,797
│ 9.38 15.51 │ 12.08
────────────────┼───────────────────────┼──────────
Total │ 36,050 28,475 │ 64,525
│ 100.00 100.00 │ 100.00. tabulate norelig male, col nokey
Religious | Male?
none? | female male | Total
--------------+----------------------+----------
Some religion | 32,670 24,058 | 56,728
| 90.62 84.49 | 87.92
--------------+----------------------+----------
No religion | 3,380 4,417 | 7,797
| 9.38 15.51 | 12.08
--------------+----------------------+----------
Total | 36,050 28,475 | 64,525
| 100.00 100.00 | 100.00
Stata tip: By default Stata produces this “key” with the output to \(\texttt{tabulate}\) indicating what the cell values represent. You can turn this off with the option \(\texttt{nokey}\).
From the above, we can compute the odds of being a religious none separately for men and women:
Recall that odds are \(\Pr(y=1)/\Pr(y=0)\).
When we compare odds between groups, we do not take the difference. Instead, we use the ratio of the two odds. This is the odds ratio.
The odds ratio for men compared to women in this example is the odds for men divided by the odds for women. That is: .1835/.1035 = 1.773
If \(x\) is a binary variable, then the odds ratio is: \[\frac{\textrm{Odds}(x=1)}{\textrm{Odds}(x=0)}\]
In our example: \[\frac{\textrm{Odds}(\textrm{men})}{\textrm{Odds}(\textrm{women})}\]
Note that:
With linear model coefficients, we are used to thinking about a positively signed coefficient as meaning an increase and a negatively signed coefficient as meaning a decrease. Instead, with odds ratios, an odds ratio above 1 means increasing odds and an odds ratio below 1 means decreasing odds.
Important: Odds and odds ratios are not the same thing. This is particularly confusing because the odds is defined as a fraction, \(\frac{\Pr(y=1)}{\Pr(y=0)}\), and so the odds is a ratio. But, the odds ratio is a ratio between two odds.
(If you think this is confusing, in the abstruse world of contingency table analysis there are also odds ratio ratios, which is the ratio between two odds ratios, so be thankful that at least you don’t have to keep that straight.)
Let’s fit a logit model in which the outcome is \(\texttt{relig\_none}\) and the only explanatory variable is \(\texttt{male}\).
model1 <- glm(relig_none ~ male, data=gss, family="binomial")
tula(model1)Family: binomial / Link: logit
AIC = 47010.951 Number of obs = 64525
BIC = 47029.101 McFadden R-sq = 0.01176
Log likelihood = -23503.476 Nagelkerke R-sq = 0.01655
─────────────────────────────────────────────────────────────────────────────
relig_none │ Coef Std. Err. z P>|z| [95% Conf Interval]
─────────────────────────────────────────────────────────────────────────────
male │ .5736 .02438 23.53 <.0001 .5258 .6214
(Intercept) │ -2.269 .01807 -125.6 <.0001 -2.304 -2.233
─────────────────────────────────────────────────────────────────────────────. logit norelig male, nolog
Logistic regression Number of obs = 64,525
LR chi2(1) = 559.47
Prob > chi2 = 0.0000
Log likelihood = -23503.476 Pseudo R2 = 0.0118
------------------------------------------------------------------------------
norelig | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
male | .5735747 .024381 23.53 0.000 .5257889 .6213605
_cons | -2.268581 .0180684 -125.56 0.000 -2.303995 -2.233168
------------------------------------------------------------------------------
The coefficient .574 is the difference in log odds between men and women.
\[ \beta_\mathtt{male} = \ln\mathrm{odds}\left(y|\mathtt{male}=1\right) - \ln\mathrm{odds}\left(y|\mathtt{male}=0\right) = .574 \]
Above, we said the odds ratio for men vs. women was:
\[ \frac{\mathrm{odds}\left(y|\mathtt{male}=1\right)}{\mathrm{odds}\left(y|\mathtt{male}=0\right)} = 1.773 \\ \]
A property of logs is that \(\ln(a) - \ln(b) = \ln(a/b)\). Applying this, the log of the odds ratio is:
\[ \ln\left[\frac{\mathrm{odds}\left(y|\mathtt{male}=1\right)}{\mathrm{odds}\left(y|\mathtt{male}=0\right)}\right] = \ln\mathrm{odds}\left(y|\mathtt{male}=1\right) - \ln\mathrm{odds}\left(y|\mathtt{male}=0\right) \]
that is, the difference in the log odds for men vs. women. If we check, we can see that \(\ln(1.773) \approx .574\) (small difference due to rounding).
log(1.773)[1] 0.572673To get the odds ratio for men vs. women, we just take \(\exp(\beta_{male})\). More generally, the coefficients we estimate in our logit model can be interpreted as logged odds ratios. To get the actual odds ratio, we exponentiate the coefficient: \(\exp(\beta)\).
Now let’s look at the example of cohort as an explanatory variable using the same data.
model2 <- glm(relig_none ~ cohort, data=gss, family="binomial")
tula(model2)Family: binomial / Link: logit
AIC = 44444.177 Number of obs = 64319
BIC = 44462.320 McFadden R-sq = 0.0628
Log likelihood = -22220.089 Nagelkerke R-sq = 0.08675
─────────────────────────────────────────────────────────────────────────────
relig_none │ Coef Std. Err. z P>|z| [95% Conf Interval]
─────────────────────────────────────────────────────────────────────────────
cohort │ .03401 .0006659 51.07 <.0001 .0327 .03531
(Intercept) │ -68.44 1.305 -52.46 <.0001 -71 -65.89
─────────────────────────────────────────────────────────────────────────────. logit norelig cohort, nolog
Logistic regression Number of obs = 64,319
LR chi2(1) = 2977.93
Prob > chi2 = 0.0000
Log likelihood = -22220.089 Pseudo R2 = 0.0628
------------------------------------------------------------------------------
norelig | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cohort | .0340069 .0006659 51.07 0.000 .0327018 .035312
_cons | -68.44223 1.304575 -52.46 0.000 -70.99915 -65.88531
------------------------------------------------------------------------------
The logit model is linear in the log odds. For continuous explanatory variable \(x\), the coefficient \(\beta_x\) in the logit model indicates the expected change in the log odds as \(x\) increases by one unit.
If we use \(c\) to indicate some particular value of \(x\), then \(\beta_x\) is the difference in the log odds when \(x=c+1\) vs. when \(x=c\).
Similar to before, we can write \(\beta_x\) as:
\[ \beta_x = \ln\mathrm{odds}\left(y|x = c+1\right) - \ln\mathrm{odds}\left(y|x = c\right) \]
The fact that \(\ln(a) - \ln(b) = \ln(a/b)\) means that:
\[ \begin{align} \beta_x & = \ln\mathrm{odds}\left(y|x = c+1\right) - \ln\mathrm{odds}\left(y|x = c\right) \\ \\ & = \ln \left[ \frac{\mathrm{odds}\left(y|x = c+1\right)}{\mathrm{odds}\left(y|x = c\right)} \right] \end{align} \]
\(\beta_x\) is the log of the odds ratio for a one-unit increase in \(x\).
Odds ratios are multiplicative changes: to get the odds after a multiplicative change, you multiply the odds before the change by the odds ratio. This is like how, if somebody gets a 10% raise, their new salary is 1.10 times their old salary.
\[\textrm{New odds} = \textrm{Old odds} \times \textrm{Odds ratio}\]
that is:
\[\textrm{Odds}(y|x=c+1) = \textrm{Odds}(y|x=c)\exp({\beta_x})\]
Some general features of the odds ratio and the coefficients of the logit model:
The closer a logit coefficient is to zero, the closer the odds ratio will be to \((1+\beta)\). For instance, in our example above, \(\beta\) for cohort was .034, and the odds ratio was 1.035.
The graph below shows the relationship between the odds ratio (\(\exp(\beta)\)) and \(\beta\). A dashed line is shown for reference at (\(\beta\), \(1+\beta\)). That the two lines are effectively on top of one another when \(\beta\) is between -0.1 and 0.1 shows that adding 1 to the logit coefficient is a very close approximation of the odds ratios when the coefficient is small. The approximation works less and less well as the magnitude of the coefficient increases.
plotdata <- data.frame(x = seq(-1, 1, by = 0.01)) %>%
mutate(y = exp(x))
plot <- ggplot(plotdata, aes(x = x, y = y)) +
geom_line(color = "red3", linewidth = 1) +
geom_abline(intercept = 1, slope = 1,
color = "blue",
linetype = "dashed") +
labs(x = "Logit coefficient",
y = "Odds ratio") +
scale_x_continuous(breaks = seq(-1, 1, by = 0.2)) +
theme_tula() +
theme(
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
plot.title = element_text(size = 16, hjust = 0.5)
) +
coord_cartesian(xlim = c(-1, 1), ylim = c(0, 3))
We will fit a model that includes region as an explanatory variable along with sex and birth year.
After fitting a logit model, the tula() function from the tulaverse package will exponentiate the coefficients for us if we use the exp=TRUE option.
gss <- gss %>%
mutate(reg = as_factor(region4cat))
model <- glm(relig_none ~ male + cohort + reg, data=gss, family="binomial")
tula(model, exp=TRUE)Family: binomial / Link: logit
AIC = 43204.372 Number of obs = 64319
BIC = 43258.802 McFadden R-sq = 0.08912
Log likelihood = -21596.186 Nagelkerke R-sq = 0.1219
─────────────────────────────────────────────────────────────────────────────
relig_none │ exp(b) DMSE z P>|z| [95% Conf Interval]
─────────────────────────────────────────────────────────────────────────────
male │ 1.764 .04448 22.5 <.0001 1.679 1.853
cohort │ 1.035 .0007016 50.72 <.0001 1.034 1.036
reg │
Midwest │ .9984 .03588 -.04508 .9640 .9305 1.071
South │ .6265 .02349 -12.47 <.0001 .5821 .6742
Pacific │ 1.688 .06756 13.09 <.0001 1.561 1.826
(Intercept) │ 7.268e-31 9.647e-31 -52.28 <.0001 5.389e-32 9.801e-30
─────────────────────────────────────────────────────────────────────────────Stata will exponentiate the coefficients for us if we include the \(\texttt{or}\) option for the \(\texttt{logit}\) command. Notice in the output below that the results are labeled “Odds Ratio.”
. logit norelig i.male cohort i.region, or nolog
Logistic regression Number of obs = 64,319
LR chi2(5) = 4225.74
Prob > chi2 = 0.0000
Log likelihood = -21596.186 Pseudo R2 = 0.0891
------------------------------------------------------------------------------
norelig | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
male |
male | 1.763798 .0444756 22.50 0.000 1.678746 1.853158
cohort | 1.034978 .0007016 50.72 0.000 1.033604 1.036354
|
region4cat |
Midwest | .9983811 .0358794 -0.05 0.964 .9304783 1.071239
South | .6264698 .0234915 -12.47 0.000 .5820787 .6742464
Pacific | 1.688234 .0675599 13.09 0.000 1.560879 1.825981
|
_cons | 7.27e-31 9.65e-31 -52.28 0.000 5.39e-32 9.80e-30
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.
Stata: \(\texttt{logit}\) vs. \(\texttt{logistic}\). Stata has two commands that both fit the \(\texttt{logit}\) model: \(\texttt{logit}\) and \(\texttt{logistic}\). The difference is just in the output: by default, \(\texttt{logistic}\) provides odds ratios and \(\texttt{logit}\) provides coefficients. But, you can obtain odds ratios using the \(\texttt{logit}\) command with the \(\texttt{or}\) option. There is nothing wrong with using the \(\texttt{logistic}\) command, but personally for my own work I never do.
The odds ratio for male is 1.76. We can interpret this as:
For increases in the odds, we can talk about odds increasing by a factor of \(F\), or we can talk about the odds increasing by \((F-1)\times100\) percent.
We can interpret the odds ratio of 1.035 for cohort as:
In the case of region, the reference category is living in the Northeast. The odds ratio of .626 for South can be interpreted as:
For decreases in the odds, we can describe odds being lower by a factor of \(F\). In percentage terms, we talk about this change as being \((1-F)\times100\) percent lower.
Elaboration: Decreases in odds: Converting an odds ratio below 1 to a percentage change sometimes confuses people. An analogy: think about the price of goods that go on sale.
Say a store is having a sale on widgets and they are 20% off. A widget is normally $1.50. To figure out the sale price of the widget, you would multiply $1.50 by .8 ($1.20). In this example, the .8 is the factor change, and that it is less than one means the new price is lower than the original price. To get the percentage difference, we would subtract .8 from 1, getting .2, meaning 20% off.
Logit model coefficients may be interpreted directly as changes in log odds, but this is not a very clear interpretation to many audiences.
If we exponentiate logit model coefficients, we get odds ratios, which can be interpreted as a factor change in the odds that y=1.
Factor changes may be easily converted to percentage changes. A factor change of 1.30 implies a 30% increase in odds, whereas a factor change of .70 implies a 30% decrease.