Skip to content
Licensed Unlicensed Requires Authentication Published by De Gruyter January 15, 2019

Hazard Ratio Estimators after Terminating Observation within Matched Pairs in Sibling and Propensity Score Matched Designs

  • Tomohiro Shinozaki ORCID logo EMAIL logo and Mohammad Ali Mansournia


Similar to unmatched cohort studies, matched cohort studies may suffer from the censoring of events prior to the end of follow-up. Moreover, in some matched-pair cohort studies, observation time is prematurely terminated immediately after the follow-up of his/her matched member is completed by an event or censoring. Although the follow-up termination within matched pairs may or may not change the hazard ratio estimators, when and how the change occurs has not been clarified. We study the change in the estimates of the hazard ratio conditional on matched pairs and/or covariates by considering two types of matched-pair designs in cohort studies—sibling pair matching and propensity score matching—in which termination can be naturally considered. If all possible confounders are shared within the matched pairs, after termination, a wide range of hazard ratio estimators coincides with that obtained from a stratified Cox model. If unshared confounders should be adjusted for in the analysis, however, such coincidence is not observed. Simulation studies on sibling designs with unshared confounders suggested that the pair-stratified covariate-adjusted Cox model for the hazard ratio conditional on matched pairs and covariates is generally preferred, for which termination does not deteriorate the estimation. Conversely, the comparison between stratifying or not stratifying on pair is a more subtle issue in propensity score matching which targets a marginal or covariate-conditional hazard ratio. Based on simulation studies considering Cox models after matching based on estimated propensity scores, we discourage pair-stratified analysis and termination, particularly after data collection.

Funding statement: This work was supported by the Japan Society for the Promotion of Science (Funder Id: 10.13039/501100001691, Grant-in-Aid for Young Scientists B 16K16015).


A Random-effect models without covariates fitted to terminated data (Section 3.2)

Maximum likelihood estimates of random-effect models are obtained for joint likelihood of data and random effects marginalized over random effects. For random-effect Cox models without covariates, the marginal (with respect to α = (α1, …, αn)) partial likelihood is


where R(Yi*) is a risk set (of pairs) at time Yi* and f(α;θ) is a user-specified parametric distribution of random effects α. Estimation of β is only relevant to the components outside the integral, which is equal to the partial likelihood of the stratified Cox model without covariates (3) [8].

For random-effect Poisson models without covariates, marginal likelihood is


where g(α, β, θ) = ∑iαi(di0* + di1*) – Yi*exp(αi)(1 + eβ) + log f(α;θ). From generalized linear mixed models theory, random effects αi should be subject to ∂g(α,β,θ)/∂αi = (di0* + di1*) – Yi*exp(αi)(1 + eβ) + Qi = 0 for all i at any β (where Qi = ∂log f(α,θ)/∂αi is a penalty term): it implies exp(αi) = {Yi*(1 + β)}–1(di0* + di1* + Qi). Rearranging the marginal likelihood yields


The estimation of β is only relevant to the parts outside the integral, which is equal to the partial likelihood of the stratified Cox model without covariates (3).

B Conditional models with shared covariates fitted to terminated data (Section 3.3)

For Cox models conditional on shared covariates Xi, partial likelihood is


where R(Yi*) is a risk set (of pairs) at time Yi*. Note that di1*di0* = 0. The contribution for the partial likelihood from pair i can be decomposed into eβ/(1 +eβ)di1*, 1/(1 +eβ)di0*, and exp(γTXi)/kR(Yi)exp(γTXi)di1+di0; from the first two factors of the likelihood including β, score equation of β is ∑i{di1*di0*eβ} = 0, which is independent of γ. The partial likelihood including β is the same as that from the model (3), so are its first and second derivatives.

For Poisson models conditional on shared covariates Xi, the score equations are


The first two equations jointly imply ∑i{di0*Yi*exp(α + γTXi)} = 0 and ∑i {di1*Yi*exp(α + γTXi)eβ} = 0; hence a maximum likelihood estimate of β is the solution of ∑i {di1*di0*eβ} = 0, which is independent of α and γ.

β is estimated independent of γ in both equations; therefore, robust variance of β is estimated as {∑i mi(β)/∂β}–1i mi(β)2{∑imi(β)/∂β}–1, where mi(β) = di1*di0*eβ and β is substituted by its estimate. After some algebra, this robust variance estimator becomes 1/i= 1ndi0* + 1/i= 1ndi1*, which is also the β-component of the inverse Fisher information in both Cox and Poisson models considered here (the latter calculation is somewhat complex and deferred to Appendix C). This is true whether Xi are adjusted for or not.

