# Exact logistic regression

Exact logistic regression is used to model binary outcome variables in which the log odds of the outcome is modeled as a linear combination of the predictor variables. It is used when the sample size is too small for a regular logistic regression (which uses the standard maximum-likelihood-based estimator) and/or when some of the cells formed by the outcome and categorical predictor variable have no observations. The estimates given by exact logistic regression do not depend on asymptotic results.

```require(elrm)
```
```## Loading required package: elrm
```

## Example of exact logistic regression

Suppose that we are interested in the factors that influence whether or not a high school senior is admitted into a very competitive engineering school. The outcome variable is binary (0/1): admit or not admit. The predictor variables of interest include student gender and whether or not the student took Advanced Placement calculus in high school. Because the response variable is binary, we need to use a model that handles 0/1 outcome variables correctly. Also, because of the number of students involved is small, we will need a procedure that can perform the estimation with a small sample size.

## Admitted students data

The data for this exact logistic data analysis include the number of students admitted, the total number of applicants broken down by gender (the variable `female`), and whether or not they had taken AP calculus (the variable `apcalc`). Since the dataset is so small, we will read it in directly.

```dat <- read.table(text = "
female apcalc admit num
0 0 0 7
0 0 1 1
0 1 0 3
0 1 1 7
1 0 0 5
1 0 1 1
1 1 0 0
1 1 1 6",
```

The `num` variable indicates frequency weight. We use this to expand the dataset and then look at some frequency tables.

```## expand dataset by repeating each row num times and drop the num
## variable
dat <- dat[rep(1:nrow(dat), dat\$num), -4]
```

```## look at various tables
xtabs(~female + apcalc, data = dat)
```
```##       apcalc
## female  0  1
##      0  8 10
##      1  6  6
```

```xtabs(~female + admit, data = dat)
```
```##       admit
## female  0  1
##      0 10  8
##      1  5  7
```

```xtabs(~apcalc + admit, data = dat)
```
```##       admit
## apcalc  0  1
##      0 12  2
##      1  3 13
```

```xtabs(~female + apcalc + admit, data = dat)
```
```## , , admit = 0
##
##       apcalc
## female 0 1
##      0 7 3
##      1 5 0
##
## , , admit = 1
##
##       apcalc
## female 0 1
##      0 1 7
##      1 1 6
```

The tables reveal that 30 students applied for the Engineering program. Of those, 15 were admitted and 15 were denied admission. There were 18 male and 12 female applicants. Sixteen of the applicants had taken AP calculus and 14 had not. Note that all of the females who took AP calculus were admitted, versus only 70% the males.

## Analysis methods you might consider

Below is a list of some analysis methods you may have encountered. Some of the methods listed are quite reasonable, while others have either fallen out of favor or have limitations.

• Exact logistic regression - This technique is appropriate because the outcome variable is binary, the sample size is small, and some cells are empty.
• Regular logistic regression - Due to the small sample size and the presence of cells with no subjects, regular logistic regression is not advisable, and it might not even be estimable.
• Two-way contingency tables - with Fisher's exact test due to small expected values.

## (Approximate) Exact logistic regression

Let's run an (approximate) exact logistic analysis using the `elrm` command in the `elrm` package. This is based on MCMC sampling. It requires a collapsed data set with number of trials and number of successes, so we make that first.

```(x <- xtabs(~admit + interaction(female, apcalc), data = dat))
```
```##      interaction(female, apcalc)
## admit 0.0 1.0 0.1 1.1
##     0   7   5   3   0
##     1   1   1   7   6
```
```# view cross tabs
```

```(cdat <- cdat <- data.frame(female = rep(0:1, 2), apcalc = rep(0:1, each = 2),
admit = x[2, ], ntrials = colSums(x)))
```
```##     female apcalc admit ntrials
## 0.0      0      0     1       8
## 1.0      1      0     1       6
## 0.1      0      1     7      10
## 1.1      1      1     6       6
```
```# view collapsed data set
```

Now we can estimate the approximate logistic regression using `elrm` and MCMC sampling. We will do 22,000 iterations with a 2,000 burnin for a final chain of 20,000. Note that for the combined model of female and apcalc, we use a chain of 5 million. This is because for inference, each effect needs at least 1,000, but because the conditional joint distribution is degenerate, for the female effect the ratio of useable trials is low, meaning that to achieve over 1,000, the total iterations must be extremely high.

```## model with female predictor only
m.female <- elrm(formula = admit/ntrials ~ female, interest = ~female, iter = 22000,
dataset = cdat, burnIn = 2000)
```
```## Generating the Markov chain ...
```
```## Progress:   0%  Progress:   5%  Progress:  10%  Progress:  15%  Progress:  20%  Progress:  25%  Progress:  30%  Progress:  35%  Progress:  40%  Progress:  45%  Progress:  50%  Progress:  55%  Progress:  60%  Progress:  65%  Progress:  70%  Progress:  75%  Progress:  80%  Progress:  85%  Progress:  90%  Progress:  95%  Progress: 100%
```
```## Generation of the Markov Chain required 2 secs
## Conducting inference ...
## Inference required 0 secs
```
```## summary of model including estimates and CIs
summary(m.female)
```
```##
## Call:
## [[1]]
## elrm(formula = admit/ntrials ~ female, interest = ~female, iter = 22000,
##     dataset = cdat, burnIn = 2000)
##
##
## Results:
##        estimate p-value p-value_se mc_size
## female  0.56478   0.484    0.00284   20000
##
## 95% Confidence Intervals for Parameters
##
##            lower   upper
## female -1.129697 2.59968
```

