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

Tomohiro Shinozaki; Mohammad Ali Mansournia

doi:10.1515/ijb-2017-0103

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 and Mohammad Ali Mansournia

From the journal The International Journal of Biostatistics

https://doi.org/10.1515/ijb-2017-0103

Showing a limited preview of this publication:

Abstract

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.

Keywords: censoring; cohort study; Cox regression; hazard ratio; matching

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).

Appendix

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 α = (α₁, …, α_n)) partial likelihood is

∫∏i=1neαieβ∑k∈R(Yi∗)eαk(1+eβ)di1∗eαi∑k∈R(Yi∗)eαk(1+eβ)di0∗f(α;θ)dα=∏i=1neβ1+eβdi1∗11+eβdi0∗∫∏i=1neαi∑k∈R(Yi∗)eαkdi1∗+di0∗f(α;θ)dα,

where R(Y_i^*) is a risk set (of pairs) at time Y_i^* 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

∫∏i=1n∏j=0,1Yi∗eαi+βZijdij∗exp−Yij∗eαi+βZijf(α;θ)dα=∏i=1nYi∗di0∗+di1∗eβdi1∗∫∏i=1nexp{αi(di0∗+di1∗)−Yi∗eαi(1+eβ)+logf(α;θ)}dα=∏i=1nYi∗di0∗+di1∗eβdi1∗∫exp{g(α,β,θ)}dα,

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

∏i=1nYij∗di0∗+di1∗eβdi1∗∫∏i=1ndi0∗+di1∗+QiYi∗(1+eβ)di0∗+di1∗exp{−(di0∗+di1∗+Qi)+logf(α;θ)}dα=∏i=1neβ1+eβdi1∗11+eβdi0∗∫∏i=1n(di0∗+di1∗+Qi)di0∗+di1∗exp{−di0∗−di1∗−Qi+logf(α;θ)}dα.

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 X_i, partial likelihood is

∏i=1n∏j=0,1exp(βZij+γTXi)∑k∈R(Yi∗)(1+eβ)exp(γTXk)dij∗,

where R(Y_i^*) is a risk set (of pairs) at time Y_i^*. Note that d_i₁^*d_i₀^* = 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)/∑k∈R(Yi∗)exp(γTXi)di1∗+di0∗; from the first two factors of the likelihood including β, score equation of β is ∑_i{d_i₁^* – d_i₀^*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 X_i, the score equations are

∑ijdij∗−Yi∗exp(α+βZij+γTXi)=0,∑ijZijdij∗−Yi∗exp(α+βZij+γTXi)=0,∑ijXidij∗−Yi∗exp(α+βZij+γTXi)=0.

The first two equations jointly imply ∑_i{d_i₀^* – Y_i^*exp(α + γ^TX_i)} = 0 and ∑_i {d_i₁^* – Y_i^*exp(α + γ^TX_i)e^β} = 0; hence a maximum likelihood estimate of β is the solution of ∑_i {d_i₁^* – d_i₀^*e^β} = 0, which is independent of α and γ.

β is estimated independent of γ in both equations; therefore, robust variance of β is estimated as {∑_i m_i(β)/∂β}^–1∑_i m_i(β)²{∑_i ∂m_i(β)/∂β}^–1, where m_i(β) = d_i₁^* – d_i₀^*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 X_i 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 X_i as a scalar X_i (though the following result holds if X_i is a vector). As shown in the Appendix B, the score equations for the Poisson model conditional on shared covariates X_i are

Sα=∑ijdij∗−Yi∗exp(α+βZij+γXi)=0,Sβ=∑ijZijdij∗−Yi∗exp(α+βZij+γXi)=0,Sγ=∑ijXidij∗−Yi∗exp(α+βZij+γXi)=0.

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

(13)eβ∑idi0∗eβ∑idi0∗eβ1+eβ∑iXi(di0∗+di1∗)eβ∑idi0∗(1+eβ)∑idi0∗∑iXi(di0∗+di1∗)eβ1+eβ∑iXi(di0∗+di1∗)∑iXi(di0∗+di1∗)(1+eβ)∑iXi2Yi∗(di0∗+di1∗).

We partition matrix (13) as

A=eβ∑idi0∗B=eβ∑idi0∗,eβ1+eβ∑iXi(di0∗+di1∗)C=eβ∑idi0∗,eβ1+eβ∑iXi(di0∗+di1∗)TD=(1+eβ)∑idi0∗∑iXi(di0∗+di1∗)∑iXi(di0∗+di1∗)(1+eβ)∑iXi2Yi∗(di0∗+di1∗)

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β)2∑idi0∗∑iXi2Yi∗(di0∗+di1∗)−∑iXi(di0∗+di1∗)2]−1,

