Control Sampling in Stata and R: One Cohort, Three Control Series, Three Estimates
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
A case-control study measures exposure in the cases and in a sample of controls, and what its cross-product (the ratio printed as an odds ratio) estimates depends on where the controls came from. Three control series are drawn from one simulated cohort of 10,000 people with 1,560 cases. Controls from people still free of disease at the end of follow-up estimate the cohort odds ratio, 3.14. Controls from everyone at the start, including people who later become cases, estimate the risk ratio, 2.50. Controls from people still at risk at each case's event time, analysed as matched sets, estimate the hazard ratio, a ratio of instantaneous event rates: 2.78. One saved draw of each gives an odds ratio of 2.99, a risk ratio (cross-product) of 2.53 and a hazard ratio of 2.61, each interval covering its target. This article concludes that a case-control estimate can be read only once its control series and analysis are named.
A biobank budget and three ways to choose controls
A research team holds a cohort of 10,000 adults followed for 5 years with no loss to follow-up. Over that time, 1,560 of them developed the disease under study. Blood stored at enrolment could show who carried a suspected exposure, but every assay costs money. The budget covers 4,680 assays: all 1,560 cases plus two controls per case, 3,120 controls in all.
A case-control study works this way. It measures the exposure in the cases and in a sample of other cohort members, the controls, instead of in everyone. The rule that picks the controls, together with the controls it yields, makes up the control series. Three colleagues propose three rules:
- draw controls from the people still free of the disease when follow-up ends, called exclusive or cumulative sampling;
- draw controls from the whole cohort as it stood at enrolment, including people who later become cases, called inclusive or case-base sampling;
- each time a case occurs, draw controls from the people still at risk at that moment, called concurrent or risk-set sampling.
All three will assay the same 1,560 cases. They will not get the same number, and each will be right about a different quantity. Published case-control studies often leave unclear how their controls were sampled, and so what their odds ratio estimates [1]. This article runs the three designs on one simulated cohort in Stata and R, and checks each estimate against its target.
One simulated cohort
The cohort here is simulated data. Of 2,000 exposed people, 600 became cases over 8,500 person-years. Of 8,000 unexposed people, 960 became cases over 37,600 person-years. Person-years add up the time each person spent under follow-up and free of the disease, so a person who never becomes a case adds 5 years.
Every control series here is drawn from one defined cohort, so every design is nested in that cohort. In standard usage, the nested case-control design means the third rule. Its controls come from each case's risk set, the people still at risk when that case occurs. The cohort totals match the worked example of the article on the nested case-control design, which sets out the design reasoning in full.
This article concentrates on the commands and on reading their output. Each data file holds one row per person or per sampled person, with an identifier, id, the exposure, exposed (1 or 0), and the case status. For each design, Stata and R analyse the same saved draw of controls, called the frozen sample here. The files are not published; the code panes show the commands and their real output.
What the controls stand for
Every analysis below starts from a 2 by 2 table of exposure among the cases and the controls. Write $a$ and $b$ for the exposed and unexposed cases, and $c$ and $d$ for the exposed and unexposed controls. The table's cross-product, also called the cross-product ratio, is
$$\text{cross-product} = \frac{a \times d}{b \times c} = \frac{a / b}{c / d}$$The second form shows what it is: the exposure odds of the cases divided by the exposure odds of the controls. The cases' side is 600 to 960 in every design here, so the answer turns on what the controls stand for. Stata's cc and logistic, with one binary exposure, both return this number and print it as an odds ratio. In R, the exponentiated coefficient of the same logistic regression gives the same number.
What the cross-product estimates depends on what the controls mirror [2, 3, 4]:
- the people still free of the disease at the end: the odds ratio;
- everyone at the start: the risk ratio;
- the person-time at risk: the rate ratio, and, analysed as matched sets (each case with the controls drawn for it), the hazard ratio, a ratio of instantaneous event rates.
The first two readings assume a closed cohort, in which everyone enters at the start and is followed for the same fixed period with no loss to follow-up. This cohort is closed, so those two targets are its 5-year odds ratio and 5-year risk ratio. Risk-set sampling does not need that condition, because each risk set holds only the people still under follow-up when its case occurs. Its target, the rate or hazard ratio, stays clear when people enter at different times or are lost to follow-up.
None of these targets needs the disease to be rare. Rarity matters only when an odds ratio is read as a risk ratio.
Hand example: one cohort, three denominators
Hand example, from the cohort totals above (simulated data). Subscript 1 marks the exposed group and 0 the unexposed group. Steps 1 to 3 use the whole cohort. Steps 4 to 6 use rounded, illustrative tables of 1,000 controls each, the same tables as in the article on the nested case-control design.
-
Risks and the risk ratio
\[ \text{risk}_1 = \frac{600}{2000} = 0.300, \;\; \text{risk}_0 = \frac{960}{8000} = 0.120 \]
The risk ratio is 0.300 / 0.120 = 2.50. Its denominators are the people at the start.
-
Odds and the odds ratio
\[ \text{odds}_1 = \frac{600}{1400} = 0.429, \;\; \text{odds}_0 = \frac{960}{7040} = 0.136 \]
The odds divide the cases by the 1,400 and 7,040 people still free of disease at the end. The odds ratio is (600 × 7,040) / (1,400 × 960) = 3.14.
-
Rates and the rate ratio
\[ \text{rate}_1 = \frac{600}{8500} = 0.0706, \;\; \text{rate}_0 = \frac{960}{37600} = 0.0255 \]
Rates are cases per person-year. From the unrounded rates, the rate ratio is 2.76.
-
Exclusive controls
\[ \frac{600 \times 834}{960 \times 166} = 3.14 \]
Draw 1,000 controls from the 8,440 people still free of disease, split like the 1,400 and 7,040: 166 exposed, 834 unexposed. The cross-product returns the odds ratio.
-
Inclusive controls
\[ \frac{600 \times 800}{960 \times 200} = 2.50 \]
Draw 1,000 controls from all 10,000 people at the start, split like the 2,000 and 8,000: 200 exposed, 800 unexposed. The cross-product returns the risk ratio.
-
Concurrent controls
\[ \frac{600 \times 816}{960 \times 184} = 2.77 \]
Draw 1,000 controls in proportion to the 8,500 and 37,600 person-years: 184 exposed, 816 unexposed after rounding to whole people. The cross-product, 2.77, differs from the rate ratio of 2.76 only through that rounding.
Result: The same 1,560 cases give 3.14, 2.50 or 2.77, depending only on which denominator the controls stand for.
These pooled tables are illustrative. A real risk-set sample keeps its matched sets and is analysed set by set, as the section on concurrent controls shows.
The full cohort: four ratios as yardsticks
The table below gives the four full-cohort ratios, with 95% confidence intervals (CI) from the models in the code panes. They are the yardsticks for the three samples. Exclusive controls aim at the 5-year odds ratio, inclusive controls at the 5-year risk ratio and risk-set controls at the hazard ratio.
The hazard is the instantaneous event rate among people still free of the disease. A Cox model, the standard regression for time-to-event data, summarises the ratio of the two groups' hazards over follow-up as one hazard ratio. When each group's rate is constant over time, the hazard ratio and the rate ratio are the same quantity.
In this simulated cohort the exposed group's hazard rises over follow-up, while the unexposed group's stays almost flat. So the ratio of the two hazards itself drifts upward, and 2.78 is the Cox model's summary of it over follow-up. The rate ratio, 2.76, pools all cases and person-years instead, so the two summaries differ slightly. Ties play no part in that gap: no two cases share an event time, as the Stata header 'Cox regression with no ties' confirms.
| Measure (unit) | Exposed | Unexposed | Ratio (95% CI) |
|---|---|---|---|
| Risk (cases per person at the start) | 0.300 | 0.120 | 5-year risk ratio 2.50 (2.29 to 2.73) |
| Odds (cases per person still free of disease at the end) | 0.429 | 0.136 | 5-year odds ratio 3.14 (2.80 to 3.53) |
| Rate (cases per person-year) | 0.0706 | 0.0255 | Rate ratio 2.76 (2.50 to 3.06) |
| Hazard (Cox model) | Rises over follow-up | Almost flat over follow-up | Hazard ratio 2.78 (2.51 to 3.07) |
Stata: the full-cohort hazard ratio
* 5. Full-cohort hazard ratio: the target of risk-set sampling
stset time, failure(case) id(id)
stcox exposed, nolog. * 5. Full-cohort hazard ratio: the target of risk-set sampling
. stset time, failure(case) id(id)
Survival-time data settings
ID variable: id
Failure event: case!=0 & case<.
Observed time interval: (time[_n-1], time]
Exit on or before: failure
--------------------------------------------------------------------------
10,000 total observations
0 exclusions
--------------------------------------------------------------------------
10,000 observations remaining, representing
10,000 subjects
1,560 failures in single-failure-per-subject data
46,100 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 5
. stcox exposed, nolog
Failure _d: case
Analysis time _t: time
ID variable: id
Cox regression with no ties
No. of subjects = 10,000 Number of obs = 10,000
No. of failures = 1,560
Time at risk = 46,100
LR chi2(1) = 343.79
Log likelihood = -14067.77 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
exposed | 2.776434 .1445672 19.61 0.000 2.507066 3.074743
------------------------------------------------------------------------------R: the four full-cohort ratios
# 3. Full-cohort risk ratio (log-binomial) and odds ratio (logistic)
fit_rr <- glm(case ~ exposed, family = binomial(link = "log"), data = cohort)
canon_ratio("cohort.rr", coef(fit_rr), vcov(fit_rr))
fit_or <- glm(case ~ exposed, family = binomial, data = cohort)
canon_ratio("cohort.or", coef(fit_or), vcov(fit_or))
# 4. Full-cohort rate ratio: Poisson model with log person-years as the offset
fit_irr <- glm(case ~ exposed + offset(log(time)), family = poisson, data = cohort)
canon_ratio("cohort.irr", coef(fit_irr), vcov(fit_irr))
# 5. Full-cohort hazard ratio: the target of risk-set sampling (Breslow ties, as in Stata)
fit_hr <- coxph(Surv(time, case) ~ exposed, data = cohort, ties = "breslow")
canon_ratio("cohort.hr", coef(fit_hr), vcov(fit_hr))> fit_rr <- glm(case ~ exposed, family = binomial(link = "log"),
+ data = cohort)
> canon_ratio("cohort.rr", coef(fit_rr), vcov(fit_rr))
CANON w7.cohort.rr 2.5000
CANON w7.cohort.rr_lb 2.2861
CANON w7.cohort.rr_ub 2.7340
> fit_or <- glm(case ~ exposed, family = binomial, data = cohort)
> canon_ratio("cohort.or", coef(fit_or), vcov(fit_or))
CANON w7.cohort.or 3.1429
CANON w7.cohort.or_lb 2.7958
CANON w7.cohort.or_ub 3.5330
> fit_irr <- glm(case ~ exposed + offset(log(time)),
+ family = poisson, data = cohort)
> canon_ratio("cohort.irr", coef(fit_irr), vcov(fit_irr))
CANON w7.cohort.irr 2.7647
CANON w7.cohort.irr_lb 2.4966
CANON w7.cohort.irr_ub 3.0616
> fit_hr <- coxph(Surv(time, case) ~ exposed, data = cohort,
+ ties = "breslow")
> canon_ratio("cohort.hr", coef(fit_hr), vcov(fit_hr))
CANON w7.cohort.hr 2.7764
CANON w7.cohort.hr_lb 2.5071
CANON w7.cohort.hr_ub 3.0747Exclusive controls: people still free of disease at the end
Cumulative (exclusive) sampling draws controls only from people still free of the disease when follow-up ends. In this cohort that pool holds 8,440 people, 1,400 exposed and 7,040 unexposed. The controls' exposure odds estimate the exposure odds in that pool, so the cross-product estimates the cohort's 5-year odds ratio, 3.14 [2, 4]. That odds ratio approaches the risk ratio only when the disease is rare, and a risk of 0.300 in the exposed is far from rare.
The frozen sample drew 3,120 controls from the pool: 539 exposed and 2,581 unexposed. With the 1,560 cases, the cross-product is (600 × 2,581) / (960 × 539) = 2.99. In the sample files, y is 1 for a case and 0 for a control, and logistic fits
Here $\operatorname{logit}(p) = \ln\{p/(1 - p)\}$ is the log odds and $i$ indexes the rows. $\beta_0$ is the intercept and $\beta_1$ the exposure coefficient. With one binary exposure, $\exp(\beta_1)$ equals the cross-product, 2.99. Its 95% CI, 2.61 to 3.44, covers the target of 3.14.
Stata: the exclusive sample
cc y exposed
canon4 exclusive.or_lb_exact.stata r(lb_or)
canon4 exclusive.or_ub_exact.stata r(ub_or)
logistic y exposed, nolog. cc y exposed
Proportion
| Exposed Unexposed | Total exposed
-----------------+------------------------+------------------------
Cases | 600 960 | 1560 0.3846
Controls | 539 2581 | 3120 0.1728
-----------------+------------------------+------------------------
Total | 1139 3541 | 4680 0.2434
| |
| Point estimate | [95% conf. interval]
|------------------------+------------------------
Odds ratio | 2.992811 | 2.600961 3.443775 (exact)
Attr. frac. ex. | .6658659 | .6155268 .709621 (exact)
Attr. frac. pop | .2561023 |
+-------------------------------------------------
chi2(1) = 253.49 Pr>chi2 = 0.0000
. canon4 exclusive.or_lb_exact.stata r(lb_or)
CANON w7.exclusive.or_lb_exact.stata 2.6010
. canon4 exclusive.or_ub_exact.stata r(ub_or)
CANON w7.exclusive.or_ub_exact.stata 3.4438
. logistic y exposed, nolog
Logistic regression Number of obs = 4,680
LR chi2(1) = 243.62
Prob > chi2 = 0.0000
Log likelihood = -2857.0778 Pseudo R2 = 0.0409
------------------------------------------------------------------------------
y | Odds ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
exposed | 2.992811 .2105856 15.58 0.000 2.607267 3.435366
_cons | .3719489 .014061 -26.16 0.000 .345386 .4005546
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.Inclusive controls: everyone at the start
Case-base (inclusive) sampling draws controls from the whole cohort as it stood at the start, whatever happens to them later. The controls then mirror the 2,000 exposed and 8,000 unexposed people at baseline, the denominators of the two 5-year risks. The cross-product therefore estimates the 5-year risk ratio, 2.50 [2, 4].
Some controls become cases later, and they stay in. In the frozen sample, 493 of the 3,120 controls drawn from all 10,000 people became cases during follow-up. They appear twice, once in each role. The controls hold 619 exposed and 2,501 unexposed people, and the risk ratio (cross-product) is 2.53.
People who appear twice make the rows dependent, which affects the standard error but not the estimate. The script therefore clusters the standard error on the person: vce(cluster id) in Stata and vcovCL() from the sandwich package in R. The 95% CI, 2.25 to 2.83, covers the target of 2.50. Stata still heads the column 'Odds ratio', because the command computes a cross-product; the design is what makes it a risk ratio.
Keeping a random baseline sample as a subcohort and following it over time gives the case-cohort design, which has its own weighted analysis [5, 6]. This article analyses the inclusive sample only as a 2 by 2 table.
Stata: the inclusive sample
cc y exposed
* A person can appear as a case and as a control, so the standard error is clustered on the person
logistic y exposed, vce(cluster id) nolog. cc y exposed
Proportion
| Exposed Unexposed | Total exposed
-----------------+------------------------+------------------------
Cases | 600 960 | 1560 0.3846
Controls | 619 2501 | 3120 0.1984
-----------------+------------------------+------------------------
Total | 1219 3461 | 4680 0.2605
| |
| Point estimate | [95% conf. interval]
|------------------------+------------------------
Odds ratio | 2.525242 | 2.201884 2.896123 (exact)
Attr. frac. ex. | .6039984 | .5458434 .6547108 (exact)
Attr. frac. pop | .2323071 |
+-------------------------------------------------
chi2(1) = 187.22 Pr>chi2 = 0.0000
. * A person can appear as a case and as a control, so the standard error is cl
> ustered on the person
. logistic y exposed, vce(cluster id) nolog
Logistic regression Number of obs = 4,680
Wald chi2(1) = 250.45
Prob > chi2 = 0.0000
Log pseudolikelihood = -2888.3749 Pseudo R2 = 0.0304
(Std. err. adjusted for 4,187 clusters in id)
------------------------------------------------------------------------------
| Robust
y | Odds ratio std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
exposed | 2.525242 .1478116 15.83 0.000 2.251537 2.83222
_cons | .3838465 .0132611 -27.72 0.000 .3587156 .4107379
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.Concurrent controls: the risk set at each case's event time
Risk-set (concurrent) sampling works through the cases in time order. When a case occurs, controls are drawn from its risk set, everyone still under follow-up and free of the disease at that moment. The case and its controls form a matched set. A control may become a case later, and one person may be drawn for several sets.
People who stay free of the disease longer sit in more risk sets, so the controls, taken together, roughly mirror the person-time at risk. The matched analysis, not a pooled table, carries the estimate. Analysed set by set, the design estimates the hazard ratio without a rare-disease assumption [3, 7].
In this cohort the target is the full-cohort Cox hazard ratio, 2.78. The matched analysis aims at exactly that value when the ratio of the hazards stays constant over time, called proportional hazards. When the ratio drifts, as here, the matched analysis aims at essentially the same summary as the full-cohort Cox model.
Conditional logistic regression is that set-by-set analysis. It is a logistic model with its own intercept for every matched set:
$$\operatorname{logit} \Pr(y_{ij} = 1) = \alpha_j + \beta\,\mathrm{exposed}_{ij}$$Here $j$ indexes the matched sets, $i$ the people within a set, and $y_{ij}$ is 1 for the case. The term $\alpha_j$ absorbs everything the members of set $j$ share, including the time of the case. Conditioning on exactly one case per set removes every $\alpha_j$, so each case is compared only with its own controls. What remains has the same form as a Cox model's contribution from that risk set, so $\exp(\beta)$ estimates the hazard ratio [7].
The frozen risk-set sample has 1,560 matched sets of one case and two controls, with 609 exposed and 2,511 unexposed controls in all. Of those controls, 244 became cases later. clogit y exposed, group(set) or nolog, where set numbers the matched sets, gives a hazard ratio of 2.61 (95% CI 2.27 to 3.01), covering the target of 2.78. R's clogit() on the same sets gives the same value.
Stata: the frozen risk-set sample
clogit y exposed, group(set) or nolog
canon4 riskset.hr exp(_b[exposed])
canon4 riskset.hr_lb exp(_b[exposed]-invnormal(0.975)*_se[exposed])
canon4 riskset.hr_ub exp(_b[exposed]+invnormal(0.975)*_se[exposed])
* For contrast only: pooling the risk sets into one unmatched 2 by 2 table ignores the sampling design
logistic y exposed, nolog. clogit y exposed, group(set) or nolog
Conditional (fixed-effects) logistic regression Number of obs = 4,680
LR chi2(1) = 188.91
Prob > chi2 = 0.0000
Log likelihood = -1619.3805 Pseudo R2 = 0.0551
------------------------------------------------------------------------------
y | Odds ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
exposed | 2.613864 .186032 13.50 0.000 2.273536 3.005136
------------------------------------------------------------------------------
. canon4 riskset.hr exp(_b[exposed])
CANON w7.riskset.hr 2.6139
. canon4 riskset.hr_lb exp(_b[exposed]-invnormal(0.975)*_se[exposed])
CANON w7.riskset.hr_lb 2.2735
. canon4 riskset.hr_ub exp(_b[exposed]+invnormal(0.975)*_se[exposed])
CANON w7.riskset.hr_ub 3.0051
. * For contrast only: pooling the risk sets into one unmatched 2 by 2 table ig
> nores the sampling design
. logistic y exposed, nolog
Logistic regression Number of obs = 4,680
LR chi2(1) = 188.17
Prob > chi2 = 0.0000
Log likelihood = -2884.8011 Pseudo R2 = 0.0316
------------------------------------------------------------------------------
y | Odds ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
exposed | 2.57697 .1775796 13.74 0.000 2.251402 2.949619
_cons | .3823178 .0145075 -25.34 0.000 .3549152 .4118361
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.Stata's own draw: stset, sttocc and clogit
Stata can draw the risk sets itself, in three steps. First, stset time, failure(case) id(id) declares the follow-up time, the event and the person. Then sttocc exposed, number(2) nodots draws two controls from each case's risk set and keeps the variable exposed in the new data. Finally, clogit _case exposed, group(_set) or nolog analyses the matched sets.
sttocc adds three variables: _case (1 for the case, 0 for a control), _set (the matched set) and _time (the case's event time). The analysis must use _case and _set, because the cohort's own case indicator would mark a control who later becomes a case as a case. Stata's draw gives a hazard ratio of 2.64 (95% CI 2.30 to 3.03). It differs from the frozen sample's 2.61 because it is a different random draw from the same risk sets.
Stata: sttocc, then clogit
* 6. Stata's own risk-set draw: sttocc keeps each case and samples 2 controls from its risk set
sttocc exposed, number(2) nodots
count
display "CANON w7.sttocc.n.stata " r(N)
tabulate _case exposed
* Conditional logistic regression keeps each matched set intact and estimates the hazard ratio
clogit _case exposed, group(_set) or nolog. tabulate _case exposed
0 for |
controls; |
1 for | exposed
cases | 0 1 | Total
-----------+----------------------+----------
0 | 2,534 586 | 3,120
1 | 960 600 | 1,560
-----------+----------------------+----------
Total | 3,494 1,186 | 4,680
. * Conditional logistic regression keeps each matched set intact and estimates
> the hazard ratio
. clogit _case exposed, group(_set) or nolog
Conditional (fixed-effects) logistic regression Number of obs = 4,680
LR chi2(1) = 198.56
Prob > chi2 = 0.0000
Log likelihood = -1614.5561 Pseudo R2 = 0.0579
------------------------------------------------------------------------------
_case | Odds ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
exposed | 2.642912 .1856319 13.84 0.000 2.303012 3.032977
------------------------------------------------------------------------------Three control series side by side
| Control series | Controls drawn from | What the analysis estimates | Target (simulated truth) | Frozen sample (95% CI) | Stata | R |
|---|---|---|---|---|---|---|
| Cumulative (exclusive) | The 8,440 people free of disease at the end | Odds ratio | 3.14 | Odds ratio 2.99 (2.61 to 3.44) | cc y exposed; logistic y exposed, nolog | glm(y ~ exposed, family = binomial, data = ex) |
| Case-base (inclusive) | All 10,000 people at the start | Risk ratio | 2.50 | Risk ratio (cross-product) 2.53 (2.25 to 2.83) | logistic y exposed, vce(cluster id) nolog | glm(y ~ exposed, family = binomial, data = inc), then vcovCL() clustered on id |
| Risk-set (concurrent) | People still at risk at each case's event time | Hazard ratio (the rate ratio when rates are constant) | 2.78 | Hazard ratio 2.61 (2.27 to 3.01) | clogit y exposed, group(set) or nolog | clogit(y ~ exposed + strata(set), data = rs) |
R: the three control series
# 6. Exclusive controls (non-cases at the end of follow-up): the cross-product estimates the odds ratio
ex <- dat("W7_exclusive.csv")
canon_n("exclusive.n", nrow(ex))
canon_n("exclusive.controls_exposed", sum(ex$y == 0 & ex$exposed == 1))
canon_n("exclusive.controls_unexposed", sum(ex$y == 0 & ex$exposed == 0))
tab_ex <- table(exposed = factor(ex$exposed, 1:0), y = factor(ex$y, 1:0))
print(tab_ex)
fit_ex <- glm(y ~ exposed, family = binomial, data = ex)
canon_ratio("exclusive.or", coef(fit_ex), vcov(fit_ex))
# Exact limits from the noncentral hypergeometric distribution (point value here is the conditional estimate)
fx <- fisher.test(tab_ex)
canon("exclusive.or_cmle.r", unname(fx$estimate))
canon("exclusive.or_lb_exact.r", fx$conf.int[1])
canon("exclusive.or_ub_exact.r", fx$conf.int[2])
# 7. Inclusive controls (everyone at baseline): the cross-product estimates the risk ratio
inc <- dat("W7_inclusive.csv")
canon_n("inclusive.n", nrow(inc))
canon_n("inclusive.controls_exposed", sum(inc$y == 0 & inc$exposed == 1))
canon_n("inclusive.controls_unexposed", sum(inc$y == 0 & inc$exposed == 0))
canon_n("inclusive.controls_later_cases", sum(inc$y == 0 & inc$case == 1))
fit_inc <- glm(y ~ exposed, family = binomial, data = inc)
# A person can appear as a case and as a control, so the variance is clustered on the person
# (HC0 with the G/(G-1) cluster factor, the same small-sample factor Stata's logistic uses)
v_inc <- vcovCL(fit_inc, cluster = ~id, type = "HC0", cadjust = TRUE)
canon_ratio("inclusive.or", coef(fit_inc), v_inc)
# 8. Concurrent (risk-set) controls, frozen draw shared with Stata: matched analysis by set
rs <- dat("W7_riskset.csv")
canon_n("riskset.n", nrow(rs))
canon_n("riskset.sets", length(unique(rs$set)))
canon_n("riskset.controls_exposed", sum(rs$y == 0 & rs$exposed == 1))
canon_n("riskset.controls_unexposed", sum(rs$y == 0 & rs$exposed == 0))
canon_n("riskset.controls_later_cases", sum(rs$y == 0 & rs$case == 1))
fit_rs <- clogit(y ~ exposed + strata(set), data = rs)
canon_ratio("riskset.hr", coef(fit_rs), vcov(fit_rs))
# For contrast only: pooling the risk sets into one unmatched 2 by 2 table ignores the sampling design
fit_crude <- glm(y ~ exposed, family = binomial, data = rs)
canon("riskset.crude_or", exp(coef(fit_crude)[["exposed"]])) > ex <- dat("W7_exclusive.csv")
> canon_n("exclusive.n", nrow(ex))
CANON w7.exclusive.n 4680
> canon_n("exclusive.controls_exposed", sum(ex$y ==
+ 0 & ex$exposed == 1))
CANON w7.exclusive.controls_exposed 539
> canon_n("exclusive.controls_unexposed", sum(ex$y ==
+ 0 & ex$exposed == 0))
CANON w7.exclusive.controls_unexposed 2581
> tab_ex <- table(exposed = factor(ex$exposed, 1:0),
+ y = factor(ex$y, 1:0))
> print(tab_ex)
y
exposed 1 0
1 600 539
0 960 2581
> fit_ex <- glm(y ~ exposed, family = binomial, data = ex)
> canon_ratio("exclusive.or", coef(fit_ex), vcov(fit_ex))
CANON w7.exclusive.or 2.9928
CANON w7.exclusive.or_lb 2.6073
CANON w7.exclusive.or_ub 3.4354
> fx <- fisher.test(tab_ex)
> canon("exclusive.or_cmle.r", unname(fx$estimate))
CANON w7.exclusive.or_cmle.r 2.9920
> canon("exclusive.or_lb_exact.r", fx$conf.int[1])
CANON w7.exclusive.or_lb_exact.r 2.6007
> canon("exclusive.or_ub_exact.r", fx$conf.int[2])
CANON w7.exclusive.or_ub_exact.r 3.4436
> inc <- dat("W7_inclusive.csv")
> canon_n("inclusive.n", nrow(inc))
CANON w7.inclusive.n 4680
> canon_n("inclusive.controls_exposed", sum(inc$y ==
+ 0 & inc$exposed == 1))
CANON w7.inclusive.controls_exposed 619
> canon_n("inclusive.controls_unexposed", sum(inc$y ==
+ 0 & inc$exposed == 0))
CANON w7.inclusive.controls_unexposed 2501
> canon_n("inclusive.controls_later_cases", sum(inc$y ==
+ 0 & inc$case == 1))
CANON w7.inclusive.controls_later_cases 493
> fit_inc <- glm(y ~ exposed, family = binomial, data = inc)
> v_inc <- vcovCL(fit_inc, cluster = ~id, type = "HC0",
+ cadjust = TRUE)
> canon_ratio("inclusive.or", coef(fit_inc), v_inc)
CANON w7.inclusive.or 2.5252
CANON w7.inclusive.or_lb 2.2515
CANON w7.inclusive.or_ub 2.8322
> rs <- dat("W7_riskset.csv")
> canon_n("riskset.n", nrow(rs))
CANON w7.riskset.n 4680
> canon_n("riskset.sets", length(unique(rs$set)))
CANON w7.riskset.sets 1560
> canon_n("riskset.controls_exposed", sum(rs$y == 0 &
+ rs$exposed == 1))
CANON w7.riskset.controls_exposed 609
> canon_n("riskset.controls_unexposed", sum(rs$y ==
+ 0 & rs$exposed == 0))
CANON w7.riskset.controls_unexposed 2511
> canon_n("riskset.controls_later_cases", sum(rs$y ==
+ 0 & rs$case == 1))
CANON w7.riskset.controls_later_cases 244
> fit_rs <- clogit(y ~ exposed + strata(set), data = rs)
> canon_ratio("riskset.hr", coef(fit_rs), vcov(fit_rs))
CANON w7.riskset.hr 2.6139
CANON w7.riskset.hr_lb 2.2735
CANON w7.riskset.hr_ub 3.0051
> fit_crude <- glm(y ~ exposed, family = binomial, data = rs)
> canon("riskset.crude_or", exp(coef(fit_crude)[["exposed"]]))
CANON w7.riskset.crude_or 2.5770 Common misreadings and their fixes
-
Dropping inclusive controls who later become cases
In the frozen inclusive sample, 493 of the 3,120 controls became cases later. Removing them leaves only controls who stayed free of the disease, the exclusive pool, and pushes the cross-product from the risk ratio toward the odds ratio.
Fix: Keep every baseline control, even when the same person also appears as a case, and cluster the standard error on the person.
-
Analysing a risk-set sample with one pooled, unmatched logistic model
Pooling the sets into one 2 by 2 table throws away when each control was drawn. In this draw the pooled value, 2.58, happens to land near the matched 2.61. That is a feature of this cohort and this draw, not a guarantee. A pooled unmatched odds ratio is not generally valid for a risk-set sample.
Fix: Analyse the matched sets with conditional logistic regression,
clogitin both Stata and R. -
Running
ccon matched dataccbuilds one unmatched 2 by 2 table. On a risk-set sample it returns the same pooled cross-product as the unmatched logistic model, with the same problem.Fix: Use
ccorlogisticfor an unmatched exclusive sample,logisticwithvce(cluster id)for an inclusive sample (cccannot cluster), andclogitfor matched sets. -
Analysing sttocc output with the cohort's own variables
After
sttocc, the case of each matched set is flagged by_caseand the set by_set. The cohort's case indicator would flag a control who becomes a case later, and a model withoutgroup(_set)ignores the matching.Fix: Follow
sttocc exposed, number(2) nodotswithclogit _case exposed, group(_set) or nolog. -
"Every case-control odds ratio estimates the odds ratio."
The same 1,560 cases gave an odds ratio of 2.99, a risk ratio (cross-product) of 2.53 and a hazard ratio of 2.61. Their targets were 3.14, 2.50 and 2.78. Stata labels each of these estimates an odds ratio, and in R the exponentiated logistic or clogit coefficient is conventionally read as one, whatever the design. Many published studies leave the sampling rule unstated [1].
Fix: What the cross-product estimates depends on how controls were sampled: the odds ratio for controls from end-of-follow-up non-cases, the risk ratio for controls from the baseline cohort, and the rate or hazard ratio for controls from the risk set at each case's event time.
What to do in your own analysis
- Write the control-sampling rule into the protocol: the pool, the time of sampling, the number of controls per case, whether a control may later become a case, and whether everyone was followed for the same period without loss.
- Name the measure that rule and that follow-up target (odds ratio, risk ratio or hazard ratio) before any model is fitted.
- For a risk-set sample, keep the matched sets and analyse them with conditional logistic regression; a pooled table can describe the sample but is usually not the estimate to report.
- For an inclusive sample, keep the controls who later become cases; the standard error is then usually clustered on the person.
- Report each estimate under the name of the measure it estimates, whatever label the software prints, and state the sampling rule beside it.
Glossary
- case-control study
- A study that measures the exposure in the cases and in a sample of controls drawn to represent their source population, instead of in everyone.
- control series
- The controls a study draws, together with the rule that drew them.
- closed cohort
- A cohort in which everyone enters at the start and is followed for the same fixed period with no loss to follow-up.
- cumulative (exclusive) sampling (การเลือกกลุ่มควบคุมแบบสะสม)
- Drawing controls from people still free of the disease at the end of follow-up; in a closed cohort, the cross-product estimates the odds ratio.
- case-base (inclusive) sampling (การเลือกกลุ่มควบคุมจากฐานประชากรทั้งหมด)
- Drawing controls from the whole cohort at the start, including people who later become cases; in a closed cohort, the cross-product estimates the risk ratio.
- risk-set (concurrent) sampling (การสุ่มกลุ่มควบคุมจากกลุ่มเสี่ยง)
- Drawing controls at each case's event time from the people still at risk; analysed as matched sets, it estimates the hazard ratio.
- risk set
- Everyone still under follow-up and free of the disease when a case occurs.
- nested case-control design (การศึกษาแบบ nested case-control)
- A case-control study drawn from a defined cohort; in standard usage, with controls drawn from each case's risk set.
- cross-product
- Exposed cases times unexposed controls, divided by unexposed cases times exposed controls; the number cc and logistic print as an odds ratio.
- conditional logistic regression
- A logistic model with its own intercept for each matched set, fitted so that each case is compared only with its own controls.
- matched set
- One case together with the controls drawn for it.
- hazard ratio (อัตราส่วนฮาซาร์ด)
- The ratio of instantaneous event rates between two groups among people still free of the event.
References
- Knol MJ, Vandenbroucke JP, Scott P, Egger M. What do case-control studies estimate? Survey of methods and assumptions in published case-control research. Am J Epidemiol. 2008;168(9):1073-1081. doi:10.1093/aje/kwn217 https://doi.org/10.1093/aje/kwn217
- Rodrigues L, Kirkwood BR. Case-control designs in the study of common diseases: updates on the demise of the rare disease assumption and the choice of sampling scheme for controls. Int J Epidemiol. 1990;19(1):205-213. doi:10.1093/ije/19.1.205 https://doi.org/10.1093/ije/19.1.205
- Greenland S, Thomas DC. On the need for the rare disease assumption in case-control studies. Am J Epidemiol. 1982;116(3):547-553. doi:10.1093/oxfordjournals.aje.a113439 https://doi.org/10.1093/oxfordjournals.aje.a113439
- Pearce N. What does the odds ratio estimate in a case-control study? Int J Epidemiol. 1993;22(6):1189-1192. doi:10.1093/ije/22.6.1189 https://doi.org/10.1093/ije/22.6.1189
- Prentice RL. A case-cohort design for epidemiologic cohort studies and disease prevention trials. Biometrika. 1986;73(1):1-11. doi:10.1093/biomet/73.1.1 https://doi.org/10.1093/biomet/73.1.1
- Barlow WE, Ichikawa L, Rosner D, Izumi S. Analysis of case-cohort designs. J Clin Epidemiol. 1999;52(12):1165-1172. doi:10.1016/S0895-4356(99)00102-X https://doi.org/10.1016/S0895-4356(99)00102-X
- Langholz B, Goldstein L. Risk set sampling in epidemiologic cohort studies. Stat Sci. 1996;11(1):35-53. doi:10.1214/ss/1032209663 https://doi.org/10.1214/ss/1032209663
Key takeaways
- The cases are the same in every design; the controls decide what the cross-product estimates, because they stand in for a denominator.
- In a closed cohort, controls from people still free of disease at the end estimate the 5-year odds ratio, 3.14 here, which approaches the risk ratio only when the disease is rare.
- In a closed cohort, controls from the whole cohort at the start, including those who later become cases, estimate the 5-year risk ratio, 2.50 here; the standard error is then usually clustered on the person.
- Controls from each case's risk set, analysed as matched sets with conditional logistic regression, estimate the hazard ratio, 2.78 here, with no need for a closed cohort; a pooled unmatched odds ratio is not generally valid.
- Report each estimate under the name of the measure it estimates, whatever label Stata or R prints.
Related in the wiki: [[nested-case-control-design]]