C Variance estimates based on Fisher information in a Poisson model conditional on shared covariates with terminated data (Section 3.3)

For simplicity, assume shared covariates Xi as a scalar Xi (though the following result holds if Xi is a vector). As shown in the Appendix B, the score equations for the Poisson model conditional on shared covariates Xi are


where (Sα, Sβ, Sγ)T is the score function of (α, β, γ)T, i. e. the first derivative of the log-likelihood with respect to (α, β, γ)T. To derive the variance of the maximum likelihood estimator of β, we reorder the score as (Sβ, Sα, Sγ)T. Further differentiating it with (β, α, γ) yields the negative of the Fisher information matrix ∂(Sβ, Sα, Sγ)T/∂(β, α, γ) as


We partition matrix (13) as


Following basic matrix rules, the (1, 1) part of the inverse of matrix (13) is obtained as (A – BD–1C)–1. Denoting K =[(1+eβ)2idi0iXi2Yi(di0+di1)iXi(di0+di1)2]1,


Substituting maximum likelihood estimator exp(βˆ)=idi1/idi0, (A – BD–1C)–1 is equal to idi1(idi1)2idi1+idi01=1/idi0+1/idi1, as desired.

D Fixed-effect models fitted to terminated data (Section 4.2)

Fixed-effect Cox models’ partial likelihood is


Equating the first derivative of the log-partial likelihood with respect to αi (i = 1, …, n) to 0, fitting the fixed-effect Cox model imposes the following conditions for all i = 1, …, n: kR(Yi),kieαkexp(β+γTXk1)+exp(γTXk0)=0. Thus, the partial likelihood reduces to


which is the partial likelihood (2) of the stratified Cox model (1) by applying Result 1:


Fixed-effect Poisson models’ score equations are further augmented with

jin pairidijYiexp(αi+βZij+γTXij)=0fori=1,,n.

Substituting Yi*exp(αi) = (di0* + di1*)/{exp(β + γTXi1) + exp(γTXi0)} from the above conditions, ∑ijZij{dij*Yi*exp(αi + βZij + γTXij)} and ∑ijXij{dij*Yi*exp(αi + βZij + γTXij)} reduce to the first derivative of log-partial likelihood (2) with respect to β and γ, respectively. To be specific, score equations of the fixed-effect Poisson model are

Sαi=jin pairidijYiexp(αi+βZij+γTXij)=0fori=1,,n,Sβ=ijZijdijYiexp(αi+βZij+γTXij)=0,Sγ=ijXijdijYiexp(αi+βZij+γTXij)=0.

The first equation of Sαi=0 imply Yi*exp(αi) = (di0* + di1*)/{exp(β + γTXi1) + exp(γTXi0)}.

Deleting Yi*exp(αi) from the second equation, Sβ = ∑i [di1* – (di0* + di1*)exp(β + γTXi1)/{exp(β + γTXi1) + exp(γTXi0)}] = ∑i [{di1*exp(γTXi0) + di0* exp(β + γTXi1)}/{exp(β + γTXi1) + exp(γTXi0)}]. This is the same as the first derivative of the log-partial likelihood of the stratified Cox model (1)


with respect to β.

Similarly, deleting Yi*exp(αi) from the third equation, Sγ = ∑i [Xi1{di1*exp(β + γTXi1) + di1*exp(γTXi0) – di1*exp(β + γTXi1) – di0*exp(β + γTXi1)}/{exp(β + γTXi1) + exp(γTXi0)}] + ∑i [Xi0{di0*exp(β + γTXi1) + di0*exp(γTXi0) – di1*exp(γTXi0) – di0*exp(γTXi0)}/{exp(β + γTXi1) + exp(γTXi0)}] = ∑i [{Xi1di1*exp(γTXi0) – Xi1di0*exp(β + γTXi1) + Xi0di0*exp(β + γTXi1) – Xi0di1*exp(γTXi0)}/{exp(β + γTXi1) + exp(γTXi0)}]. As before, the first derivative of the log-partial likelihood of the stratified Cox model (1) with respect to γ provides the same function.


We are grateful to Dr Takahiro Tabuchi (Osaka International Cancer Institute, Japan) for discussing a statistical analysis plan of the Longitudinal Survey of Middle-aged and Elderly Persons.


[1] Rothman KJ, Greenland S, Lash TL. Design strategies to improve study accuracy. In: Greenland S, Rothman KJ, Lash TL, editors. Modern epidemiology. 3rd ed. Philadelphia, PA: Lippincott Williams and Wilkins, 2008: 168–82.Search in Google Scholar

