UKB is a population-based cohort including more than 500 000 United Kingdom (UK) residents [23], aged 37–72 years at the recruitment in 2006–2010. Baseline assessments occurred at enrollment across 22 UK centers and included questionnaires and interviews on medical history, lifestyle and previous exposures as well as physical health measures. Diagnoses were self-reported during initial assessment and collected (retrospectively and prospectively) from hospital records, cause of death registers and primary care data [24]. The UKB study has received general ethical approval from the National Research Ethics Committee (REC reference 11/NW/0382). The present use of UKB data was specifically approved under application number 41,143, with data released on February 4, 2021. In addition, the analysis conducted in the current investigation was approved by the Swedish Ethical Review Authority (dnr: 2020–04415).

Osteoporosis diagnoses were retrieved from a curated data set (UKB field id: 131962 and 131965) comprising the date of the first osteoporosis diagnosis for each participant. These diagnoses were based on self-reported medical history (collected through a verbal-interview with a trained nurse at the visit to the UKB assessment center) in combination with data from national health registries, including primary care, hospital admissions and cause-of death records until 2020, mapped to the ICD-10 codes: M80 (osteoporosis with pathological fracture) and M81 (osteoporosis without pathological fracture). The BMD measures were not used to define osteoporosis diagnoses in the cohort. BMD, based on quantitative ultrasound (QUS) measurement of the calcaneus (heel bone), was determined at recruitment for the majority of the UKB participants (57%). BMD is reported as a T-score, which compares the participant’s bone density with that normally expected in someone of the same sex. BMD is therefore a continuous trait, expressed in standard deviations (SD).

Information on OCP use, including age at initiating and discontinuing exposure, was obtained from touchscreen questionnaires during the initial assessment. The relevant questions include: “Have you ever taken the contraceptive pill?“, “About how old were you when you first went on the contraceptive pill?”, and “How old were you when you last used the contraceptive pill?”. For participants who reported that they were still using OCP at the recruitment, the age of last use was set to the age at the recruitment. The difference between the age of last use and the age of first use was used as a proxy for duration of OCP use. For women who were current users at the time of assessment (when BMD measures were taken), the difference between the age at assessment and the age of first use was used as a proxy for the duration of OCP use. Covariates for the analysis include age, year of birth, menopausal status (premenopausal, postmenopausal, uncertain), reproductive surgeries (hysterectomy yes/no, bilateral oophorectomy yes/no), number of life births, BMI (as a continuous variable) and current smoking (yes/no). Covariate information used in the present analysis was obtained at the time of recruitment to the UKB, which is also when BMD measurements were conducted.

The UKB includes data from 273,375 women, and in this study, we included 257,185 women who self-identified as white Irish, white British, or other white backgrounds. Women with missing information on any of the covariates were excluded from all analyses. For osteoporosis, we used time-dependent Cox regression analyses, which require information on the age at initiation and discontinuation of OCP. In these analyses, participants were excluded if they lacked information on the age when they started using OCPs and/or the age of last use, resulting in 184,291 participants in the adjusted time-dependent Cox regression analyses.

For the BMD analyses, the cohort was limited to women who had undergone quantitative ultrasound (QUS) of the calcaneus (heel bone), resulting in 146,012 women (Supplementary Table S1). Women were first classified as either never users or ever users. Women who did not provide information on OCP use were excluded, resulting in 143 283 participants in the multivariate linear regression analyses. For the duration of use analyses, women who lacked information to estimate the duration of use were also excluded, resulting in 130 301 in the duration-stratified multivariate linear regression analyses.

All statistical analyses was carried out in R version 4.1.1. Chi-square tests (frequency) and t-tests (mean) were used to evaluate overall differences between subgroups. Participants lacking complete information on outcome, exposure and covariates were excluded from the analysis.

When analysing osteoporosis, a prospective study design was employed, and data were analysed using time-dependent Cox regression. This approach was used regardless of whether the outcome data were collected retrospectively (e.g., self-reported at the time of recruitment to the UKB) or from register-based diagnoses including events after the recruitment. We compiled a data frame that included the age at which changes in exposures occurred (i.e., when women started and stopped using OCPs) and age at which changes time-varying covariates occurred. Additionally, the data frame included whether participants had received an osteoporosis diagnosis by the end of the follow-up, and if so, the age at which they were first diagnosed with osteoporosis. To assess the effect of OCP use on osteoporosis diagnosis, adjusted hazard ratios (HRs) were estimated using time-dependent Cox proportional hazards models. Age was used as the primary time scale [25], as this has been suggested to be the optimal way to account for age-related confounding, especially for outcomes such as osteoporosis, which increase sharply with age.

Since OCP use varies over time, it was modelled as a time-varying exposure using Cox models for counting processes [26]. In this framework, women were classified as non-users before initiating OCP use, current users during active use, and previous users after discontinuation. Non-users served as the reference group in the analyses. A time-varying representation of OCP use avoids misclassification of users’ survival time before OCP initiation and defines a natural start of follow-up for never users [26, 27]. This approach accurately represents the OCP status and classifies the “event-free” person-time of OCP users before their age at initiation as the unexposed follow-up time [28]. Note that age is naturally adjusted for non-parametrically, as age is used as the primary timescale in the Cox regression modelling. Results are presented as relative change in hazard rate ratios between current OCP users and never users, or previous users and never users. In the Cox regression, smoking and menopause were also modelled as time-varying covariates, while other covariates were treated as time-fixed.

We used the coxph () function from the R survival package, applying a counting process approach to appropriately handle time-varying exposures and covariates. This method allows individuals to enter follow-up at any time point, with exposure status updated accordingly. Women were followed until the first diagnosis of osteoporosis or the end of study follow-up in 2020, whichever occurred first.

The analyses of BMD were performed using a retrospective study design and multiple linear regression with OCP first analyzed as a binary exposure (ever use vs. never use). Further analysis stratified the duration of OCP use into different duration intervals: never, \(\:\le\:\)1 year, 2–5 years, 6–10 years, 11–15 years and > 15 years. Duration groups were compared with never users as the reference using linear regression. Sensitivity analyses were performed in premenopausal, postmenopausal and women of uncertain menopausal status separately. To evaluate the effect of initiating OCP use in adolescence, we stratified women based on age at OCP initiation, those who started before age 18 and those who started at age 18 or older and compared both groups to never-users. Results were presented as differences in BMD in units of T-score in standard deviation units.