BD−1C=Keβ(1+eβ)∑idi0∗∑iXi2Yi∗(di0∗+di1∗)−11+eβ∑iXi(di0∗+di1∗)∑i(di0∗+di1∗)−∑idi0∗∑i(di0∗+di1∗)+∑iXi(di0∗+di1∗)∑idi0∗TC=Ke2β∑idi0∗1+eβ[(1+eβ)2∑idi0∗∑iXi2Yi∗(di0∗+di1∗)−∑i(di0∗+di1∗)∑iXi(di0∗+di1∗)−∑iXi(di0∗+di1∗)∑iXi(di0∗+di1∗)+∑iXi(di0∗+di1∗)2]=Ke2β∑idi0∗1+eβK−1=e2β∑idi0∗1+eβ.

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

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

Fixed-effect Cox models’ partial likelihood is

∏i=1n∏j=0,1eαiexp(βZij+γTXij)∑k∈R(Yi∗)eαkexp(β+γTXk1)+exp(γTXk0)dij∗.

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: ∑k∈R(Yi∗),k≠ieαkexp(β+γTXk1)+exp(γTXk0)=0. Thus, the partial likelihood reduces to

∏i=1n∏j=0,1eαiexp(βZij+γTXij)eαiexp(β+γTXi1)+exp(γTXi0)dij∗,

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

∏i=1nexp(β+γTXi1)exp(β+γTXi1)+exp(γTXi0)di1∗exp(γTXi0)exp(β+γTXi1)+exp(γTXi0)di0∗.

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

∑jin pairidij∗−Yi∗exp(αi+βZij+γTXij)=0fori=1,…,n.

Substituting Y_i^*exp(α_i) = (d_i₀^* + d_i₁^*)/{exp(β + γ^TX_i₁) + exp(γ^TX_i₀)} from the above conditions, ∑_ijZ_ij{d_ij^* – Y_i^*exp(α_i + βZ_ij + γ^TX_ij)} and ∑_ijX_ij{d_ij^* – Y_i^*exp(α_i + βZ_ij + γ^TX_ij)} 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 pairidij∗−Yi∗exp(αi+βZij+γTXij)=0fori=1,…,n,Sβ=∑ijZijdij∗−Yi∗exp(αi+βZij+γTXij)=0,Sγ=∑ijXijdij∗−Yi∗exp(αi+βZij+γTXij)=0.

The first equation of Sαi=0 imply Y_i^*exp(α_i) = (d_i₀^* + d_i₁^*)/{exp(β + γ^TX_i₁) + exp(γ^TX_i₀)}.

Deleting Y_i^*exp(α_i) from the second equation, S_β = ∑_i [d_i₁^* – (d_i₀^* + d_i₁^*)exp(β + γ^TX_i₁)/{exp(β + γ^TX_i₁) + exp(γ^TX_i₀)}] = ∑_i [{d_i₁^*exp(γ^TX_i₀) + d_i₀^* exp(β + γ^TX_i₁)}/{exp(β + γ^TX_i₁) + exp(γ^TX_i₀)}]. This is the same as the first derivative of the log-partial likelihood of the stratified Cox model (1)

∑idi1∗[β+γTXi1−logexp(β+γTXi1)+exp(γTXi0)]+∑idi0∗[γTXi0−logexp(β+γTXi1)+exp(γTXi0)]

with respect to β.

Similarly, deleting Y_i^*exp(α_i) from the third equation, S_γ = ∑_i [X_i₁{d_i₁^*exp(β + γ^TX_i₁) + d_i₁^*exp(γ^TX_i₀) – d_i₁^*exp(β + γ^TX_i₁) – d_i₀^*exp(β + γ^TX_i₁)}/{exp(β + γ^TX_i₁) + exp(γ^TX_i₀)}] + ∑_i [X_i₀{d_i₀^*exp(β + γ^TX_i₁) + d_i₀^*exp(γ^TX_i₀) – d_i₁^*exp(γ^TX_i₀) – d_i₀^*exp(γ^TX_i₀)}/{exp(β + γ^TX_i₁) + exp(γ^TX_i₀)}] = ∑_i [{X_i₁d_i₁^*exp(γ^TX_i₀) – X_i₁d_i₀^*exp(β + γ^TX_i₁) + X_i₀d_i₀^*exp(β + γ^TX_i₁) – X_i₀d_i₁^*exp(γ^TX_i₀)}/{exp(β + γ^TX_i₁) + exp(γ^TX_i₀)}]. As before, the first derivative of the log-partial likelihood of the stratified Cox model (1) with respect to γ provides the same function.

Acknowledgements

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.

References

[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 https://j.mp/1sexgVw. 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

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

Abstract

Appendix

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

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

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

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

Acknowledgements

References

Journal and Issue

Articles in the same Issue