3. Logit Models

A 2 by 2 Table

Consider the data on contraceptive use by desire for more children on Table 3.2 (page 14 of the notes). We can read these data into our software as 2 binomial observations. To make life easier, I will enter desire for more children as a dummy variable that takes the value 1 for women who want no more children and 0 otherwise.

> library(dplyr)
> cuse <- data.frame(matrix(c(
+ 0, 219, 753,
+ 1, 288, 347), 2, 3, byrow=TRUE))
> names(cuse) <- c("nomore", "using","notUsing")
> cuse <- mutate(cuse, n = using + notUsing)
> cuse

  nomore using notUsing   n
1      0   219      753 972
2      1   288      347 635

Testing Homogeneity

Let us start by fitting the null model. With grouped binomial data the outcome is a matrix with two columns, representing the number of “successes” and “failures”, in this case users and non-users of contraception. The binomial denominator is calculated internally as the sum of these two columns.

> cuse <- mutate(cuse, Y = cbind(using, notUsing))
> m0 <- glm(Y ~ 1, family=binomial, data = cuse)
> summary(m0)

Call:
glm(formula = Y ~ 1, family = binomial, data = cuse)

Deviance Residuals: 
     1       2  
-6.240   7.262  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.77455    0.05368  -14.43   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 91.674  on 1  degrees of freedom
Residual deviance: 91.674  on 1  degrees of freedom
AIC: 107.54

Number of Fisher Scoring iterations: 4

> sum(residuals(m0, type="pearson")^2) # Pearson's chi-squared

[1] 92.64424

The estimate of the constant is simply the logit of the overall proportion using contraception, say p=y/n, and the standard error is the square root of 1/y + 1/(n-y). You may want to check these results by hand. To compute logits you may use qlogis(), with inverse plogis(), for quantiles and probabilities of the standard logistic distribution. These specialized functions are more accurate than just computing log(p/(1-p)).

The deviance is 91.67 on one d.f., providing ample evidence that the null model does not fit the data. Thus, we reject the hypothesis that the probability of using contraception is the same in the two groups.

An alternative test is Pearson’s chi-squared, which is 92.64 on one d.f., and leads to the same conclusions. These two tests are asymptotically equivalent, so they tend to give similar results in large samples.

The Odds Ratio

Let us now fit the model with “want no more” children as the predictor. This model is saturated for this dataset, using two parameters to model two probabilities:

> m1 <- glm(Y ~ nomore, family=binomial, data=cuse)
> summary(m1)

Call:
glm(formula = Y ~ nomore, family = binomial, data = cuse)

Deviance Residuals: 
[1]  0  0

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.23499    0.07677 -16.086   <2e-16 ***
nomore       1.04863    0.11067   9.475   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9.1674e+01  on 1  degrees of freedom
Residual deviance: 1.7986e-14  on 0  degrees of freedom
AIC: 17.87

Number of Fisher Scoring iterations: 2

> exp(coef(m1)["nomore"])

  nomore 
2.853737 

The constant corresponds to the log-odds of using contraception among whomen who do want more children, and the coefficient of nomore is the difference in log-odds between the two groups.

Exponentiating this coefficient we get an odds ratio of about three. Contrary to popular belief, this does not mean that “women who want no more children are three times more likely to use contraception”. There are two errors in this interpretation.

First, and more importantly, it is the odds of using contraception among women who want no more children that are three times those of women who want more, not the probability, which is what’s usually understood by “likelihood”. The interpretation would be approximately correct if the event under study was rare, because if p is small then 1-p is close to one and the odds ratio is approximately the same as the relative risk. Here the observed proportions are 0.454 and 0.225, and the ratio is 2.01, so women who want no more children are twice as likely to use contraception as those who want more.

Second, even if the probability was tripled, that would make the women three times as likely, or two times more likely, to use contraception, not three times more likely. In this case the probability is doubled, and that makes women twice as likely, not two times more likely.

Testing Significance

The z-statistic is as reported on page 16 of the notes. Let us square it:

> b <- coef(m1)
> se <- sqrt(diag(vcov(m1)))
> (b[2]/se[2])^2

  nomore 
89.77765 

This is Wald’s chi-squared statistic for the hypothesis that the coefficient of nomore is zero, or equivalently that the odds-ratio is one.

We can also compare these models using a likelihood ratio test.

> -2*(logLik(m0) - logLik(m1))

'log Lik.' 91.6744 (df=1)

Can you explain why we get 91.67, which is the deviance of the null model? Hint: What’s the deviance of the current model?

A third test of the effect of want no more is given by Pearson’s chi-squared statistic, which we calculated earlier as 92.64. This is equivalent to the standard z-test for comparing two proportions if you use the pooled proportion to estimate the standard error.

All three statistics are different, but they are asymptotically equivalent. In our example they are also close in value and lead to the same overwhelming rejection of the hypothesis that the probability of using contraception is the same in the two groups.

Confidence Intervals

You can obtain a confidence interval in R by calling the confint() function, which uses a profile log-likelihood. You can obtain the more conventional confidence intervals by calling confint.default(). Let us obtain a confidence interval for the odds ratio using both methods.

> exp(confint(m1,"nomore"))

Waiting for profiling to be done...

   2.5 %   97.5 % 
2.298942 3.548111 

> exp(confint.default(m1,"nomore"))

          2.5 %   97.5 %
nomore 2.297258 3.545015

With 95% confidence, the odds of using contraception among women who want no more kids are between 2.3 and 3.5 times as high as for women who want more children.

Note that the standard confidence bounds for the odds ratio are not calculated by adding and subtracting twice its standard error. Instead, the calculation is done in the logit scale and the results are then exponentiated. This is done because the normal approximation is more accurate (and makes more sense) in the logit scale, which has no range restrictions.

Exercise. Calculate the conventional z-test for comparing the proportions using contraception in the two groups and verify that the square coincides with Pearson’s chi-squared statistic.

Updated fall 2022