Variations on Logistic Regression (ELMR Chapter 4)
Biostat 200C
Author
Dr. Jin Zhou @ UCLA
Published
April 14, 2026
Display system information and load tidyverse and faraway packages
sessionInfo()
R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.7.4
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] htmlwidgets_1.6.4 compiler_4.5.2 fastmap_1.2.0 cli_3.6.5
[5] tools_4.5.2 htmltools_0.5.8.1 rstudioapi_0.17.1 yaml_2.3.10
[9] rmarkdown_2.29 knitr_1.50 jsonlite_2.0.0 xfun_0.56
[13] digest_0.6.37 rlang_1.1.7 evaluate_1.0.5
library(tidyverse)library(faraway)
faraway package contains the datasets in the ELMR book.
1 Latent variables
Consider a latent variable \(T\) and let \[
Y = \begin{cases}
1 & \text{if } T \le t \\
0 & \text{if } T > t
\end{cases}.
\] Then \[
p = \mathbb{P}(Y = 1) = \mathbb{P}(T \le t).
\]
If \(T\) follows the logistic distribution\[
\mathbb{P}(T \le t) = \frac{e^{(t - \mu) / \sigma}}{1 + e^{(t-\mu)/\sigma}},
\] then we recover the logistic regression by setting \[
\frac{t-\mu}{\sigma} = \eta = \mathbf{x}^T \boldsymbol{\beta}.
\]
tibble(t =seq(-6, 6, 0.1), pdf =dlogis(t, location =0, scale =1)) %>%ggplot() +geom_line(mapping =aes(x = t, y = pdf)) +labs(x ="t", y ="Density", title ="Logistic(0, 1) distribution")
2 Link functions
The logit link function\[
\eta = g(p) = \log \frac{p}{1-p}
\] is not the only choice for Bernoulli and binomial regression.
Any functions \(g: [0, 1] \mapsto \mathbb{R}\) that is smooth and monotone qualifies as a link function.
Probit link function, corresponding a normal latent variable. \[
\eta = g(p) = \Phi^{-1}(p),
\] where \(\Phi\) is the cumulative distribution function (cdf) of a standard normal. The corresponding inverse link function is \[
p = \Phi(\eta).
\]
Complementary log-log link functionm, corresponding to a Gumbel-distributed latent variable. \[
\eta = g(p) = \log ( - \log(1-p)).
\] The corresponding inverse link function is \[
p = 1 - e^{-e^{\eta}}.
\]
Cauchit link function, corresponding to a Cauchy-distributed latent variable. \[
\eta = g(p) = \tan((p - 1/2) \pi).
\] The corresponding inverse link function is \[
p = \frac{\arctan(\eta)}{\pi} + \frac 12.
\]
3 Bliss data
bliss data records the number of insects dying at different levels of insecticide concentration.
tibble(conc = bliss$conc,logit =predict(mlogit , type ="response"),probit =predict(mprobit , type ="response"),cloglog =predict(mcloglog, type ="response"),cauchit =predict(mcauchit, type ="response"))
we don’t see vast difference in the predictions from 4 link functions.
Let’s compare the predictions at a wider range [-4, 8]. We observe wider differences at extreme values of linear predictors. Logit and probit are very close. Cauchit approaches 0 and 1 at slower rate, while cloglog at a faster rate.
df <-tibble(conc =seq(-4, 8, 0.2))tibble(dose = df$conc,logit =predict(mlogit , type ="response", newdata = df),probit =predict(mprobit , type ="response", newdata = df),cloglog =predict(mcloglog, type ="response", newdata = df),cauchit =predict(mcauchit, type ="response", newdata = df)) %>%pivot_longer(logit:cauchit, names_to ="link", values_to ="pred") %>%ggplot() +geom_line(mapping =aes(x = dose, y = pred, color = link)) +labs(x ="Dose", y ="Predicted probability", title ="Predicions from different link functions")
Logit link is the default choice because of
simpler math,
easier interpretation using odds, and
and easier analysis of retrospectively sampled data.
4 Prospective and retrospective sampling
In prospective sampling or cohort study or follow-up study, the predictors are fixed and then the outcome is observed.
In retrospective sampling or case-control study, the outcome is fixed and then the predictors are observed. It is required that the probability of inclusion in the study is independent of the predictor values.
Baby food example. We want to study the effect of sex and feeding method on whether a baby gets respiratory disease in the first year.
babyfood
disease nondisease sex food
1 77 381 Boy Bottle
2 19 128 Boy Suppl
3 47 447 Boy Breast
4 48 336 Girl Bottle
5 16 111 Girl Suppl
6 31 433 Girl Breast
xtabs(disease / (disease + nondisease) ~ sex + food, babyfood)
food
sex Bottle Breast Suppl
Boy 0.16812227 0.09514170 0.12925170
Girl 0.12500000 0.06681034 0.12598425
Assuming prospective sampling (collect babies first then follow up whether they develop respiratory disease in the 1 year), we fit a logistic regression
library(gtsummary)babyfood %>%# change reference level to Breast and Girlmutate(food =relevel(food, ref ="Breast"), sex =relevel(sex, ref ="Girl")) %>%glm(cbind(disease, nondisease) ~ sex + food, family = binomial, data = .) %>%tbl_regression(intercept =TRUE, exponentiate = F)
Characteristic
log(OR)
95% CI
p-value
(Intercept)
-2.6
-2.9, -2.3
<0.001
sex
Girl
—
—
Boy
0.31
0.04, 0.59
0.027
food
Breast
—
—
Bottle
0.67
0.37, 0.97
<0.001
Suppl
0.50
0.06, 0.91
0.022
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
The regression coefficient 0.67 for Bottle respresents the increased risk of developing respiratory disease incurred by bottle feeding relative to breast feeding. Similarly the regression coefficient 0.31 represents the increased risk of boys relative to girls.
Assuming retrospective sampling (inspect the medical history of many babies at 1 year old), it seems more sensible compute the log-odds of sex given the respiratory disease status. We actually observe the same regression coefficient 0.31 for dis!
babyfood %>%# change reference level to Breast and Girlmutate(food =relevel(food, ref ="Breast"), sex =relevel(sex, ref ="Girl")) %>%pivot_longer(disease:nondisease, names_to ="dis", values_to ="count") %>%mutate(dis =relevel(as.factor(dis), ref ="nondisease")) %>%glm(sex ~ dis + food, family = binomial, weights = count, data = .) %>%tbl_regression(intercept =TRUE, exponentiate = F)
Characteristic
log(OR)
95% CI
p-value
(Intercept)
0.04
-0.09, 0.17
0.6
dis
nondisease
—
—
disease
0.31
0.04, 0.59
0.027
food
Breast
—
—
Bottle
0.09
-0.09, 0.28
0.3
Suppl
0.07
-0.20, 0.34
0.6
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
We see that using logistic regression for the retrospective design is as effective as a prospective design for estimating a relative risk. Retrospective design is more convenient because
it is easier to collect cases and controls (we don’t need to wait for long time to even not able to collect any rare cases),
is is simpler to collect many predictors.
Prospective design has its unique advantage though. Let \[\begin{eqnarray*}
\pi_0 &=& \text{probability of being included in the study if they do not have disease} \\
\pi_1 &=& \text{probability of being included in the study if they do have disease}.
\end{eqnarray*}\] In prosective studies, \(\pi_0 = \pi_1\). In retrospective studies, typically \(\pi_1 \gg \pi_0\). Let \[\begin{eqnarray*}
p^\star(\mathbf{x}) &=& \text{conditional probability that an individual has the disease given that he or she was included in the study} \\
p(\mathbf{x}) &=& \text{unconditional probability that an individual has the disease as we would obtain from a prospective study}.
\end{eqnarray*}\] By Bayes theorem \[
p^\star(\mathbf{x}) = \frac{\pi_1 p(\mathbf{x})}{\pi_1 p(\mathbf{x}) + \pi_0 (1 - p(\mathbf{x}))},
\] which after rearrangement shows \[
\text{logit}(p^\star(\mathbf{x})) = \log \frac{\pi_1}{\pi_0} + \text{logit}(p(\mathbf{x})).
\] In general \(\pi_1/\pi_0\) is unknown, therefore we would not be able to estimate intercept \(\beta_0\) in a retrospective study. However we can still estimate the relative effect, e.g., odds ratio, of non-intercept predictors.
5 Prediction
Given covariate \(\mathbf{x}_0\), the predicted response on the link scale is \[
\widehat{\eta} = \mathbf{x}_0^T \widehat{\boldsymbol{\beta}}
\] with variance \[
\operatorname{Var} \widehat{\eta} = \mathbf{x}_0^T \cdot \operatorname{Var} \widehat{\boldsymbol{\beta}} \cdot \mathbf{x}_0 = \mathbf{x}_0^T \left( \mathbf{X}^T \widehat{\mathbf{W}} \mathbf{X} \right)^{-1} \mathbf{x}_0.
\]
For the insecticide data bliss, to predict the probability of killing an insect at dose 2.5:
# prediction at linear predictor scalepredict(mlogit, newdata =tibble(conc =2.5), type ="link", se.fit =TRUE)
and a 95% confidence interval on the probability scale is
ilogit(c(eta0 -1.96* se, eta0 +1.96* se))
[1] 0.5342962 0.7358471
6 Effective doses, ED50 and LD50
When there is a single continuous predictor or when other predictors are held fixed, we may wish to estimate the value of \(x\) corresponding to a chosen \(p\).
ED50 stands for the effective dose for which there will be a 50% chance of “success”. When a “success” is to kill the subjects or determine toxicity, the term LD50 (lethal dose) would be used.
Solving \[
p = \frac{e^{\beta_0 + x \beta_1}}{1 + e^{\beta_0 + x \beta_1}} = 0.5
\] yields \[
\widehat{\text{ED}_{50}} = - \frac{\widehat{\beta}_0}{\widehat{\beta}_1}.
\]
In the bliss data,
(ld50 <-- mlogit$coef[1] / mlogit$coef[2])
(Intercept)
2
To quantify the uncertainty in estimating ED50 or LD50, we can use the delta method\[
\operatorname{Var} g\left(\widehat{\boldsymbol{\theta}}\right) \approx \nabla g\left(\widehat{\boldsymbol{\theta}}\right)^T \cdot \operatorname{Var} \widehat{\boldsymbol{\theta}} \cdot \nabla g\left(\widehat{\boldsymbol{\theta}}\right).
\] For ED50 or LD50 \[
\nabla g(\beta_0, \beta_1) = \begin{pmatrix} - \frac{1}{\beta_1} \\ \frac{\beta_0}{\beta_1^2} \end{pmatrix}.
\]
The MASS package has a convenience function dose.p for calculating the effective doses and their standard errors
library(MASS)dose.p(mlogit, p =c(0.25, 0.5, 0.75))
Dose SE
p = 0.25: 1.054465 0.2315932
p = 0.50: 2.000000 0.1784367
p = 0.75: 2.945535 0.2315932
7 Conditional logistic regression for matched case-control studies
In case-control studies, we try to determine the effect of certain risk factors on the outcome. Ideally we hope to collect all confouding variables and model them in the correct way in the logistic regression. In reality, it can be difficult.
In a matched case-control study, we match each case with one or more controls that have the same or similar values of some set of potential confouding variables. For example, if we have a 56-year-old, Hispanic male case, we try to match him with a few controls who are also 56-year-old Hispanic males. Matching also gvies us the possibility of adjusting for confounders that are difficult to measure, e.g., diet, environmental exposure, etc.
Disadvantages of matched case-control study include
it can be difficult to form the matched sets,
one cannot estimate the effects of the variables used to determine the matches,
it may be difficult to generalize to the population.
In a 1:M design, we match \(M\) controls to each case. Suppose we have \(n\) matched sets and we take \(i=0\) to represent the case and \(i=1,\ldots,M\) to represent the controls. We propose a logistic regression \[
\text{logit}(p_j(\mathbf{x}_{ij})) = \alpha_j + \mathbf{x}_{ij}^T \boldsymbol{\beta},
\] where \(\alpha_j\) models the effect of the confounding variables in the \(j\)-th matched set. Thus \[
p_j(\mathbf{x}_{ij}) = \frac{e^{\alpha_j + \mathbf{x}_{ij}^T \boldsymbol{\beta}}}{1 + e^{\alpha_j + \mathbf{x}_{ij}^T \boldsymbol{\beta}}}, \quad i = 0, 1, \ldots, M.
\]
Given a matched set \(j\) of \(M+1\) subjects known to have one case and \(M\) controls, the conditional probability of the observed outcome, or, in other words, that subject \(i=0\) is the case and the rest are controls is \[\begin{eqnarray*}
& & \mathbb{P} \left( Y_{0j}=1, Y_{1j}=\cdots = Y_{Mj}=0 \mid \sum_{i=0}^M Y_{ij} = 1 \right) \\
&=& \frac{p_j(\mathbf{x}_{0j}) \prod_{i=1}^M (1 - p_j(\mathbf{x}_{ij}))}{\sum_{i=0}^M p_j(\mathbf{x}_{ij}) \prod_{i'\ne i} (1 - p_j(\mathbf{x}_{i'j}))} \\
&=& \frac{\exp \mathbf{x}_{0j}^T \boldsymbol{\beta}}{\sum_{i=0}^M \exp \mathbf{x}_{ij}^T \boldsymbol{\beta}} \\
&=& \frac{1}{1 + \sum_{i=1}^M \exp (\mathbf{x}_{ij} - \mathbf{x}_{0j})^T \boldsymbol{\beta}}.
\end{eqnarray*}\]\(\alpha_j\) conveniently cancels out in the final expression. We can form the conditional likelihood function \[
L(\boldsymbol{\beta}) = \prod_{j=1}^n \frac{1}{1 + \sum_{i=1}^M \exp (\mathbf{x}_{ij} - \mathbf{x}_{0j})^T \boldsymbol{\beta}},
\] which is identical to that from a Cox proportional hazards mdoel.
The amlxray data concerns the x-ray exposure and childhood acute myeloid leukemia. The sets are matched on age, race, and county of residence.
(amlxray <-as_tibble(amlxray))
# A tibble: 238 × 11
ID disease Sex downs age Mray MupRay MlowRay Fray Cray CnRay
<fct> <dbl> <fct> <fct> <int> <fct> <fct> <fct> <fct> <fct> <ord>
1 7004 1 F no 0 no no no no no 1
2 7004 0 F no 0 no no no no no 1
3 7006 1 M no 6 no no no no yes 3
4 7006 0 M no 6 no no no no yes 2
5 7009 1 F no 8 no no no no no 1
6 7009 0 F no 8 no no no no no 1
7 7010 1 M yes 1 no no no no no 1
8 7010 0 M no 1 no no no no no 1
9 7013 1 M no 4 yes no yes no no 1
10 7013 0 M no 4 no no no no yes 2
# ℹ 228 more rows
Downs syndrome is a known risk factor. In this data set, all Downs syndrome children are cases. The coefficient for downs will be infinity.
# A tibble: 7 × 11
ID disease Sex downs age Mray MupRay MlowRay Fray Cray CnRay
<fct> <dbl> <fct> <fct> <int> <fct> <fct> <fct> <fct> <fct> <ord>
1 7010 1 M yes 1 no no no no no 1
2 7018 1 F yes 0 no no no no no 1
3 7066 1 F yes 1 no no no no no 1
4 7077 1 M yes 1 yes yes no no no 1
5 7146 1 F yes 1 yes no yes yes no 1
6 7176 1 F yes 1 no no no no no 1
7 7189 1 F yes 1 no no no no no 1
So we exclude Donws syndrome cases and their matched subjects.
Now we can fit a conditional logit model using predictors Sex, Mray (mother exposure to x-ray), Fray (father exposure to x-ray), and CnRay (child exposure to x-ray). Since CnRray is an ordered factor, linear, quadratic, and cubic contrasts (from orthogonal polynomial coding) are estimated. Only the linear effect is significant.
library(survival)cmod <-clogit(disease ~ Sex + Mray + Fray + CnRay +strata(ID), data = ramlxray)summary(cmod)
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
8 Separation (ELMR 2.7)
Let’s start with the famours Fisher Iris data. We restrict to species setosa and versicolor and classify them using predictors Sepal.Width and Sepal.Length.
lmod <-glm(Species ~ Sepal.Width + Sepal.Length, family = binomial, data = irisr)
Warning: glm.fit: algorithm did not converge
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(lmod)
Call:
glm(formula = Species ~ Sepal.Width + Sepal.Length, family = binomial,
data = irisr)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -360.6 195972.9 -0.002 0.999
Sepal.Width -110.1 55361.5 -0.002 0.998
Sepal.Length 131.8 64577.0 0.002 0.998
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1.3863e+02 on 99 degrees of freedom
Residual deviance: 7.1185e-09 on 97 degrees of freedom
AIC: 6
Number of Fisher Scoring iterations: 25
glm function is complaining there is some issue in algorithmic convergence. The residual deviance is essential 0, indicating a perfect fit, but none of the predictors are significant.
Graphical summary of the data reveals that the two species setosa and versicolor are perfectly separable by a straight line \(\eta = \beta_0 + x_{\text{Sepal.Width}} \beta_{\text{Sepal.Width}} + x_{\text{Sepal.Length}} \beta_{\text{Sepal.Length}} = 0\).
irisr %>%ggplot() +geom_point(mapping =aes(x = Sepal.Width, y = Sepal.Length, color = Species)) +geom_abline(intercept =- lmod$coef[1] / lmod$coef[3], slope =- lmod$coef[2] / lmod$coef[3]) +labs(title ="Perfect separation in Iris data")
For correct inference, we can use the bias reduction method implemented in brglm package. It gives finite estimates and correct inference on the coefficients.
library(brglm)bmod <-brglm(Species ~ Sepal.Width + Sepal.Length, family = binomial, data = irisr)summary(bmod)
Call:
brglm(formula = Species ~ Sepal.Width + Sepal.Length, family = binomial,
data = irisr)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -24.508 12.493 -1.962 0.04979 *
Sepal.Width -8.897 2.748 -3.237 0.00121 **
Sepal.Length 9.732 3.334 2.919 0.00351 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 130.638 on 99 degrees of freedom
Residual deviance: 3.323 on 97 degrees of freedom
Penalized deviance: 6.60971
AIC: 9.323