[2] Mansournia MA, Hernán MA, Greenland S. Matched designs and causal diagrams. Int J Epidemiol. 2013;42:860–9.10.1093/ije/dyt083Search in Google Scholar PubMed PubMed Central

[3] Sjölander A, Greenland S. Ignoring the matching variables in cohort studies: when is it valid and why?. Stat Med. 2013;32:4696–708.10.1002/sim.5879Search in Google Scholar PubMed

[4] Sjölander A, Zetterqvist J. Confounders, mediators, or colliders: what types of shared covariates does a sibling comparison design control for?. Epidemiology. 2017;28:540–7.10.1097/EDE.0000000000000649Search in Google Scholar PubMed

[5] Sjölander A, Frisell T, Kuja-Halkola R, Öberg S, Zetterqvist J. Carryover effects in sibling comparison designs. Epidemiology. 2016;27:852–8.10.1097/EDE.0000000000000541Search in Google Scholar PubMed

[6] Sjölander A, Johansson ALV, Lundholm C, Altman D, Almqvist C, Pawitan Y. Analysis of 1:1 matched cohort studies and twin studies, with binary exposures and binary outcomes. Stat Sci. 2012;27:395–411.10.1214/12-STS390Search in Google Scholar

[7] Holt JD, Prentice RL. Survival analyses in twin studies and matched pair experiments. Biometrika. 1974;61:17–30.10.1093/biomet/61.1.17Search in Google Scholar

[8] Shinozaki T, Mansournia MA, Matsuyama Y. On hazard ratio estimators by proportional hazards models in matched-pair cohort studies. Emerg Themes Epidemiol. 2017;14:6.10.1186/s12982-017-0060-8Search in Google Scholar PubMed PubMed Central

[9] Sjölander A, Lichtenstein P, Larsson H, Pawitan Y. Between-within models for survival analysis. Stat Med. 2013;32:3067–76.10.1002/sim.5767Search in Google Scholar PubMed

[10] Austin PC. A critical appraisal of propensity-score matching in the medical literature between 1996 and 2003 (with Discussion). Stat Med. 2008;27:2037–69.10.1002/sim.3243Search in Google Scholar

[11] Rubin DB. Matching to remove bias in observational studies. Biometrics. 1973;29:159–84.10.1017/CBO9780511810725.007Search in Google Scholar

[12] Rubin DB. The use of matched sampling and regression adjustment to remove bias in observational studies. Biometrics. 1973;29:185–203.10.1017/CBO9780511810725.008Search in Google Scholar

[13] Stuart EA. Matching methods for causal inference: a review and a look forward. Stat Sci. 2010;25:1–21.10.1214/09-STS313Search in Google Scholar PubMed PubMed Central

[14] King G, Nielsen R Why propensity scores should not be used for matching. Copy at Export BibTex tagged XML download paper, 481, 2016.Search in Google Scholar

[15] Sutradhar R, Baxter NN, Austin PC. Terminating observation within matched pairs of subjects in a matched cohort analysis: a Monte Carlo simulation study. Stat Med. 2016;35:294–304.10.1002/sim.6621Search in Google Scholar PubMed

[16] Richardson DP, Sutradhar R, Daly C, Paszat LF, Wilton AS, Rabeneck L, et al. Hospitalization rates in survivors of young adult malignancies. J Clin Oncol. 2014;33:2655–9.10.1200/JCO.2014.60.1914Search in Google Scholar PubMed

[17] Cronin-Fenton DP, Antonsen S, Cetin K, Acquavella J, Daniels A, Lash TL. Methods and rationale used in a matched cohort study of the incidence of new primary cancers following prostate cancer. Clin Epidemiol. 2013;5:429–37.10.2147/CLEP.S49713Search in Google Scholar PubMed PubMed Central

[18] Oshio T. The association between involvement in family caregiving and mental health among middle-aged adults in Japan. Soc Sci Med. 2014;115:121–9.10.1016/j.socscimed.2014.06.016Search in Google Scholar PubMed

[19] Tabuchi T, Fujiwara T, Shinozaki T. Tobacco price increase and smoking behaviour changes in various subgroups: a nationwide longitudinal 7-year follow-up study among a middle-aged Japanese population. Tob Control. 2017;26:69–77.10.1136/tobaccocontrol-2015-052804Search in Google Scholar PubMed PubMed Central

[20] Allison P. Fixed effects regression models, quantitative applications in the social sciences, Volume 160. Los Angeles: SAGE, 2009.10.4135/9781412993869Search in Google Scholar

