Logit, ologit and clogit: One Link Function, Three Different Questions
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
Logistic, ordinal and conditional logistic regression share the logit link, so exponentiated coefficients are odds ratios. Outcome type and sampling design choose the model. In 500 simulated newborns, inadequate antenatal care has an odds ratio of 2.01 for a higher complication grade. It is adjusted for preterm birth, a possible mediator, so 2.01 leaves out any effect that runs through preterm birth. The ordinal model assumes one odds ratio at every cut between grades. A Brant test finds evidence of a departure for preterm birth (P = 0.003) but none for antenatal care. In 20 matched case-control pairs, ordinary logistic regression gives an odds ratio of 1.83, conditional logistic regression 4.00. They answer different questions, pooled versus within pairs; the design decides which to report. Only 5 pairs differ in exposure, so 4.00 is imprecise (95% confidence interval 0.45 to 35.79). This article works both estimates by hand, then gives Stata and R code and a decision table.
Two designs, three commands, two odds ratios
A neonatal unit wants to know whether inadequate antenatal care (ANC) goes with more complications in newborns. Inadequate here means too few or too late check-ups during pregnancy. The team studies the question with two designs at once. The unit, its newborns and every number below are simulated.
In a birth cohort, each of 500 newborns has a complication grade within 28 days: none, mild, moderate or severe. The grades are ordered, but the steps between them need not be equal. Separately, a case-control study matches each of 20 newborns with a complication (cases) to one control. Each control was a newborn still free of complications on the day its case had one.
The analyst runs three Stata commands: logistic, ologit and clogit. On the matched pairs, logistic gives an odds ratio of 1.83 and clogit gives 4.00. The newborns, the exposure and the outcome are the same, so a reviewer asks which number is right.
Each number is correct arithmetic, but the two make different comparisons, and only the matched analysis follows the design. All three models share one piece, the logit link. They differ in what they assume about the outcome and about how the people in the data were sampled. Before each model, this article names the outcome type and the sampling design, because those two facts choose the command.
The logit link that all three models share
Every model in this article describes a probability $p$, such as the probability that a newborn has a complication. A probability is confined between 0 and 1, while a regression equation can take any value. A link function joins the two: it transforms $p$ onto a scale on which a linear model can be fitted.
The logit link is the log of the odds:
$$\operatorname{logit}(p) = \ln\left(\frac{p}{1 - p}\right)$$Here $p/(1 - p)$ is the odds and $\ln$ is the natural logarithm. The logit runs from minus infinity to plus infinity. A sum of coefficients times covariates, the linear predictor, can take any value and still map to a valid probability.
With one binary exposure $x$, coded 1 for inadequate ANC and 0 for adequate ANC, ordinary logistic regression fits
$$\operatorname{logit} P(Y = 1 \mid x) = \beta_0 + \beta_1 x$$Here $Y = 1$ marks a newborn with the outcome, and $\beta_0$ is the log odds when $x = 0$. The slope $\beta_1$ is the change in log odds when $x$ goes from 0 to 1. A difference of two log odds is the log of their ratio. So $\beta_1$ is a log odds ratio, and $e^{\beta_1}$ is the odds ratio [1].
Stata's logit command prints $\beta_1$; its logistic command fits the same model and prints $e^{\beta_1}$. R's glm() with family = binomial uses the logit link by default and prints $\beta_1$.
Ordinary logistic regression assumes two things about the data. The outcome is binary, and observations are independent: one person's outcome says nothing about another's, given the covariates. The ordered grades break the first assumption, and the matched pairs break the second. Matching also ties each control to its case, so it has to enter the comparison itself, not only the standard error.
Name the outcome type and the sampling design first
The three analyses in this article share the logit link but call for different models. The choice rests on two facts known before any data are analysed. The first is the outcome type: two categories, or several ordered ones. The second is the sampling design: independent people, or matched sets made of a case and its controls.
- Binary outcome, independent people: ordinary logistic regression; Stata
logistic, Rglm(..., family = binomial). - Ordered categories, independent people: ordinal logistic regression; Stata
ologit, RMASS::polr(). - Binary outcome (case or control), matched sets: conditional logistic regression; Stata
clogit, Rsurvival::clogit().
The logit appears in a slightly different place in each model. Logistic regression models the logit of one probability. Ordinal logistic regression models one logit per cut between grades. Conditional logistic regression models the logit with a separate intercept for each matched set, then removes those intercepts.
Each model section below opens by naming these two facts. The decision table near the end sets the three models side by side, with what each exponentiated coefficient means.
The simulated birth cohort
The cohort holds 500 newborns: 200 whose mothers had inadequate ANC and 120 born preterm, with some in both groups. Each newborn has a grade $Y$ within 28 days: 0 none, 1 mild, 2 moderate, 3 severe. In all, 256 newborns have no complication, 101 a mild one, 71 a moderate one and 72 a severe one.
The counts were built from a known model, so every estimate can be checked against it. For each cut between grades, the model sets the log odds of being at grade $j$ or worse:
$$\operatorname{logit} P(Y \ge j) = \alpha_j + \beta x + \gamma_j z$$Here $x$ is 1 for inadequate ANC, $z$ is 1 for preterm birth, and $j = 1, 2, 3$ indexes the three cuts. The intercepts are $\alpha_1 = -0.4$, $\alpha_2 = -1.4$ and $\alpha_3 = -2.5$. The ANC slope is $\beta = \ln 2$, an odds ratio of 2.00 at every cut. The preterm slopes are $\gamma_1 = 0.3$, $\gamma_2 = 0.7$ and $\gamma_3 = 1.2$: odds ratios of 1.35, 2.01 and 3.32, growing with severity.
This simulated truth gives ANC one odds ratio at every cut and preterm birth a different one at each. The counts are the expected numbers under the model, rounded to whole newborns, not random draws. So the fitted odds ratios below land close to the values used to build the data.
Stata: building the cohort's counts
* expected number at or above each cut, rounded to whole newborns
forvalues j = 1/3 {
generate int atleast`j' = floor(newborns*invlogit(`alpha`j'' + ln(2)*inadequate_anc + `gamma`j''*preterm) + 0.5)
}
generate int sev0 = newborns - atleast1
generate int sev1 = atleast1 - atleast2
generate int sev2 = atleast2 - atleast3
generate int sev3 = atleast3
list inadequate_anc preterm newborns sev0 sev1 sev2 sev3, noobs. * expected number at or above each cut, rounded to whole newborns
. forvalues j = 1/3 {
2. generate int atleast`j' = floor(newborns*invlogit(`alpha`j'' + ln(2)*inadequate_anc + `gamma`j''*preterm) + 0.5)
3. }
. generate int sev0 = newborns - atleast1
. generate int sev1 = atleast1 - atleast2
. generate int sev2 = atleast2 - atleast3
. generate int sev3 = atleast3
. list inadequate_anc preterm newborns sev0 sev1 sev2 sev3, noobs
+-----------------------------------------------------------+
| inadeq~c preterm newborns sev0 sev1 sev2 sev3 |
|-----------------------------------------------------------|
| 0 0 240 144 49 29 18 |
| 0 1 60 31 9 7 13 |
| 1 0 140 60 34 26 20 |
| 1 1 60 21 9 9 21 |
+-----------------------------------------------------------+First model: a binary outcome in independent newborns
Outcome type: binary, any complication (grade 1 or worse) against none. Sampling design: independent newborns from one birth cohort. Model: ordinary logistic regression, with inadequate ANC and preterm birth as covariates. Preterm birth may itself lie on the path from ANC to complications, a point taken up near the end.
Of the 500 newborns, 244 had a complication of some grade. The fitted ANC coefficient is 0.691 on the log odds scale. Exponentiated, it gives an odds ratio of 2.00 (95% confidence interval 1.38 to 2.88), adjusted for preterm birth. The preterm odds ratio is 1.40 (0.92 to 2.13).
Both odds ratios describe one cut only, none against any complication. The model cannot tell a mild complication from a severe one. Collapsing four grades into two discards that information, and the next model keeps it.
Stata: logit and logistic on any complication
generate byte any_morbidity = severity >= 1
quietly count if any_morbidity == 1
canon bin.n_any 0 r(N)
logit any_morbidity i.inadequate_anc i.preterm, nolog
canon logit.b_anc 4 _b[1.inadequate_anc]
canon logit.se_anc 4 _se[1.inadequate_anc]
canon logit.b_preterm 4 _b[1.preterm]
canon logit.se_preterm 4 _se[1.preterm]
canon logit.b_cons 4 _b[_cons]
canon logit.ll 4 e(ll)
logistic any_morbidity i.inadequate_anc i.preterm, nolog. generate byte any_morbidity = severity >= 1
. quietly count if any_morbidity == 1
. canon bin.n_any 0 r(N)
CANON c7.bin.n_any 244
. logit any_morbidity i.inadequate_anc i.preterm, nolog
Logistic regression Number of obs = 500
LR chi2(2) = 17.79
Prob > chi2 = 0.0001
Log likelihood = -337.53237 Pseudo R2 = 0.0257
----------------------------------------------------------------------------------
any_morbidity | Coefficient Std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc | .6912976 .1868748 3.70 0.000 .3250298 1.057565
1.preterm | .3354302 .2149308 1.56 0.119 -.0858265 .756687
_cons | -.404775 .1254744 -3.23 0.001 -.6507003 -.1588497
----------------------------------------------------------------------------------
. canon logit.b_anc 4 _b[1.inadequate_anc]
CANON c7.logit.b_anc 0.6913
. canon logit.se_anc 4 _se[1.inadequate_anc]
CANON c7.logit.se_anc 0.1869
. canon logit.b_preterm 4 _b[1.preterm]
CANON c7.logit.b_preterm 0.3354
. canon logit.se_preterm 4 _se[1.preterm]
CANON c7.logit.se_preterm 0.2149
. canon logit.b_cons 4 _b[_cons]
CANON c7.logit.b_cons -0.4048
. canon logit.ll 4 e(ll)
CANON c7.logit.ll -337.5324
. logistic any_morbidity i.inadequate_anc i.preterm, nolog
Logistic regression Number of obs = 500
LR chi2(2) = 17.79
Prob > chi2 = 0.0001
Log likelihood = -337.53237 Pseudo R2 = 0.0257
----------------------------------------------------------------------------------
any_morbidity | Odds ratio Std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc | 1.996304 .3730589 3.70 0.000 1.384072 2.879352
1.preterm | 1.398542 .3005898 1.56 0.119 .9177535 2.131204
_cons | .6671269 .0837074 -3.23 0.001 .5216803 .8531246
----------------------------------------------------------------------------------
Note: _cons estimates baseline odds.R: glm() on any complication
o$any_morbidity <- as.integer(o$severity >= 1)
canon("bin.n_any", sum(o$any_morbidity), 0)
fl <- glm(any_morbidity ~ inadequate_anc + preterm, family = binomial, data = o)
print(summary(fl)$coefficients)
se_l <- sqrt(diag(vcov(fl)))
canon("logit.b_anc", coef(fl)[["inadequate_anc"]])
canon("logit.se_anc", se_l[["inadequate_anc"]])
canon("logit.b_preterm", coef(fl)[["preterm"]])
canon("logit.se_preterm", se_l[["preterm"]])
canon("logit.b_cons", coef(fl)[["(Intercept)"]])
canon("logit.ll", as.numeric(logLik(fl)))
# the same fit shown as odds ratios
print(exp(cbind(OR = coef(fl), confint.default(fl)))) > o$any_morbidity <- as.integer(o$severity >= 1)
> canon("bin.n_any", sum(o$any_morbidity), 0)
CANON c7.bin.n_any 244
> fl <- glm(any_morbidity ~ inadequate_anc + preterm,
+ family = binomial, data = o)
> print(summary(fl)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4047750 0.1254744 -3.225956 0.0012555244
inadequate_anc 0.6912976 0.1868748 3.699256 0.0002162321
preterm 0.3354302 0.2149308 1.560643 0.1186080947
> se_l <- sqrt(diag(vcov(fl)))
> canon("logit.b_anc", coef(fl)[["inadequate_anc"]])
CANON c7.logit.b_anc 0.6913
> canon("logit.se_anc", se_l[["inadequate_anc"]])
CANON c7.logit.se_anc 0.1869
> canon("logit.b_preterm", coef(fl)[["preterm"]])
CANON c7.logit.b_preterm 0.3354
> canon("logit.se_preterm", se_l[["preterm"]])
CANON c7.logit.se_preterm 0.2149
> canon("logit.b_cons", coef(fl)[["(Intercept)"]])
CANON c7.logit.b_cons -0.4048
> canon("logit.ll", as.numeric(logLik(fl)))
CANON c7.logit.ll -337.5324
> print(exp(cbind(OR = coef(fl), confint.default(fl))))
OR 2.5 % 97.5 %
(Intercept) 0.6671269 0.5216803 0.8531246
inadequate_anc 1.9963043 1.3840720 2.8793525
preterm 1.3985420 0.9177535 2.1312037 Second model: ordered grades in the same newborns
Outcome type: ordered categories, the four complication grades. Sampling design: the same independent newborns. Model: ordinal logistic regression, which uses all the grades at once and keeps their order [2].
An outcome with four grades has three cut points, the boundaries between neighbouring grades. Each cut splits the newborns into two groups: grade $j$ or worse, against better than grade $j$. The log odds of being at grade $j$ or worse is the cumulative logit for that cut.
Ordinal logistic regression fits one cumulative logit per cut, each with its own intercept, but with one shared set of slopes:
$$\operatorname{logit} P(Y \ge j) = \alpha_j + \beta x + \gamma z$$The intercepts $\alpha_j$ fall as $j$ rises, because fewer newborns reach a higher grade. The slopes $\beta$ and $\gamma$ carry no $j$, so the ANC odds ratio $e^{\beta}$ is the same at every cut. This is the proportional odds assumption. Unlike the model that built the data, it forces preterm birth to have one slope at every cut.
Stata's ologit and R's polr() print the intercepts with the opposite sign, as cut points $\kappa_j = -\alpha_j$ that rise with $j$. The slopes keep their sign: a positive coefficient makes higher grades more likely.
In the simulated cohort, the ANC odds ratio adjusted for preterm birth is 2.01 (95% confidence interval 1.43 to 2.82). That is close to the 2.00 used to build the data. It is the odds ratio for being in a higher grade, wherever the line between grades is drawn. Its interval is narrower than the binary model's 1.38 to 2.88, because the ordinal model also uses the grades above the first cut.
Read concretely, the odds ratio applies to every split of the grades at once. Among newborns of the same preterm status, inadequate ANC about doubles the odds of any complication. It does the same for the odds of a moderate or severe complication, and of a severe one. Holding preterm birth fixed may omit part of ANC's effect, as What the three odds ratios have in common explains.
With a common outcome, doubling the odds does not double the risk, as the probabilities below show. The cut points are 0.458, 1.361 and 2.267. Take a newborn born at term, not preterm, with adequate ANC. The model gives this newborn a probability of 0.613 of no complication and 0.094 of a severe one.
With inadequate ANC, these become 0.441 and 0.172. One odds ratio thus shifts the whole distribution of grades toward the severe end. That single shift is what the proportional odds assumption buys, and what it risks.
The preterm odds ratio is 1.84 (1.24 to 2.74). Unlike the ANC odds ratio, it summarises odds ratios that, in the simulated truth, rise from 1.35 at the first cut to 3.32 at the third. Whether one number is a fair summary is the question the proportional odds check asks.
Stata: ologit on the four grades
ologit severity i.inadequate_anc i.preterm, nolog. ologit severity i.inadequate_anc i.preterm, nolog
Ordered logistic regression Number of obs = 500
LR chi2(2) = 27.91
Prob > chi2 = 0.0000
Log likelihood = -597.08706 Pseudo R2 = 0.0228
----------------------------------------------------------------------------------
severity | Coefficient Std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc | .6966207 .1732433 4.02 0.000 .3570702 1.036171
1.preterm | .6124173 .2012158 3.04 0.002 .2180416 1.006793
-----------------+----------------------------------------------------------------
/cut1 | .4581421 .1220414 .2189454 .6973388
/cut2 | 1.360579 .1350555 1.095875 1.625283
/cut3 | 2.266819 .1627571 1.947821 2.585817
----------------------------------------------------------------------------------R: polr() on the four grades
fo <- polr(factor(severity, levels = 0:3, ordered = TRUE) ~ inadequate_anc + preterm, data = o,
method = "logistic", Hess = TRUE)
print(summary(fo))> fo <- polr(factor(severity, levels = 0:3, ordered = TRUE) ~
+ inadequate_anc + preterm, data = o, method = "logistic",
+ Hess = TRUE)
> print(summary(fo))
Call:
polr(formula = factor(severity, levels = 0:3, ordered = TRUE) ~
inadequate_anc + preterm, data = o, Hess = TRUE, method = "logistic")
Coefficients:
Value Std. Error t value
inadequate_anc 0.6966 0.1732 4.021
preterm 0.6124 0.2012 3.044
Intercepts:
Value Std. Error t value
0|1 0.4581 0.1220 3.7540
1|2 1.3606 0.1351 10.0742
2|3 2.2668 0.1628 13.9276
Residual Deviance: 1194.174
AIC: 1204.174Checking proportional odds one cut at a time
The proportional odds assumption can be looked at directly. Fit an ordinary logistic regression at each of the three cuts. The three fits compare grade 1 or worse with none, grade 2 or worse with less, and grade 3 with less. If the assumption holds, each covariate's odds ratio should be similar in all three.
The table below shows that ANC behaves as the assumption requires, with odds ratios of 2.00, 2.01 and 2.01 at the three cuts. Preterm birth does not: 1.40, 2.05 and 3.31, close to the simulated truth of 1.35, 2.01 and 3.32. The first cut's fit is the binary model above, which is why its numbers match that model's.
| Cut | ANC odds ratio, fitted | ANC odds ratio, simulated truth | Preterm odds ratio, fitted | Preterm odds ratio, simulated truth |
|---|---|---|---|---|
| Grade 1 or worse versus none | 2.00 | 2.00 | 1.40 | 1.35 |
| Grade 2 or worse versus grades 0 to 1 | 2.01 | 2.00 | 2.05 | 2.01 |
| Grade 3 versus grades 0 to 2 | 2.01 | 2.00 | 3.31 | 3.32 |
| Ordinal model, one odds ratio for all cuts | 2.01 | 2.00 | 1.84 | No single value |
Stata: one logistic fit per cut
* the proportional-odds check, step 1: a separate binary logit at each cut (severity >= 1, >= 2, >= 3);
* plain 0/1 covariates, so that e(b) holds inadequate_anc, preterm and the constant in that order
forvalues j = 1/3 {
generate byte ge`j' = severity >= `j'
logit ge`j' inadequate_anc preterm, nolog
matrix bcut`j' = e(b)
predict double pi`j', pr
canon cut`j'.b_anc 4 _b[inadequate_anc]
canon cut`j'.or_anc 4 exp(_b[inadequate_anc])
canon cut`j'.b_preterm 4 _b[preterm]
canon cut`j'.or_preterm 4 exp(_b[preterm])
}. * the proportional-odds check, step 1: a separate binary logit at each cut (severity >= 1, >= 2, >= 3);
. * plain 0/1 covariates, so that e(b) holds inadequate_anc, preterm and the constant in that order
. forvalues j = 1/3 {
2. generate byte ge`j' = severity >= `j'
3. logit ge`j' inadequate_anc preterm, nolog
4. matrix bcut`j' = e(b)
5. predict double pi`j', pr
6. canon cut`j'.b_anc 4 _b[inadequate_anc]
7. canon cut`j'.or_anc 4 exp(_b[inadequate_anc])
8. canon cut`j'.b_preterm 4 _b[preterm]
9. canon cut`j'.or_preterm 4 exp(_b[preterm])
10. }
Logistic regression Number of obs = 500
LR chi2(2) = 17.79
Prob > chi2 = 0.0001
Log likelihood = -337.53237 Pseudo R2 = 0.0257
--------------------------------------------------------------------------------
ge1 | Coefficient Std. err. z P>|z| [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc | .6912976 .1868748 3.70 0.000 .3250298 1.057565
preterm | .3354302 .2149308 1.56 0.119 -.0858265 .756687
_cons | -.404775 .1254744 -3.23 0.001 -.6507003 -.1588497
--------------------------------------------------------------------------------
CANON c7.cut1.b_anc 0.6913
CANON c7.cut1.or_anc 1.9963
CANON c7.cut1.b_preterm 0.3354
CANON c7.cut1.or_preterm 1.3985
Logistic regression Number of obs = 500
LR chi2(2) = 24.29
Prob > chi2 = 0.0000
Log likelihood = -287.11966 Pseudo R2 = 0.0406
--------------------------------------------------------------------------------
ge2 | Coefficient Std. err. z P>|z| [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc | .6964978 .2038964 3.42 0.001 .2968681 1.096127
preterm | .7169957 .2238947 3.20 0.001 .2781702 1.155821
_cons | -1.411917 .1515253 -9.32 0.000 -1.708901 -1.114933
--------------------------------------------------------------------------------
CANON c7.cut2.b_anc 0.6965
CANON c7.cut2.or_anc 2.0067
CANON c7.cut2.b_preterm 0.7170
CANON c7.cut2.or_preterm 2.0483
Logistic regression Number of obs = 500
LR chi2(2) = 29.04
Prob > chi2 = 0.0000
Log likelihood = -191.55985 Pseudo R2 = 0.0705
--------------------------------------------------------------------------------
ge3 | Coefficient Std. err. z P>|z| [95% conf. interval]
---------------+----------------------------------------------------------------
inadequate_anc | .6983656 .2648017 2.64 0.008 .1793638 1.217367
preterm | 1.1974 .2683201 4.46 0.000 .6715027 1.723298
_cons | -2.501102 .2175031 -11.50 0.000 -2.9274 -2.074803
--------------------------------------------------------------------------------
CANON c7.cut3.b_anc 0.6984
CANON c7.cut3.or_anc 2.0105
CANON c7.cut3.b_preterm 1.1974
CANON c7.cut3.or_preterm 3.3115The Brant test: equal slopes across cuts?
Comparing three odds ratios by eye is not a test, because each fitted value carries sampling error. The Brant test asks whether the slopes of the separate cut fits differ by more than that error allows [3]. It is a Wald test: it weighs the differences between slopes against their estimated sampling variation and returns a chi-squared statistic.
The separate fits use the same newborns, so their slopes are correlated. Brant's method estimates that correlation from the three fits, where a naive comparison would ignore it [3]. The test can be run for all covariates at once or for one covariate at a time. Its degrees of freedom count the slope differences tested, two per covariate here.
For all covariates together, the statistic is 11.72 on 4 degrees of freedom (P = 0.020). The test finds evidence against the assumption for the model as a whole. For preterm birth alone, the statistic is 11.65 on 2 degrees of freedom (P = 0.003).
For ANC alone, the statistic is 0.001 on 2 degrees of freedom (P = 0.999). Because the counts were built without sampling noise, the ANC statistic sits almost at zero. In a real sample in which the assumption holds, the statistic varies from sample to sample. Its P value can then fall anywhere between 0 and 1.
The overall test says that something is wrong; the per-covariate tests say where. Add-on commands automate the test: brant in Stata, from the user-written SPost package, and the brant package in R. The code below uses neither: it computes the test directly from the three cut fits, so every step is visible.
Stata: the Brant test computed in Mata
* step 2, the Brant (1990) Wald test that the slopes are equal across the three cuts. For cuts j and l the
* covariance of the two slope vectors is A_j^-1 [X' diag(pi_hi (1 - pi_lo)) X] A_l^-1, with
* A_j = X' diag(pi_j (1 - pi_j)) X, lo = min(j, l) and hi = max(j, l); the contrasts compare cut 1 with
* cuts 2 and 3 (4 df for both covariates, 2 df for each covariate)
mata:
X = st_data(., ("inadequate_anc", "preterm")), J(st_nobs(), 1, 1)
P = st_data(., ("pi1", "pi2", "pi3"))
B1 = st_matrix("bcut1")
B2 = st_matrix("bcut2")
B3 = st_matrix("bcut3")
b = (B1[|1, 1 \ 1, 2|], B2[|1, 1 \ 1, 2|], B3[|1, 1 \ 1, 2|])'
V = J(6, 6, 0)
for (j = 1; j <= 3; j++) {
for (l = 1; l <= 3; l++) {
lo = min((j, l))
hi = max((j, l))
Ajl = invsym(cross(X, P[., j]:*(1 :- P[., j]), X)) * cross(X, P[., hi]:*(1 :- P[., lo]), X) * invsym(cross(X, P[., l]:*(1 :- P[., l]), X))
V[|2*j - 1, 2*l - 1 \ 2*j, 2*l|] = Ajl[|1, 1 \ 2, 2|]
}
}
D = (1, 0, -1, 0, 0, 0 \ 0, 1, 0, -1, 0, 0 \ 1, 0, 0, 0, -1, 0 \ 0, 1, 0, 0, 0, -1)
Da = D[(1 \ 3), .]
Dp = D[(2 \ 4), .]
st_numscalar("sc_brant_all", (D*b)' * invsym(D*V*D') * (D*b))
st_numscalar("sc_brant_anc", (Da*b)' * invsym(Da*V*Da') * (Da*b))
st_numscalar("sc_brant_pre", (Dp*b)' * invsym(Dp*V*Dp') * (Dp*b))
end
canon brant.chi2_all 4 scalar(sc_brant_all)
canon brant.df_all 0 4
canon brant.p_all 4 chi2tail(4, scalar(sc_brant_all))
canon brant.chi2_anc 4 scalar(sc_brant_anc)
canon brant.p_anc 4 chi2tail(2, scalar(sc_brant_anc))
canon brant.chi2_preterm 4 scalar(sc_brant_pre)
canon brant.p_preterm 4 chi2tail(2, scalar(sc_brant_pre))
canon brant.df_each 0 2. * step 2, the Brant (1990) Wald test that the slopes are equal across the three cuts. For cuts j and l the
. * covariance of the two slope vectors is A_j^-1 [X' diag(pi_hi (1 - pi_lo)) X] A_l^-1, with
. * A_j = X' diag(pi_j (1 - pi_j)) X, lo = min(j, l) and hi = max(j, l); the contrasts compare cut 1 with
. * cuts 2 and 3 (4 df for both covariates, 2 df for each covariate)
. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
: X = st_data(., ("inadequate_anc", "preterm")), J(st_nobs(), 1, 1)
: P = st_data(., ("pi1", "pi2", "pi3"))
: B1 = st_matrix("bcut1")
: B2 = st_matrix("bcut2")
: B3 = st_matrix("bcut3")
: b = (B1[|1, 1 \ 1, 2|], B2[|1, 1 \ 1, 2|], B3[|1, 1 \ 1, 2|])'
: V = J(6, 6, 0)
: for (j = 1; j <= 3; j++) {
> for (l = 1; l <= 3; l++) {
> lo = min((j, l))
> hi = max((j, l))
> Ajl = invsym(cross(X, P[., j]:*(1 :- P[., j]), X)) * cross(X, P[., hi]:*(1 :- P[., lo]), X) * invsym(cross(X, P[., l]:*(1 :- P[., l]), X))
> V[|2*j - 1, 2*l - 1 \ 2*j, 2*l|] = Ajl[|1, 1 \ 2, 2|]
> }
> }
: D = (1, 0, -1, 0, 0, 0 \ 0, 1, 0, -1, 0, 0 \ 1, 0, 0, 0, -1, 0 \ 0, 1, 0, 0, 0, -1)
: Da = D[(1 \ 3), .]
: Dp = D[(2 \ 4), .]
: st_numscalar("sc_brant_all", (D*b)' * invsym(D*V*D') * (D*b))
: st_numscalar("sc_brant_anc", (Da*b)' * invsym(Da*V*Da') * (Da*b))
: st_numscalar("sc_brant_pre", (Dp*b)' * invsym(Dp*V*Dp') * (Dp*b))
: end
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
. canon brant.chi2_all 4 scalar(sc_brant_all)
CANON c7.brant.chi2_all 11.7179
. canon brant.df_all 0 4
CANON c7.brant.df_all 4
. canon brant.p_all 4 chi2tail(4, scalar(sc_brant_all))
CANON c7.brant.p_all 0.0196
. canon brant.chi2_anc 4 scalar(sc_brant_anc)
CANON c7.brant.chi2_anc 0.0011
. canon brant.p_anc 4 chi2tail(2, scalar(sc_brant_anc))
CANON c7.brant.p_anc 0.9994
. canon brant.chi2_preterm 4 scalar(sc_brant_pre)
CANON c7.brant.chi2_preterm 11.6476
. canon brant.p_preterm 4 chi2tail(2, scalar(sc_brant_pre))
CANON c7.brant.p_preterm 0.0030
. canon brant.df_each 0 2
CANON c7.brant.df_each 2R: the cut fits and the Brant test
# the proportional-odds check, step 1: a separate binary logit at each cut (severity >= 1, >= 2, >= 3)
fits <- lapply(1:3, function(j) {
o$ge <- as.integer(o$severity >= j)
glm(ge ~ inadequate_anc + preterm, family = binomial, data = o)
})
for (j in 1:3) {
f <- fits[[j]]
canon(sprintf("cut%d.b_anc", j), coef(f)[["inadequate_anc"]])
canon(sprintf("cut%d.or_anc", j), exp(coef(f)[["inadequate_anc"]]))
canon(sprintf("cut%d.b_preterm", j), coef(f)[["preterm"]])
canon(sprintf("cut%d.or_preterm", j), exp(coef(f)[["preterm"]]))
}
# step 2, the Brant (1990) Wald test that the slopes are equal across the three cuts. For cuts j and l the
# covariance of the two slope vectors is A_j^-1 [X' diag(pi_hi (1 - pi_lo)) X] A_l^-1, with
# A_j = X' diag(pi_j (1 - pi_j)) X, lo = min(j, l) and hi = max(j, l); the contrasts compare cut 1 with
# cuts 2 and 3 (4 df for both covariates, 2 df for each covariate)
X <- cbind(inadequate_anc = o$inadequate_anc, preterm = o$preterm, cons = 1)
P <- sapply(fits, fitted)
b <- unlist(lapply(fits, function(f) coef(f)[c("inadequate_anc", "preterm")]))
V <- matrix(0, 6, 6)
for (j in 1:3) {
for (l in 1:3) {
lo <- min(j, l)
hi <- max(j, l)
Ajl <- solve(crossprod(X, X * (P[, j] * (1 - P[, j])))) %*%
crossprod(X, X * (P[, hi] * (1 - P[, lo]))) %*%
solve(crossprod(X, X * (P[, l] * (1 - P[, l]))))
V[(2 * j - 1):(2 * j), (2 * l - 1):(2 * l)] <- Ajl[1:2, 1:2]
}
}
D <- rbind(c(1, 0, -1, 0, 0, 0), c(0, 1, 0, -1, 0, 0), c(1, 0, 0, 0, -1, 0), c(0, 1, 0, 0, 0, -1))
wald <- function(Dm) {
db <- Dm %*% b
as.numeric(t(db) %*% solve(Dm %*% V %*% t(Dm)) %*% db)
}
chi_all <- wald(D)
chi_anc <- wald(D[c(1, 3), ])
chi_pre <- wald(D[c(2, 4), ])
canon("brant.chi2_all", chi_all)
canon("brant.df_all", 4, 0)
canon("brant.p_all", pchisq(chi_all, 4, lower.tail = FALSE))
canon("brant.chi2_anc", chi_anc)
canon("brant.p_anc", pchisq(chi_anc, 2, lower.tail = FALSE))
canon("brant.chi2_preterm", chi_pre)
canon("brant.p_preterm", pchisq(chi_pre, 2, lower.tail = FALSE))
canon("brant.df_each", 2, 0) > fits <- lapply(1:3, function(j) {
+ o$ge <- as.integer(o$severity >= j)
+ glm(ge ~ inadequate_anc + preterm, family = binomial, data = o)
+ })
> for (j in 1:3) {
+ f <- fits[[j]]
+ canon(sprintf("cut%d.b_anc", j), coef(f)[["inadequate_anc"]])
+ canon(sprintf("cut%d.or_anc", j), exp(coef(f)[["inadequate_anc"]]))
+ canon(sprintf("cut%d.b_preterm", j), coef(f)[["preterm"]])
+ canon(sprintf("cut%d.or_preterm", j), exp(coef(f)[["preterm"]]))
+ }
CANON c7.cut1.b_anc 0.6913
CANON c7.cut1.or_anc 1.9963
CANON c7.cut1.b_preterm 0.3354
CANON c7.cut1.or_preterm 1.3985
CANON c7.cut2.b_anc 0.6965
CANON c7.cut2.or_anc 2.0067
CANON c7.cut2.b_preterm 0.7170
CANON c7.cut2.or_preterm 2.0483
CANON c7.cut3.b_anc 0.6984
CANON c7.cut3.or_anc 2.0105
CANON c7.cut3.b_preterm 1.1974
CANON c7.cut3.or_preterm 3.3115
> X <- cbind(inadequate_anc = o$inadequate_anc, preterm = o$preterm,
+ cons = 1)
> P <- sapply(fits, fitted)
> b <- unlist(lapply(fits, function(f) coef(f)[c("inadequate_anc",
+ "preterm")]))
> V <- matrix(0, 6, 6)
> for (j in 1:3) {
+ for (l in 1:3) {
+ lo <- min(j, l)
+ hi <- max(j, l)
+ Ajl <- solve(crossprod(X, X * (P[, j] * (1 - P[, j])))) %*%
+ crossprod(X, X * (P[, hi] * (1 - P[, lo]))) %*% solve(crossprod(X,
+ X * (P[, l] * (1 - P[, l]))))
+ V[(2 * j - 1):(2 * j), (2 * l - 1):(2 * l)] <- Ajl[1:2,
+ 1:2]
+ }
+ }
> D <- rbind(c(1, 0, -1, 0, 0, 0), c(0, 1, 0, -1, 0,
+ 0), c(1, 0, 0, 0, -1, 0), c(0, 1, 0, 0, 0, -1))
> wald <- function(Dm) {
+ db <- Dm %*% b
+ as.numeric(t(db) %*% solve(Dm %*% V %*% t(Dm)) %*% db)
+ }
> chi_all <- wald(D)
> chi_anc <- wald(D[c(1, 3), ])
> chi_pre <- wald(D[c(2, 4), ])
> canon("brant.chi2_all", chi_all)
CANON c7.brant.chi2_all 11.7179
> canon("brant.df_all", 4, 0)
CANON c7.brant.df_all 4
> canon("brant.p_all", pchisq(chi_all, 4, lower.tail = FALSE))
CANON c7.brant.p_all 0.0196
> canon("brant.chi2_anc", chi_anc)
CANON c7.brant.chi2_anc 0.0011
> canon("brant.p_anc", pchisq(chi_anc, 2, lower.tail = FALSE))
CANON c7.brant.p_anc 0.9994
> canon("brant.chi2_preterm", chi_pre)
CANON c7.brant.chi2_preterm 11.6476
> canon("brant.p_preterm", pchisq(chi_pre, 2, lower.tail = FALSE))
CANON c7.brant.p_preterm 0.0030
> canon("brant.df_each", 2, 0)
CANON c7.brant.df_each 2 When proportional odds fails
A flagged covariate does not make the ordinal model useless. Here the check flags preterm birth only. The ANC odds ratio of 2.01 still summarises three nearly equal cut odds ratios. The preterm odds ratio of 1.84 does not: it sits between cut odds ratios that run from 1.40 to 3.31.
Three alternatives relax the assumption; each is named here and set out in the posts linked below. A partial proportional odds model lets the flagged covariate, here preterm birth, have its own odds ratio at each cut. A generalized ordered logit model frees every covariate in this way. Multinomial logistic regression drops the order and compares each grade with a reference grade.
In this cohort, one reasonable report would keep the ordinal odds ratio for ANC and give preterm birth its three cut-specific odds ratios. A partial proportional odds model would deliver both in one fit. Either way, the report should say that the check was run and what it found.
Separate posts on ordinal logistic regression and on multinomial versus ordinal regression set out these options. Whichever is chosen, it may help readers to see the cut-specific odds ratios beside any summary. A Brant P value also depends on sample size: a large study can flag small departures, and a small one can miss large ones.
Third model: matched case-control pairs
Outcome type: binary again, case or control. Sampling design: matched sets, each a case with the control drawn for it. Model: conditional logistic regression, which compares each case only with its own control [4].
Each of the 20 cases had a complication within 28 days of birth. Its control was a newborn still free of complications on the day the case had one. The pair is therefore matched on age in days, and the control may itself become a case later. Drawing controls from those still at risk is risk-set sampling, set out in a separate post on control sampling.
A case and its control form a matched set, here a pair, with its own identifier in the data. The 20 pairs give 40 rows, one per pair member. The exposure is inadequate ANC, recorded for both members of every pair.
Matching types are covered in a post on conditional logistic regression. So is overmatching: matching on a factor tied to the exposure but not a confounder, which can cost precision.
| Case's ANC | Control's ANC | Pairs | Pair type |
|---|---|---|---|
| Inadequate | Inadequate | 8 | Concordant |
| Inadequate | Adequate | 4 | Discordant |
| Adequate | Inadequate | 1 | Discordant |
| Adequate | Adequate | 7 | Concordant |
| All pairs | Any | 20 | 15 concordant, 5 discordant |
Concordant and discordant pairs
A concordant pair is one in which case and control share exposure status. There are 15: 8 with both exposed and 7 with neither. A discordant pair is one in which they differ: only the case was exposed in 4 pairs, and only the control in 1. That split drives everything that follows.
| ANC | Cases | Controls |
|---|---|---|
| Inadequate | 12 | 9 |
| Adequate | 8 | 11 |
| Total | 20 | 20 |
Hand example: pooling the pairs
Hand example. Outcome type: binary. Sampling design: matched pairs, treated by this model as 40 independent newborns.
Ordinary logistic regression ignores which control belongs to which case. It sees 20 cases and 20 controls and pools them into the table above, exposure by case status.
-
Exposed cases
\[ 8 + 4 = 12 \]
Pairs with both exposed, plus pairs with only the case exposed. The other 7 + 1 = 8 cases had adequate ANC.
-
Exposed controls
\[ 8 + 1 = 9 \]
Pairs with both exposed, plus pairs with only the control exposed. The other 7 + 4 = 11 controls had adequate ANC.
-
Exposure odds among cases
\[ \frac{12}{8} = 1.50 \]
Inadequate against adequate ANC among the 20 cases.
-
Exposure odds among controls
\[ \frac{9}{11} = 0.818 \]
The same odds among the 20 controls.
-
Unmatched odds ratio
\[ \mathrm{OR} = \frac{1.50}{0.818} = 1.83 \]
The ratio of the two exposure odds, as in any unmatched case-control table.
Result: Pooled over the pairs, the odds of inadequate ANC are 1.83 times as high among cases as among controls.
Conditional logistic regression: comparing within pairs
Conditional logistic regression asks a narrower question. Given that exactly one member of a pair became a case, how does exposure change which member it was? Only pairs whose members differ in exposure can answer that.
In a concordant pair, both members share exposure, so exposure cannot explain why one of them is the case. Such a pair carries no information about the ANC coefficient. In a model with more covariates, concordant pairs can still inform the coefficients of covariates that differ within them.
The 5 discordant pairs carry all the information about exposure. In 4 of them the exposed member is the case; in 1, the unexposed member is. Take 1:1 pairs with one binary exposure and no other covariates. Then the conditional (matched) odds ratio, $\mathrm{OR}_c$, is the ratio of these two counts:
$$\mathrm{OR}_c = \frac{4}{1} = 4.00$$The numerator counts pairs in which only the case was exposed, and the denominator pairs in which only the control was. The next hand example shows where this ratio comes from [5].
Hand example: why the pair intercept cancels
Give pair $k$ its own intercept $\omega_k$, standing for everything the pair members share, and a common exposure slope $\beta$. Write $x_1$ for the exposure of the member who became the case and $x_2$ for the other member's. The likelihood is the probability of the observed data, read as a function of the coefficients. The conditional likelihood asks how probable it is that member 1, not member 2, is the case, given that the pair holds exactly one case [5].
-
Probability for one pair
\[ P_k = \frac{e^{\omega_k + \beta x_1}}{e^{\omega_k + \beta x_1} + e^{\omega_k + \beta x_2}} \]
Under the logistic model, $e^{\omega_k + \beta x}$ is a member's odds of becoming a case. Given one case in the pair, the case's probability is its share of the two odds.
-
The pair intercept cancels
\[ P_k = \frac{e^{\beta x_1}}{e^{\beta x_1} + e^{\beta x_2}} \]
Dividing top and bottom by $e^{\omega_k}$ removes the intercept, so nothing the pair members share needs to be estimated.
-
A concordant pair
\[ x_1 = x_2: \quad P_k = \frac{e^{\beta x_1}}{2\, e^{\beta x_1}} = 0.5 \]
Whatever $\beta$ is, the probability is 0.5, so the 15 concordant pairs carry no information about $\beta$.
-
A discordant pair
\[ q = \frac{e^{\beta}}{1 + e^{\beta}} \]
Write $q$ for the probability that the exposed member is the case. The unexposed member is the case with probability $1 - q$.
-
The conditional likelihood
\[ L(\beta) \propto q^{4}\,(1 - q)^{1} \]
In 4 discordant pairs the exposed member is the case, and in 1 the unexposed member is. The concordant pairs only multiply $L(\beta)$ by a constant.
-
The maximum
\[ q = \frac{4}{5} = 0.80 \]
The product $q^{4}(1 - q)$ is largest when $q$ equals the share of discordant pairs in which the exposed member is the case. Then $1 - q = 0.20$.
-
Back to the odds ratio
\[ e^{\beta} = \frac{q}{1 - q} = \frac{0.80}{0.20} = 4.00 \]
Solving the discordant-pair step for $e^{\beta}$ gives the odds $q/(1 - q)$, which is the matched odds ratio $\mathrm{OR}_c$.
Result: The conditional likelihood gives the matched odds ratio of 4.00, the ratio 4/1 of the two kinds of discordant pair; the concordant pairs drop out.
The same 20 pairs in Stata and R
The simulation script types in the four pair counts, expands them to 40 rows and fits both models. Ordinary logistic regression returns 1.83, with a 95% confidence interval of 0.52 to 6.43. Conditional logistic regression returns 4.00, with a 95% confidence interval of 0.45 to 35.79. The Stata code also refits the ordinary model with its standard error clustered on the pair, a variant taken up under the common misreadings near the end.
The matched estimate is imprecise. Its coefficient is 1.386, the log of 4.00. On that log scale its standard error is $\sqrt{1/4 + 1/1} = 1.118$. That formula holds for 1:1 pairs with one binary exposure and no other covariates.
The squared standard error, the variance, is the sum $1/4 + 1/1$. Most of it comes from the single pair in which only the control was exposed. With 5 informative pairs, the interval runs from below 1 to far above 4.00. The data therefore cannot rule out no association.
This Wald interval, built from the estimate and its standard error, is only approximate with so few informative pairs. An exact interval, based directly on the split of 4 against 1, would be wider still.
Stata labels the clogit output conditional (fixed-effects) logistic regression. The fixed effects are the pair intercepts $\omega_k$, which the conditional likelihood removes instead of estimating. The output therefore has no constant term.
R's clogit() prints a call to coxph(), R's Cox regression for time-to-event data, because it fits the conditional likelihood through that routine. The post on control sampling follows that link and what it means for a matched odds ratio.
Stata: entering the 20 pairs
clear
input byte anc_case byte anc_ctrl int n_pairs
1 1 8
0 0 7
1 0 4
0 1 1
end
quietly summarize n_pairs
canon pairs.n 0 r(sum)
quietly summarize n_pairs if anc_case == anc_ctrl
canon pairs.concordant 0 r(sum)
quietly summarize n_pairs if anc_case != anc_ctrl
canon pairs.discordant 0 r(sum)
quietly summarize n_pairs if anc_case == 1 & anc_ctrl == 0
canon pairs.case_only_inadequate 0 r(sum)
quietly summarize n_pairs if anc_case == 0 & anc_ctrl == 1
canon pairs.control_only_inadequate 0 r(sum)
expand n_pairs
drop n_pairs
generate int set_id = _n
reshape long anc_, i(set_id) j(role) string
rename anc_ inadequate_anc
generate byte case = role == "case"
bysort set_id: assert _N == 2
bysort set_id (case): assert case[1] == 0 & case[2] == 1
quietly count
canon pairs.rows 0 r(N)
* the pooled 2 x 2 table that ordinary logistic regression sees
tabulate inadequate_anc case. clear
. input byte anc_case byte anc_ctrl int n_pairs
anc_case anc_ctrl n_pairs
1. 1 1 8
2. 0 0 7
3. 1 0 4
4. 0 1 1
5. end
. quietly summarize n_pairs
. canon pairs.n 0 r(sum)
CANON c7.pairs.n 20
. quietly summarize n_pairs if anc_case == anc_ctrl
. canon pairs.concordant 0 r(sum)
CANON c7.pairs.concordant 15
. quietly summarize n_pairs if anc_case != anc_ctrl
. canon pairs.discordant 0 r(sum)
CANON c7.pairs.discordant 5
. quietly summarize n_pairs if anc_case == 1 & anc_ctrl == 0
. canon pairs.case_only_inadequate 0 r(sum)
CANON c7.pairs.case_only_inadequate 4
. quietly summarize n_pairs if anc_case == 0 & anc_ctrl == 1
. canon pairs.control_only_inadequate 0 r(sum)
CANON c7.pairs.control_only_inadequate 1
. expand n_pairs
(16 observations created)
. drop n_pairs
. generate int set_id = _n
. reshape long anc_, i(set_id) j(role) string
(j = case ctrl)
Data Wide -> Long
-----------------------------------------------------------------------------
Number of observations 20 -> 40
Number of variables 3 -> 3
j variable (2 values) -> role
xij variables:
anc_case anc_ctrl -> anc_
-----------------------------------------------------------------------------
. rename anc_ inadequate_anc
. generate byte case = role == "case"
. bysort set_id: assert _N == 2
. bysort set_id (case): assert case[1] == 0 & case[2] == 1
. quietly count
. canon pairs.rows 0 r(N)
CANON c7.pairs.rows 40
. * the pooled 2 x 2 table that ordinary logistic regression sees
. tabulate inadequate_anc case
inadequate | case
_anc | 0 1 | Total
-----------+----------------------+----------
0 | 11 8 | 19
1 | 9 12 | 21
-----------+----------------------+----------
Total | 20 20 | 40Stata: ordinary logistic regression on the pairs
* ordinary logistic regression ignores the matched sets
logistic case i.inadequate_anc, nolog. * ordinary logistic regression ignores the matched sets
. logistic case i.inadequate_anc, nolog
Logistic regression Number of obs = 40
LR chi2(1) = 0.91
Prob > chi2 = 0.3413
Log likelihood = -27.27303 Pseudo R2 = 0.0163
----------------------------------------------------------------------------------
case | Odds ratio Std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc | 1.833333 1.174412 0.95 0.344 .522362 6.434448
_cons | .7272727 .3379347 -0.69 0.493 .292535 1.808076
----------------------------------------------------------------------------------
Note: _cons estimates baseline odds.Stata: the pooled fit with a cluster-robust standard error
* clustering the standard error on the matched set changes the interval, not the odds ratio
logistic case i.inadequate_anc, vce(cluster set_id) nolog. * clustering the standard error on the matched set changes the interval, not the odds ratio
. logistic case i.inadequate_anc, vce(cluster set_id) nolog
Logistic regression Number of obs = 40
Wald chi2(1) = 1.82
Prob > chi2 = 0.1770
Log pseudolikelihood = -27.27303 Pseudo R2 = 0.0163
(Std. err. adjusted for 20 clusters in set_id)
----------------------------------------------------------------------------------
| Robust
case | Odds ratio std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc | 1.833333 .823118 1.35 0.177 .7604568 4.419858
_cons | .7272727 .1778604 -1.30 0.193 .4503282 1.174534
----------------------------------------------------------------------------------
Note: _cons estimates baseline odds.R: the pairs and the ordinary logistic fit
pairs <- data.frame(anc_case = c(1, 0, 1, 0), anc_ctrl = c(1, 0, 0, 1), n_pairs = c(8, 7, 4, 1))
canon("pairs.n", sum(pairs$n_pairs), 0)
canon("pairs.concordant", sum(pairs$n_pairs[pairs$anc_case == pairs$anc_ctrl]), 0)
canon("pairs.discordant", sum(pairs$n_pairs[pairs$anc_case != pairs$anc_ctrl]), 0)
canon("pairs.case_only_inadequate", sum(pairs$n_pairs[pairs$anc_case == 1 & pairs$anc_ctrl == 0]), 0)
canon("pairs.control_only_inadequate", sum(pairs$n_pairs[pairs$anc_case == 0 & pairs$anc_ctrl == 1]), 0)
pp <- pairs[rep(seq_len(nrow(pairs)), pairs$n_pairs), ]
pp$set_id <- seq_len(nrow(pp))
m <- rbind(data.frame(set_id = pp$set_id, case = 1L, inadequate_anc = pp$anc_case),
data.frame(set_id = pp$set_id, case = 0L, inadequate_anc = pp$anc_ctrl))
m <- m[order(m$set_id, m$case), ]
stopifnot(all(table(m$set_id) == 2), all(tapply(m$case, m$set_id, sum) == 1))
canon("pairs.rows", nrow(m), 0)
# the pooled 2 x 2 table that ordinary logistic regression sees
print(table(inadequate_anc = m$inadequate_anc, case = m$case))
a <- sum(m$case == 1 & m$inadequate_anc == 1)
bb <- sum(m$case == 1 & m$inadequate_anc == 0)
cc <- sum(m$case == 0 & m$inadequate_anc == 1)
dd <- sum(m$case == 0 & m$inadequate_anc == 0)
canon("pooled.cases_inadequate", a, 0)
canon("pooled.cases_adequate", bb, 0)
canon("pooled.controls_inadequate", cc, 0)
canon("pooled.controls_adequate", dd, 0)
canon("pooled.odds_cases", a / bb)
canon("pooled.odds_controls", cc / dd)
# ordinary logistic regression ignores the matched sets
fp <- glm(case ~ inadequate_anc, family = binomial, data = m)
print(summary(fp)$coefficients)
bp <- coef(fp)[["inadequate_anc"]]
sp <- sqrt(vcov(fp)["inadequate_anc", "inadequate_anc"])
canon("pooled.or", exp(bp))
canon("pooled.se", sp)
canon("pooled.or_lb", exp(bp - z * sp))
canon("pooled.or_ub", exp(bp + z * sp)) > pairs <- data.frame(anc_case = c(1, 0, 1, 0), anc_ctrl = c(1,
+ 0, 0, 1), n_pairs = c(8, 7, 4, 1))
> canon("pairs.n", sum(pairs$n_pairs), 0)
CANON c7.pairs.n 20
> canon("pairs.concordant", sum(pairs$n_pairs[pairs$anc_case ==
+ pairs$anc_ctrl]), 0)
CANON c7.pairs.concordant 15
> canon("pairs.discordant", sum(pairs$n_pairs[pairs$anc_case !=
+ pairs$anc_ctrl]), 0)
CANON c7.pairs.discordant 5
> canon("pairs.case_only_inadequate", sum(pairs$n_pairs[pairs$anc_case ==
+ 1 & pairs$anc_ctrl == 0]), 0)
CANON c7.pairs.case_only_inadequate 4
> canon("pairs.control_only_inadequate", sum(pairs$n_pairs[pairs$anc_case ==
+ 0 & pairs$anc_ctrl == 1]), 0)
CANON c7.pairs.control_only_inadequate 1
> pp <- pairs[rep(seq_len(nrow(pairs)), pairs$n_pairs),
+ ]
> pp$set_id <- seq_len(nrow(pp))
> m <- rbind(data.frame(set_id = pp$set_id, case = 1,
+ inadequate_anc = pp$anc_case), data.frame(set_id = pp$set_id,
+ case = 0, inadequate_anc = pp$anc_ctrl))
> m <- m[order(m$set_id, m$case), ]
> stopifnot(all(table(m$set_id) == 2), all(tapply(m$case,
+ m$set_id, sum) == 1))
> canon("pairs.rows", nrow(m), 0)
CANON c7.pairs.rows 40
> print(table(inadequate_anc = m$inadequate_anc, case = m$case))
case
inadequate_anc 0 1
0 11 8
1 9 12
> a <- sum(m$case == 1 & m$inadequate_anc == 1)
> bb <- sum(m$case == 1 & m$inadequate_anc == 0)
> cc <- sum(m$case == 0 & m$inadequate_anc == 1)
> dd <- sum(m$case == 0 & m$inadequate_anc == 0)
> canon("pooled.cases_inadequate", a, 0)
CANON c7.pooled.cases_inadequate 12
> canon("pooled.cases_adequate", bb, 0)
CANON c7.pooled.cases_adequate 8
> canon("pooled.controls_inadequate", cc, 0)
CANON c7.pooled.controls_inadequate 9
> canon("pooled.controls_adequate", dd, 0)
CANON c7.pooled.controls_adequate 11
> canon("pooled.odds_cases", a/bb)
CANON c7.pooled.odds_cases 1.5000
> canon("pooled.odds_controls", cc/dd)
CANON c7.pooled.odds_controls 0.8182
> fp <- glm(case ~ inadequate_anc, family = binomial,
+ data = m)
> print(summary(fp)$coefficients)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.3184537 0.4646602 -0.6853476 0.4931246
inadequate_anc 0.6061358 0.6405884 0.9462172 0.3440378
> bp <- coef(fp)[["inadequate_anc"]]
> sp <- sqrt(vcov(fp)["inadequate_anc", "inadequate_anc"])
> canon("pooled.or", exp(bp))
CANON c7.pooled.or 1.8333
> canon("pooled.se", sp)
CANON c7.pooled.se 0.6406
> canon("pooled.or_lb", exp(bp - z * sp))
CANON c7.pooled.or_lb 0.5224
> canon("pooled.or_ub", exp(bp + z * sp))
CANON c7.pooled.or_ub 6.4344 Stata: conditional logistic regression on the pairs
* conditional logistic regression compares the case with its own control inside each matched set
clogit case i.inadequate_anc, group(set_id) or nolog. * conditional logistic regression compares the case with its own control inside each matched set
. clogit case i.inadequate_anc, group(set_id) or nolog
Conditional (fixed-effects) logistic regression Number of obs = 40
LR chi2(1) = 1.93
Prob > chi2 = 0.1650
Log likelihood = -12.89922 Pseudo R2 = 0.0695
----------------------------------------------------------------------------------
case | Odds ratio Std. err. z P>|z| [95% conf. interval]
-----------------+----------------------------------------------------------------
1.inadequate_anc | 4 4.472136 1.24 0.215 .4470826 35.78757
----------------------------------------------------------------------------------R: clogit() on the pairs
# conditional logistic regression compares the case with its own control inside each matched set
fc <- clogit(case ~ inadequate_anc + strata(set_id), data = m)
print(summary(fc))> fc <- clogit(case ~ inadequate_anc + strata(set_id),
+ data = m)
> print(summary(fc))
Call:
coxph(formula = Surv(rep(1, 40L), case) ~ inadequate_anc + strata(set_id),
data = m, method = "exact")
n= 40, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
inadequate_anc 1.386 4.000 1.118 1.24 0.215
exp(coef) exp(-coef) lower .95 upper .95
inadequate_anc 4 0.25 0.4471 35.79
Concordance= 0.575 (se = 0.075 )
Likelihood ratio test= 1.93 on 1 df, p=0.2
Wald test = 1.54 on 1 df, p=0.2
Score (logrank) test = 1.8 on 1 df, p=0.2Why 1.83 and 4.00 differ
The two numbers answer different questions. The unmatched odds ratio compares the exposure odds of all cases with those of all controls. The matched odds ratio compares exposure within pairs, between newborns who share whatever the pair members share.
Pooling blurs the within-pair contrast. Every concordant pair adds the same exposure to both columns of the pooled table. Pairs with both exposed add 8 exposed members to each column; pairs with neither add 7 unexposed members to each. Here 15 of the 20 pairs are concordant.
With inadequate ANC in 12 of 20 cases and 9 of 20 controls, exposure unrelated within a pair would make about 10 pairs concordant (0.60 x 0.45 + 0.40 x 0.55 = 0.49, so 0.49 x 20 = 9.8), against 15 observed. That is compatible with the matching factor being related to exposure, though 20 pairs cannot show it.
When the matching factor, here age in days, is related to exposure, matching pulls the controls' exposure toward the cases'. An analysis that ignores the sets then tends to bias the odds ratio toward 1 [5]. The shared additions above are how that bias shows up in the pooled table.
When the matching factor is unrelated to exposure, the matched and pooled analyses estimate the same quantity. They then differ mainly in precision. With only 5 informative pairs, chance can also open a gap between them.
Non-collapsibility, defined below for the cohort models, is a gap between an odds ratio within subgroups and one averaged over them, which appears when the subgroups predict the outcome even without confounding. Here every pair holds one case and one control of the same age in days, so neither the pair nor age can predict who is a case. The gap between 1.83 and 4.00 therefore does not come from non-collapsibility. It reflects the design, a within-pair against a pooled comparison in which matching may pull the controls' exposure toward the cases', together with chance in so few informative pairs.
So which number belongs in the report? The study was designed for the within-pair comparison, so the matched odds ratio is usually the one to report. It should come with its interval, 0.45 to 35.79, which is compatible both with no association and with a strong one.
Matching deals only with the matching factors. Other confounders still enter the conditional model as covariates, and their coefficients are estimated within the matched sets. A matched odds ratio without them is adjusted for age in days and nothing else.
One link, three models: a decision table
The table collects the three models. Read each row from the left: the outcome type and the sampling design come first, and the model and the commands follow. The last column says what the exponentiated coefficient means, which is where the three models differ most.
The last row uses a hazard ratio, which compares instantaneous event rates among people still free of the event. Under proportional hazards, that ratio stays the same throughout follow-up.
| Outcome type and sampling design | Model | Stata | R | Exponentiated coefficient |
|---|---|---|---|---|
| Binary outcome; independent people, for example a cohort or an unmatched case-control study | Logistic regression | logistic y i.x | glm(y ~ x, family = binomial) | Odds ratio for exposed versus unexposed, holding the other covariates fixed; in a case-control study, what it estimates depends on how the controls were sampled (the odds ratio only for controls from people still free of the outcome at the end) |
| Ordered categories; independent people | Ordinal logistic regression (proportional odds) | ologit y i.x | MASS::polr(factor(y, ordered = TRUE) ~ x) | One odds ratio for a higher category, assumed equal at every cut point; check it with a Brant test |
| Binary (case or control); matched sets, including risk-set sampled controls | Conditional logistic regression | clogit case i.x, group(set_id) or | survival::clogit(case ~ x + strata(set_id)) | Odds ratio within matched sets; with risk-set sampling it estimates the hazard ratio, under proportional hazards; matching factors get no coefficient |
What the three odds ratios have in common
In all three models, a coefficient is a difference in log odds and its exponent is an odds ratio. What that odds ratio is conditional on differs. Logistic and ordinal logistic regression condition on the other covariates; conditional logistic regression also conditions on the matched set. Each conditioning set can define its own estimand, the precise quantity a study sets out to estimate.
In the birth cohort, the ANC odds ratios of 2.00 and 2.01 are conditional on preterm birth. Preterm birth is a prognostic factor, a characteristic that predicts the outcome. A model without it would give a marginal odds ratio, averaged over preterm status. Odds ratios are non-collapsible, so the two need not agree even if preterm birth were unrelated to ANC [6].
Conditional and marginal odds ratios, and conditional and marginal hazard ratios, can differ even when there is no confounding. In a randomised trial, adjusting for a strong prognostic factor typically moves the odds ratio away from 1. That is a change of estimand, not the removal of bias. The randomised-trial case [7] is worked through in the part of this series on conditional and marginal odds ratios.
A report should therefore name what each odds ratio is conditional on, beside the model. Inadequate ANC may lead to preterm birth, making preterm birth a mediator, a step on the path to complications. If so, the odds ratios of 2.00 and 2.01 omit that path, unlike ANC's total effect, which includes it. Which covariates to condition on is a causal question, best settled with a causal diagram before the analysis.
Common misreadings and their fixes
-
"clogit gives the larger odds ratio, so it is more accurate."
The size of an estimate says nothing about its accuracy. Here the matched 4.00 is further from 1 than the unmatched 1.83, but its interval, 0.45 to 35.79, is far wider than 0.52 to 6.43. In other data the matched odds ratio can be the smaller one.
Fix: The two numbers answer different questions: the matched odds ratio compares exposure within pairs, the unmatched one pools the pairs. Choose the model from the design and the estimand, and report the interval.
-
"ologit fits any ordered outcome."
In the simulated cohort, the ordinal model gave preterm birth one odds ratio, 1.84, while its cut odds ratios ran from 1.40 to 3.31. The Brant test flagged preterm birth (P = 0.003) and not ANC (P = 0.999).
Fix: ologit assumes one odds ratio at every cut point. Check that assumption and say what you did when it fails.
-
"clogit throws away the concordant pairs."
Every pair enters the conditional likelihood. With one binary exposure and no other covariates, a concordant pair contributes a probability of 0.5 whatever the coefficient, so it carries no information about it. With more covariates, it can still inform those that differ within the pair.
Fix: Keep every matched set in the data, and report how many sets are discordant on the exposure, since they set the precision of its odds ratio.
-
"A cluster-robust standard error repairs an analysis that ignored the matching."
A cluster-robust standard error allows observations within a cluster, here a pair, to be correlated. Clustering the pooled fit on the pair, as in the Stata output shown with the pairs, narrows its 95% interval to 0.76 to 4.42, against 0.52 to 6.43. It narrows here because each pair holds exactly one case by design. The odds ratio stays at 1.83, still the pooled comparison, so neither interval answers the within-pair question.
Fix: Analyse matched sets with conditional logistic regression. A different standard error changes the interval, never the estimate or the question it answers.
-
"Ordinary logistic regression with one indicator per pair gives the matched odds ratio."
Each pair indicator is estimated from only two newborns. Ordinary maximum likelihood, the usual fitting method, then pushes the odds ratio further from 1. For 1:1 pairs with one binary exposure, that model's odds ratio is the square of the matched one [5].
Fix: Condition on the matched sets with clogit: the conditional likelihood removes the pair intercepts instead of estimating them.
-
"An odds ratio of 4.00 means a fourfold risk of complications."
The matched odds ratio is a ratio of odds within pairs, not a ratio of risks. Here each control comes from newborns still at risk at its case's event time. Under proportional hazards, that design makes it estimate the hazard ratio, not a risk ratio at 28 days. The post on control sampling works this through.
Fix: Name the measure that the design estimates, and report 4.00 with its 95% confidence interval, 0.45 to 35.79.
What to do in your own analysis
- Write the outcome type and the sampling design into the analysis plan before choosing a command.
- For an ordered outcome, consider ordinal logistic regression rather than collapsing the grades, and report the cut-specific odds ratios beside the common one.
- Check proportional odds with a Brant test, overall and for each covariate, and state what was done for any covariate it flags.
- For matched or risk-set sampled controls, build the set identifier at the sampling step, not afterwards, and analyse the sets with conditional logistic regression.
- Report the number of discordant sets and the 95% confidence interval beside a matched odds ratio.
- If a matched and an unmatched analysis disagree, consider reporting the matched analysis as the primary estimate. The unmatched odds ratio describes the pooled table and ignores the matching. It may be biased toward 1 when a matching factor is related to exposure.
Glossary
- logit link
- The link ln{p/(1 - p)}, the log odds, which ties a linear predictor to a probability; its coefficients are log odds ratios.
- cut point (จุดตัดระหว่างหมวดของผลลัพธ์แบบลำดับ)
- A boundary between two neighbouring categories of an ordered outcome; each cut point has its own intercept in ordinal logistic regression.
- cumulative logit
- The log odds of being at or above a given category of an ordered outcome, one for each cut point.
- ordinal logistic regression
- A model for an ordered outcome that fits a cumulative logit at every cut point with one shared set of slopes.
- proportional odds assumption (ข้อสมมติ proportional odds)
- The assumption that a covariate's odds ratio is the same at every cut point of an ordered outcome.
- Brant test
- A Wald test of the proportional odds assumption that compares the slopes of separate cut fits, overall and covariate by covariate.
- matched set (ชุดที่จับคู่)
- One case together with the controls drawn or matched for it; a pair in 1:1 matching.
- concordant pair (คู่ที่สอดคล้องกัน)
- A matched pair in which case and control share the same exposure; with one binary exposure it carries no information about the exposure in conditional logistic regression.
- discordant pair (คู่ที่ไม่สอดคล้องกัน)
- A matched pair in which case and control differ in exposure; with 1:1 matching, one binary exposure and no other covariates, the matched odds ratio is the ratio of the two kinds.
- conditional likelihood
- The probability of which member of a matched set is the case, given that the set holds exactly one case; the set intercepts cancel from it.
- conditional logistic regression
- A logistic model with its own intercept for each matched set, fitted through the conditional likelihood so that each case is compared only with its own controls.
- risk-set sampling (การสุ่มกลุ่มควบคุมจากกลุ่มเสี่ยง)
- Drawing each case's controls from the people still at risk at the case's event time.
- non-collapsibility (ความไม่ยุบรวมได้)
- The property of the odds ratio and the hazard ratio whereby a marginal value can differ from the common conditional value without any confounding.
- mediator (ตัวแปรสื่อกลาง)
- A variable on the causal path from exposure to outcome through which part of the effect may pass.
- cluster-robust standard error
- A standard error that allows observations within a cluster, such as a matched pair, to be correlated; it changes the interval, not the estimate.
References
- StataCorp. Stata Base Reference Manual: [R] logit. Stata documentation (PDF) at stata.com. Accessed 11 October 2026. https://www.stata.com/manuals/rlogit.pdf
- StataCorp. Stata Base Reference Manual: [R] ologit. Stata documentation (PDF) at stata.com. Accessed 11 October 2026. https://www.stata.com/manuals/rologit.pdf
- Brant R. Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics. 1990;46(4):1171 ff. doi:10.2307/2532457 https://doi.org/10.2307/2532457
- StataCorp. Stata Base Reference Manual: [R] clogit. Stata documentation (PDF) at stata.com. Accessed 11 October 2026. https://www.stata.com/manuals/rclogit.pdf
- Breslow NE, Day NE. Statistical Methods in Cancer Research, Volume I: The Analysis of Case-Control Studies. IARC Scientific Publications No. 32. Lyon: International Agency for Research on Cancer; 1980. https://publications.iarc.who.int/Book-And-Report-Series/Iarc-Scientific-Publications/Statistical-Methods-In-Cancer-Research-Volume-I-The-Analysis-Of-Case-Control-Studies-1980
- Greenland S, Pearl J, Robins JM. Confounding and collapsibility in causal inference. Stat Sci. 1999;14(1):29-46. doi:10.1214/ss/1009211805 https://doi.org/10.1214/ss/1009211805
- Gail MH, Wieand S, Piantadosi S. Biased estimates of treatment effect in randomized experiments with nonlinear regressions and omitted covariates. Biometrika. 1984;71(3):431-444. doi:10.1093/biomet/71.3.431 https://doi.org/10.1093/biomet/71.3.431
Key takeaways
- Logistic, ordinal logistic and conditional logistic regression share the logit link, so each exponentiated coefficient is an odds ratio.
- The outcome type and the sampling design, named before any command, choose among the three models.
- Ordinal logistic regression assumes one odds ratio at every cut point; check it overall and per covariate with a Brant test.
- In 1:1 matched pairs with one binary exposure, only discordant pairs inform the exposure odds ratio: 4 against 1 gives 4.00, while pooling the pairs gives 1.83.
- The matched estimate answers the design's question, but with 5 discordant pairs its interval of 0.45 to 35.79 shows it is imprecise.
Related in the wiki: [[conditional-vs-marginal-or]] [[case-control-sampling-schemes-stata-and-r]] [[ordinal-logistic-regression]] [[conditional-logistic-regression]] [[multinomial-vs-ordinal-regression]]