```## trace plot and histogram of sampled values from the sufficient
## statistic
plot(m.female\$mc, col = "grey")
```

```## model with apcalc predictor only
m.apcalc <- elrm(formula = admit/ntrials ~ apcalc, interest = ~apcalc, iter = 22000,
dataset = cdat, burnIn = 2000)
```
```## Generating the Markov chain ...
```
```## Progress:   0%  Progress:   5%  Progress:  10%  Progress:  15%  Progress:  20%  Progress:  25%  Progress:  30%  Progress:  35%  Progress:  40%  Progress:  45%  Progress:  50%  Progress:  55%  Progress:  60%  Progress:  65%  Progress:  70%  Progress:  75%  Progress:  80%  Progress:  85%  Progress:  90%  Progress:  95%  Progress: 100%
```
```## Generation of the Markov Chain required 3 secs
## Conducting inference ...
```
```## Warning: 'apcalc' observed value of the sufficient statistic was not
## sampled
```
```## Inference required 0 secs
```

```## summary of model including estimates and CIs
summary(m.apcalc)
```
```##
## Call:
## [[1]]
## elrm(formula = admit/ntrials ~ apcalc, interest = ~apcalc, iter = 22000,
##     dataset = cdat, burnIn = 2000)
##
##
## Results:
##        estimate p-value p-value_se mc_size
## apcalc       NA       0          0   20000
##
## 95% Confidence Intervals for Parameters
##
##        lower upper
## apcalc    NA    NA
```

```## trace plot and histogram of sampled values from the sufficient
## statistic
plot(m.apcalc\$mc, col = "grey")
```

```## run not automated for time purposes
results <- elrm(formula = admit/ntrials ~ female + apcalc,
interest = ~ female + apcalc,
iter = 5005000, dataset = cdat, burnIn = 5000, r = 2)
```
```## Generating the Markov chain ...
```
```## Progress:   0%  Progress:   5%  Progress:  10%  Progress:  15%  Progress:  20%  Progress:  25%  Progress:  30%  Progress:  35%  Progress:  40%  Progress:  45%  Progress:  50%  Progress:  55%  Progress:  60%  Progress:  65%  Progress:  70%  Progress:  75%  Progress:  80%  Progress:  85%  Progress:  90%  Progress:  95%  Progress: 100%
```
```## Generation of the Markov Chain required 1.4 mins
## Conducting inference ...
## Inference required 1.8833 mins
```
```summary(results)
```
```##
## Call:
## [[1]]
## elrm(formula = admit/ntrials ~ female + apcalc, interest = ~female +
##     apcalc, r = 2, iter = 5005000, dataset = cdat, burnIn = 5000)
##
##
## Results:
##        estimate p-value p-value_se mc_size
## joint        NA 0.00047    0.00001 5000000
## female  1.27347 0.34629    0.01267    1646
## apcalc  3.95707 0.00031    0.00002 1116420
##
## 95% Confidence Intervals for Parameters
##
##            lower    upper
## female -1.157186 5.201046
## apcalc  1.107982 8.762402
```

• The first part of the output echoes the call used to create the model. The results has three rows, one for the joint distribution, and one for each parameter.
• The estimate column has the parameter estimates (joint is NA). The next column gives the p-values. For the joint row, the p-value is that both effects are simultaneously zero. For `female` and `apcalc`, it is the p-value for testing that the individual parameter estimate is zero. Next is the Monte Carlo standard errors for the p-value. Finally, mc_size is the length of the Markov chain of sampled values of sufficient statistics used for each parameter estimate. Note that all the length of the chain for the joint test was 5 million, a mere 1,739 for the female parameter. Because each Markov chain needs to be sufficiently long for stable inference, we needed to increase the total chain size to stably estimate female. For apcalc, the chain length is larger than necessary (over 1 million).
• Below the main results are the confidence intervals for each parameter estimate.
• We can also plot each model, to obtain trace plots of the Markov chains.

## Things to consider

• Firth logit may be helpful if you have separation in your data. This can be done in R using the logistf package.
• Exact logistic regression is an alternative to conditional logistic regression if you have stratification, since both condition on the number of positive outcomes within each stratum. The estimates from these two analyses will be different because a conditional logistic conditions only on the intercept term, exact logistic conditions on the sufficient statistics of the other regression parameters as well as the intercept term.

## References

• UCLA: IDRE (Institute for Digital Research and Education). Data Analysis Examples. from http://www.ats.ucla.edu/stat/dae/ (accessed January 31, 2014)
• Zamar, D., McNeney, B., and Graham, J. (2007). elrm: Software Implementing Exact-like Inference for Logistic Regression Models. http://www.jstatsoft.org/v21/i03/paper
• Collett, D. (2002). Modeling Binary Data, Second Edition. Boca Raton: Chapman and Hall.
• Cox, D. R. and Snell, E. J. (1989). Analysis of Binary Data, Second Edition. Boca Raton: Chapman and Hall.
• Hirji, K. F. (2005). Exact Analysis of Discrete Data. Boca Raton: Chapman and Hall.