[21] Lesko CR, Edwards JK, Cole SR, Moore RD, Lau B. When to censor? Am J Epidemiol. 2018;187:623–32.10.1093/aje/kwx281Search in Google Scholar PubMed PubMed Central

[22] Hernán MA, Robins JM. Causal inference. Boca Raton: Chapman & Hall/CRC, Forthcoming.Search in Google Scholar

[23] Rosenbaum PR, Rubin DB. The central role of the propensity score in observational studies for causal effects. Biometrika. 1983;70:41–55.10.21236/ADA114514Search in Google Scholar

[24] Greenland S. Introduction to categorical statistics. In: Greenland S, Rothman KJ, Lash TL, editors. Modern epidemiology, 3rd ed. Philadelphia, PA: Lippincott Williams and Wilkins, 2008: 238–57.Search in Google Scholar

[25] Greenland S, Robins JM. Estimation of a common effect parameter from sparse follow-up data. Biometrics. 1985;41:55–68.10.2307/2530643Search in Google Scholar

[26] Greenland S. Application of stratified analysis methods. In: Greenland S, Rothman KJ, Lash TL, editors. Modern epidemiology, 3rd ed. Philadelphia, PA: Lippincott Williams and Wilkins, 2008: 283–302.Search in Google Scholar

[27] Cummings P, McKnight B, Greenland S. Matched cohort methods for injury research. Epidemiol Rev. 2003;25:43–50.10.1093/epirev/mxg002Search in Google Scholar PubMed

[28] Richardson DB, Langholz B. Background stratified Poisson regression analysis of cohort data. Radiat Environ Biophys. 2012;51:15–22.10.1007/s00411-011-0394-5Search in Google Scholar PubMed PubMed Central

[29] Austin PC. A comparison of 12 algorithms for matching on the propensity score. Stat Med. 2014;33:1057–69.10.1002/sim.6004Search in Google Scholar PubMed PubMed Central

[30] Austin PC. The use of propensity score methods with survival or time-to-event outcomes: reporting measures of effect similar to those used in randomized experiments. Stat Med. 2014;33:1242–58.10.1002/sim.5984Search in Google Scholar PubMed PubMed Central

[31] Li L, Greene T. A weighting analogue to pair matching in propensity score analysis. Int J Biostat. 2013;9:215–34.10.1515/ijb-2012-0030Search in Google Scholar PubMed

[32] Martens EP, Pestman WR, Klungel OH. Re: conditioning on the propensity score can result in biased estimation of common measures of treatment effect: a Monte Carlo study. Stat Med. 2007;26:3208–10.10.1002/sim.2878Search in Google Scholar PubMed

[33] Greenland S, Mansournia MA, Altman DG. Sparse data bias: a problem hiding in plain sight. Br Med J. 2016;352:i1981.10.1136/bmj.i1981Search in Google Scholar PubMed

[34] Mansournia MA, Jewell NP, Greenland S. Case-control matching: effects, misconceptions, and recommendations. Eur J Epidemiol. 2018;33:5–14.10.1007/s10654-017-0325-0Search in Google Scholar PubMed

[35] Greenland S, Jewell NP, Mansournia MA. Theory and methodology: essential tools that can become dangerous belief systems. Eur J Epidemiol. 2018;33:503–6.10.1007/s10654-018-0395-7Search in Google Scholar PubMed PubMed Central

[36] Greenland S, Mansournia MA. Penalization, bias reduction, and default priors in logistic and related categorical and survival regressions. Stat Med. 2015;34:3133–43.10.1002/sim.6537Search in Google Scholar PubMed

[37] Mansournia MA, Etminan M, Danaei G, Kaufman JS, Collins G. Handling time varying confounding in observational research. Br Med J. 2017;359:j4587.10.1136/bmj.j4587Search in Google Scholar PubMed

[38] Hernán MA. The hazards of hazard ratios. Epidemiology. 2010;21:13–15.10.1097/EDE.0b013e3181c1ea43Search in Google Scholar PubMed PubMed Central

[39] Greenland S. Absence of confounding does not correspond to collapsibility of the rate ratio or rate difference. Epidemiology. 1996;7:498–501.10.1097/00001648-199609000-00008Search in Google Scholar

Received: 2017-12-16
Revised: 2018-12-20
Accepted: 2018-12-21
Published Online: 2019-01-15

© 2019 Walter de Gruyter GmbH, Berlin/Boston

Downloaded on 10.6.2023 from
Scroll to top button