Advertisement
Original Article|Articles in Press

Cox regression using a calendar time scale was unbiased in simulations of COVID-19 vaccine effectiveness & safety

Open AccessPublished:February 15, 2023DOI:https://doi.org/10.1016/j.jclinepi.2023.02.012

      Abstract

      Background

      Observational studies on corona virus disease 2019 (COVID-19) vaccines compare event rates in vaccinated and unvaccinated person time using Poisson or Cox regression. In Cox regression, the chosen time scale needs to account for the time-varying incidence of severe acute respiratory syndrome corona virus 2 (SARS-CoV-2) infection and COVID-19 vaccination.We aimed to quantify bias in person-time based methods, with and without adjustment for calendar time, using simulations and empirical data analysis.

      Methods

      We simulated 500,000 individuals who were followed for 365 days, and a point exposure resembling COVID-19 vaccination (cumulative incidence 80%). We generated an effectiveness outcome, emulating the incidence of severe acute respiratory syndrome corona virus 2 infection in Denmark during 2021 (risk 10%), and a safety outcome with seasonal variation (myocarditis, risk 1/5,000). Incidence rate ratios (IRRs) were set to 0.1 for effectiveness and 5.0 for safety outcomes. IRRs and hazard ratios (HRs) were estimated using Poisson and Cox regression with a time under observation scale, and a calendar time scale. Bias was defined as estimated IRR or HR−true IRR. Further, we obtained estimates for both outcomes using data from the Danish health registries.

      Results

      Unadjusted IRRs (biaseffectivenes +0.16; biassafety −2.09) and HRs estimated using a time-under-observation scale (+0.28;-2.15) were biased. Adjustment for calendar time reduced bias in Cox (+0.03; +0.33) and Poisson regression (0.00; −0.28). Cox regression using a calendar time scale was least biased (0.00, +0.12). When analyzing empirical data, adjusted Poisson and Cox regression using a calendar time scale yielded estimates in accordance with existing evidence.

      Conclusion

      Lack of adjustment for the time-varying incidence of COVID-19 related outcomes may severely bias estimates.

      Keywords

      What is new?

        Key Findings:

      • Cox regression using a calendar time scale and Poisson regression adjusted for calendar time using restricted cubic splines were unbiased in simulations of COVID-19 vaccine effectiveness and safety.

        What this adds to what was known?

      • The highly time-varying incidence of SARS-CoV-2 infection and the rapid uptake of COVID-19 vaccines may bias analyses regarding COVID-19 vaccine effectiveness
      • When using calendar time as the underlying time scale in Cox regression, the researcher does not need to make any parametric assumptions regarding the instantaneous risk for the outcome of interest.

        What is the implication and what should change now?

      • Observational analyses regarding COVID-19 vaccine effectiveness or safety should adjust for calendar time, preferably using flexible methods such as restricted cubic splines or by using calendar time as the underlying time scale in Cox regression.
      • If Cox regression is used in such analyses, the underlying time scale should be stated explicitly.

      1. Introduction

      An increasing number of observational analyses regarding corona virus disease 2019 (COVID-19) vaccine effectiveness and safety have been published recently. Often, these studies compare the occurrence of an event of interest, e.g., documented severe acute respiratory syndrome corona virus 2 (SARS-CoV-2) infection [
      • Haas E.J.
      • Angulo F.J.
      • McLaughlin J.M.
      • Anis E.
      • Singer S.R.
      • Khan F.
      • et al.
      Impact and effectiveness of mRNA BNT162b2 vaccine against SARS-CoV-2 infections and COVID-19 cases, hospitalisations, and deaths following a nationwide vaccination campaign in Israel: an observational study using national surveillance data.
      ], hospitalization for COVID-19 or adverse events related to vaccination [
      • Husby A.
      • Hansen J.V.
      • Fosbøl E.
      • Thiesson E.M.
      • Madsen M.
      • Thomsen R.W.
      • et al.
      SARS-CoV-2 vaccination and myocarditis or myopericarditis: population based cohort study.
      ], relative to vaccinated and unvaccinated person time. Due to the rapid uptake of the COVID-19 vaccines in countries with up-to-date registries such as Israel [
      • Rosen B.
      • Waitzberg R.
      • Israeli A.
      Israel’s rapid rollout of vaccinations for COVID-19.
      ], the United Kingdom [] and Scandinavia [
      • Reilev M.
      • Olesen M.
      • Kildegaard H.
      • Stovring H.
      • Andersen J.H.
      • Hallas J.
      • et al.
      Changing characteristics over time of individuals receiving COVID-19 vaccines in Denmark: A population-based descriptive study of vaccine uptake.
      ], much of the unvaccinated and vaccinated person time has not been accrued simultaneously. This increases the risk of temporal biases affecting study results [
      • Griffin B.A.
      • Anderson G.L.
      • Shih R.A.
      • Whitsel E.A.
      Use of alternative time scales in Cox proportional hazard models: implications for time-varying environmental exposures.
      ], especially with the highly time-varying incidence of SARS-CoV-2 infection or when analyzing safety outcomes with seasonal variation. The estimate of interest in such studies is the incidence rate ratio (IRR) or hazard ratio (HR) associating vaccination to risk of SARS-CoV-2 infection, or a safety outcome, and is usually estimated using a Poisson or Cox proportional hazards regression model [
      • Rothman K.J.
      • Lash T.L.
      • VanderWeele T.J.
      • Haneuse S.
      Modern epidemiology.
      ]. Adjustment for temporal differences between vaccinated and unvaccinated person time can be performed through covariate adjustment or by choosing an appropriate time scale in Cox regression analyses. We therefore aimed to compare different approaches to adjust for calendar time using Poisson and Cox regression in analyses of vaccine effectiveness and safety in the presence of strong temporal trends related to the COVID-19 pandemic waves.

      2. Methods

      We performed a simulation study, investigating the performance of Cox and Poisson regression with and without covariate adjustment for calendar time, by matching on calendar time, and with a time under observation and calendar time scale for Cox regression analyses. We tested whether these estimators yielded unbiased results by simulating the observed SARS-CoV-2 incidence and COVID-19 vaccine uptake in Denmark during the year of 2021. We further evaluated the estimators using empirical data, by investigating COVID-19 vaccine effectiveness toward the delta variant of SARS-CoV-2 and the risk of myocarditis following vaccination with mRNA-1273.

      2.1 Simulation

      We simulated data representing key aspects of the COVID-19 pandemic in Denmark during the year of 2021. We simulated 500,000 individuals who were observed from 01 January 2021 (day 1) to 31 December 2021 (day 365).

      2.2 Exposure

      To emulate key aspects of the COVID-19 immunization program (for details, see appendix 1 in the supplementary material), we simulated a point exposure mimicking administration of the first vaccine dose, with a cumulative incidence proportion of 80% during the year of observation. Among vaccinated, vaccination dates were drawn from a normal distribution with a mean of 180 days (30 June 2021), and a standard deviation of 30 days. Individuals were considered exposed from the first day of vaccination and for the rest of the year.

      2.3 Outcomes

      We simulated two different types of outcomes:
      • 1)
        An effectiveness outcome with a high cumulative incidence proportion at the end of the year (10%) resembling the outcome of documented SARS-CoV-2 infection. For this, we obtained daily case numbers of testing positive for SARS-CoV-2 in Denmark during 2021 [
        • Ritchie H.
        • Mathieu E.
        • Rodés-Guirao L.
        • Appel C.
        • Giattino C.
        • Ortiz-Ospina E.
        • et al.
        Coronavirus pandemic (COVID-19).
        ], divided these by the total size of the adult population (4.7 million) and adjusted the daily risk to sum up to a cumulative incidence proportion of 0.1 on day 365.
      • 2)
        For the safety outcome, we simulated a rare event that exhibits seasonal variation with a winter peak. Examples of such outcomes are myocarditis or Guillain-Barré syndrome which have previously been associated with COVID-19 immunization [
        • Husby A.
        • Hansen J.V.
        • Fosbøl E.
        • Thiesson E.M.
        • Madsen M.
        • Thomsen R.W.
        • et al.
        SARS-CoV-2 vaccination and myocarditis or myopericarditis: population based cohort study.
        ,
        • Patone M.
        • Handunnetthi L.
        • Saatci D.
        • Pan J.
        • Katikireddi S.V.
        • Razvi S.
        • et al.
        Neurological complications after first dose of COVID-19 vaccines and SARS-CoV-2 infection.
        ,
        • Patone M.
        • Mei X.W.
        • Handunnetthi L.
        • Dixon S.
        • Zaccardi F.
        • Shankar-Hari M.
        • et al.
        Risks of myocarditis, pericarditis, and cardiac arrhythmias associated with COVID-19 vaccination or SARS-CoV-2 infection.
        ,
        • Karlstad Ø.
        • Hovi P.
        • Husby A.
        • Härkänen T.
        • Selmer R.M.
        • Pihlström N.
        • et al.
        SARS-CoV-2 vaccination and myocarditis in a nordic cohort study of 23 million residents.
        ,
        • Li X.
        • Raventós B.
        • Roel E.
        • Pistillo A.
        • Martinez-Hernandez E.
        • Delmestri A.
        • et al.
        Association between covid-19 vaccination, SARS-CoV-2 infection, and risk of immune mediated neurological events: population based cohort and self-controlled case series analysis.
        ]. In detail, we obtained the daily risk of the outcome based on a normal distribution with a mean of 30 days and a standard deviation of 50 days, which was added to a constant baseline rate. To convert the obtained densities to a daily risk, the density for each day was standardized to sum up to one and then multiplied by the desired cumulative incidence proportion of 1/5,000. As the described normal distribution provided densities for days before day 0 (January 1st) but practically no densities for the end of the year (day 330-365), we considered the densities for days before day 0 to occur during the end of the year, i.e., the density for day −1 was set to day 364, the density for day −2 to 363, and so on. This creates a wrap around of the daily hazard. For a graphical depiction of the daily hazard and cumulative incidence of the simulated outcomes and exposure between day 1 and 365, see Figure 1.
        Figure thumbnail gr1
        Fig. 1Daily simulated hazard and simulated cumulative incidence of the effectiveness outcome, safety outcome and vaccination for the study period among unvaccinated individuals.
        Figure thumbnail gr2
        Fig. 2Mean incidence rate ratios (IRRs) and hazard ratios obtained by analyzing simulated data. Bold horizontal lines represent the point estimate ±1 empirical standard error (SE), thin lines the point estimate ±2 SE. The dashed line represents the true, simulated incidence rate ratio. ∗ Using a time under observation scale.
      Event times for each outcome were obtained by performing a Bernoulli trial for each individual and day based on the above-described probabilities.

      2.4 Follow-up

      All individuals were followed from day 1 and until day 365 or the onset of the analyzed outcome of interest. Without loss of generalizability, but to maintain simplicity, we ignored any other sources of censoring such as migration. Exposure status was considered a time-varying covariate, i.e., all individuals contributed unexposed person time in the beginning of the study and individuals contributed exposed person-time from the date of vaccination and onwards.
      We considered two different time scales in this analysis, a ‘time under observation’ scale and a calendar time scale. For both time scales, unvaccinated individuals entered the study on January 1st which was defined as day 0. When using the ‘time under observation’ scale, the time scale was reset upon a change in exposure status, i.e, vaccination, and the date of vaccination was counted as day 0 for the vaccinated person time accrued by an individual. When using the calendar time scale, the time scale was not reset for vaccinated individuals.

      2.5 Scenarios

      For each type of outcome we tested a range of true IRRs when comparing exposed to unexposed person time. For the effectiveness outcome, we tested true IRRs of 1.0, 0.5, and 0.1 corresponding, respectively, to no protection against SARS-CoV-2 infection, the presumed vaccine effectiveness against the Omicron variant of SARS-CoV-2 shortly after immunization [
      • Hansen C.H.
      • Schelde A.B.
      • Moustsen-Helm I.R.
      • Emborg H.D.
      • Krause T.G.
      • Mølbak K.
      • et al.
      Vaccine effectiveness against SARS-CoV-2 infection with the Omicron or Delta variants following a two-dose or booster BNT162b2 or mRNA-1273 vaccination series: A Danish cohort study.
      ], and vaccine effectiveness against the Delta variant [
      • Tartof S.Y.
      • Slezak J.M.
      • Fischer H.
      • Hong V.
      • Ackerson B.K.
      • Ranasinghe O.N.
      • et al.
      Effectiveness of mRNA BNT162b2 COVID-19 vaccine up to 6 months in a large integrated health system in the USA: a retrospective cohort study.
      ]. For the safety outcome we simulated true IRRs of 1.0, 2.0 and 5.0, i.e., no increased risk, and different risk increases corresponding to the event of myocarditis after COVID-19 immunization with the mRNA-1273 vaccine [
      • Husby A.
      • Hansen J.V.
      • Fosbøl E.
      • Thiesson E.M.
      • Madsen M.
      • Thomsen R.W.
      • et al.
      SARS-CoV-2 vaccination and myocarditis or myopericarditis: population based cohort study.
      ,
      • Patone M.
      • Mei X.W.
      • Handunnetthi L.
      • Dixon S.
      • Zaccardi F.
      • Shankar-Hari M.
      • et al.
      Risks of myocarditis, pericarditis, and cardiac arrhythmias associated with COVID-19 vaccination or SARS-CoV-2 infection.
      ,
      • Karlstad Ø.
      • Hovi P.
      • Husby A.
      • Härkänen T.
      • Selmer R.M.
      • Pihlström N.
      • et al.
      SARS-CoV-2 vaccination and myocarditis in a nordic cohort study of 23 million residents.
      ]. All simulated effects were constant over time.

      2.6 Estimators

      For each scenario, we estimated IRRs and HRs, with and without adjusting for calendar time using restricted cubic splines. In adjusted Poisson and Cox regression, we split the observed person time every 14 days based on calendar time and constructed restricted cubic splines with four knots placed at the fifth, 35th, 65th, and 95th quantile of the distribution of observed start dates for each segment of person time [
      • Harrell FE Jr.,
      Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis.
      ]. The resulting spline-based variables were included as independent variables in the Poisson model and Cox model using a time under observation scale. IRRs comparing exposed to unexposed person time were estimated using Poisson regression with the logarithm of the accrued person time being used as the offset.
      We estimated HRs using the semi-parametric Cox proportional hazards regression model, which allows the baseline hazard to vary freely over time. This feature inherently adjusts for any temporal trends in the underlying time scale, which typically can be age, calendar time or time since coming under observation. We therefore hypothesized that choice of the underlying time scale is crucial for obtaining unbiased results in Cox regression. We estimated HRs using time under observation as the underlying time scale, i.e., days under the current exposure status, and a calendar time scale (denoted CoxCT). For the time under observation scale, we also evaluated a simple 1:1 matching estimator: For each vaccinated individual we randomly sampled an unvaccinated individual, and the vaccination date of the vaccinated individual was assigned as the date of coming under observation for the unvaccinated individual. Only individuals who remained unvaccinated and had not experienced the outcome of interest were eligible for matching. Vaccinated and unvaccinated individuals who could not be matched were excluded.
      Time under observation is a common time scale in pharmacoepidemiological studies [
      • Lund J.L.
      • Richardson D.B.
      • Stürmer T.
      The active comparator, new user study design in pharmacoepidemiology: historical foundations and contemporary application.
      ,
      • Suissa S.
      Immortal time bias in pharmaco-epidemiology.
      ], mostly used when comparing two treatments, but the same logic can be applied to an anchoring by cohort entry, i.e., time 0 is placed at the beginning of observation under a given exposure status (”cohort entry”). We expected the analysis using this time scale to be susceptible to temporal biases and HRs estimated using a calendar time scale to be unaffected. We evaluated that all estimators yielded unbiased results in the absence of temporal trends by analyzing an outcome with a constant rate simulated using an exponential decay function. For a graphical depiction of how event times are compared under different time scales, please see supplementary figure 1.
      All simulations were repeated 1,000 times and from these we calculated the mean estimate on the log-scale, standard deviation of all estimates on the log-scale, the root mean squared error, and the coverage probability of the estimated 95% confidence intervals (CIs) for all described estimators. Bias was both calculated as the absolute difference between the estimated HR or IRR and the true IRR, i.e., exp(βexposure) – IRR, and also on the log-scale, i.e., βexposure–log(IRR). We further calculated Monte Carlo standard errors and 95% Monte Carlo uncertainty intervals for all performance metrics [
      • Morris T.P.
      • White I.R.
      • Crowther M.J.
      Using simulation studies to evaluate statistical methods.
      ].

      2.7 Empirical data analysis

      To evaluate whether results obtained through analysis of simulated data are applicable to the analysis of real-world data, we performed a similar analysis using information from the Danish civil registration system [
      • Pedersen C.B.
      The Danish civil registration system.
      ], COVID-19 test results [
      • Voldstedlund M.
      • Haarh M.
      • Mølbak K.
      Representatives the MB of. The Danish microbiology database (MiBa) 2010 to 2013.
      ] and information on COVID-19 immunization from the Danish vaccination register [
      • Grove Krause T.
      • Jakobsen S.
      • Haarh M.
      • Mølbak K.
      The Danish vaccination register.
      ]. The study period was December 27th 2020 to November 15th 2021.
      For the vaccine effectiveness analysis, we identified all adult individuals living in the Copenhagen municipality at the beginning of the study period and followed them until 15 November 2021, death, emigration or the outcome of interest. The exposure of interest was immunization with the BNT162b2 vaccine [
      • Polack F.P.
      • Thomas S.J.
      • Kitchin N.
      • Absalon J.
      • Gurtman A.
      • Lockhart S.
      • et al.
      Safety and Efficacy of the BNT162b2 mRNA covid-19 vaccine.
      ], and we considered individuals exposed from the day of administration of a second dose, i.e., individuals started contributing exposed person time on this date and during the following 60 days. For the empirical data analysis, we chose to analyze a short risk window following immunization as we expected the effect of vaccination to be constant during this period, and therefore comparable to the simulation. Individuals contributed unvaccinated person time until occurrence of the above-mentioned events or administration of a first dose of a COVID-19 vaccine. Prior evidence indicates that immunization with BNT162b2 reduced the risk of documented SARS-CoV-2 infection with the Delta variant with 85%–95% [
      • Haas E.J.
      • Angulo F.J.
      • McLaughlin J.M.
      • Anis E.
      • Singer S.R.
      • Khan F.
      • et al.
      Impact and effectiveness of mRNA BNT162b2 vaccine against SARS-CoV-2 infections and COVID-19 cases, hospitalisations, and deaths following a nationwide vaccination campaign in Israel: an observational study using national surveillance data.
      ,
      • Tartof S.Y.
      • Slezak J.M.
      • Fischer H.
      • Hong V.
      • Ackerson B.K.
      • Ranasinghe O.N.
      • et al.
      Effectiveness of mRNA BNT162b2 COVID-19 vaccine up to 6 months in a large integrated health system in the USA: a retrospective cohort study.
      ].
      For the safety outcome, we followed all adult individuals living in Denmark from the beginning of the study period. The exposure of interest was administration of a first or second dose of the mRNA-1273 vaccine [
      • Baden L.R.
      • El Sahly H.M.
      • Essink B.
      • Kotloff K.
      • Frey S.
      • Novak R.
      • et al.
      Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine.
      ], and individuals contributed exposed person time for 28 days following administration of a vaccine dose. Administration of any other COVID-19 vaccine or documented SARS-CoV-2 infection was considered a reason for censoring. The mRNA-1273 vaccine has previously been associated with an increased risk of myocarditis, with estimates of relative risk greater than 2.5 being reported in multiple studies [
      • Husby A.
      • Hansen J.V.
      • Fosbøl E.
      • Thiesson E.M.
      • Madsen M.
      • Thomsen R.W.
      • et al.
      SARS-CoV-2 vaccination and myocarditis or myopericarditis: population based cohort study.
      ,
      • Patone M.
      • Mei X.W.
      • Handunnetthi L.
      • Dixon S.
      • Zaccardi F.
      • Shankar-Hari M.
      • et al.
      Risks of myocarditis, pericarditis, and cardiac arrhythmias associated with COVID-19 vaccination or SARS-CoV-2 infection.
      ,
      • Karlstad Ø.
      • Hovi P.
      • Husby A.
      • Härkänen T.
      • Selmer R.M.
      • Pihlström N.
      • et al.
      SARS-CoV-2 vaccination and myocarditis in a nordic cohort study of 23 million residents.
      ]. Myocarditis and pericarditis are currently the only serious adverse events of the BNT162b2 and mRNA-1273 vaccine that have been reported in multiple, large observational studies. The outcome of myocarditis was defined as a first ever hospital diagnosis with an international classification of disease 10th revision (ICD-10) code of I40.0, I40.1, I40.9, I41.1, I41.8, or I51.4.
      For the effectiveness and safety outcomes, we obtained the point estimate and a 95% CI for all described estimators. Estimates were not adjusted for factors other than calendar time, as we did not aim to estimate causal effects but explore the differences between estimators.
      All statistical analyses were performed using Stata version 17 (StataCorp LLC). All source code used for the simulations and empirical analyses are available from https://gitlab.sdu.dk/lclund/vaccine-simulation/.

      3. Results

      3.1 Simulation

      In simulated analyses of vaccine effectiveness with a true IRR of 0.1, we obtained a mean point estimate of 0.26 (bias = estimated IRR–true IRR = 0.16; log bias = log(IRR) – log(true IRR) = 0.94) using unadjusted Poisson regression and 0.10 when adjusting for calendar time. In Cox regression analyses, the mean HR was 0.38 (bias 0.28; log bias 1.33) using a time under observation scale, which was reduced to 0.13 (bias 0.03; log bias 0.23) when adjusting for calendar time. Using calendar time as the underlying scale or matching on calendar time yielded an unbiased estimate (HR 0.10) (Figure 2). Adjusted Poisson, CoxCT and matching estimators provided valid 95% CIs (coverage probability of 95%). When evaluating a true IRR of 0.5, adjusted Poisson, CoxCT and matching estimators again yielded unbiased results (IRR 0.50 and HR 0.50). The other estimators were substantially biased (bias of 0.12, 0.76, and 1.44; log bias 0.93, 1.35, and 0.21). In simulated analyses of vaccine safety, with a true IRR of 5.0, no estimators yielded unbiased results. We observed the lowest bias for the CoxCT estimator (HR 5.12, bias 0.12, log bias 0.02), the adjusted Cox estimator (HR 5.33, bias 0.33; log bias 0.06), the adjusted Poisson estimator (IRR 4.72, bias −0.28; log bias −0.06), and the matching estimator (HR 5.36, bias 0.36, log bias 0.07), although the other estimators were strongly biased (log bias −0.54 and −0.56). Obtained estimates varied substantially with empirical standard errors between 0.15 (unadjusted Poisson) and 0.66 (adjusted Cox). When reducing the strength of association to a true IRR of 2.0, adjusted Poisson and CoxCT estimators were close to unbiased (0.01 and 0.06; log bias 0.01 and 0.03). Finally, we visualized the observed cumulative incidence for a single iteration of the simulated effectiveness and safety outcome using both a time under observation and calendar time scale (Figure 3). The bias that was observed when using the time under observation scale may be explained by the cumulative incidence of vaccinated individuals being shifted to the left when using the time under observation scale, leading to vaccinated individuals followed during the last half of the study period being compared to unvaccinated individuals who were followed during the early part of the year.
      Figure thumbnail gr3
      Fig. 3Observed cumulative incidence for simulated safety and effectiveness outcomes stratified by vaccination status when making use of different time scales. The depicted data show a single iteration of the simulation for the effectiveness (cumulative incidence proportion 0.1, true incidence rate ratio 0.1) and safety outcome (cumulative incidence proportion 1/5,000, true incidence rate ratio 5.0). The left column depicts the cumulative incidence for the time under observation scale (time 0 for unvaccinated individuals on January 1st, time 0 for vaccinated individuals on the date of vaccination), and the right column demonstrates a calendar time scale (time 0 on January 1st for all individuals, vaccinated individuals come under observation when vaccinated, e.g., on day 180). Depicted true IRRs are 0.1 for the effectiveness outcome and 5.0 for the safety outcome. Cumulative incidences are only depicted until day 270 for the time under observation scale, as no vaccinated individual was observed for more than 270 days.
      Figure thumbnail gr4
      Fig. 4Estimated IRRs and hazard ratios from the empirical data analysis. The effectiveness outcome was documented SARS-CoV-2 infection. The safety outcome was a first ever hospital diagnosis of myocarditis. The bold black lines represent the point estimate ±1 estimated standard error (∼68% confidence interval), the thin black lines represent the 95% confidence interval. ∗ Using a time under observation scale. Estimates for the adjusted and matched Cox regression could not be reported due to the regression not converging to meaningful results. SARS-CoV-2, severe acute respiratory syndrome corona virus 2.
      For detailed results of all scenarios, see Table 1.
      Table 1Results from simulations regarding vaccine effectiveness and safety
      EstimatorMeanBias
      Calculated as the estimated relative risk minus the simulated incidence rate ratio.
      (95% UI)
      Log bias
      Calculated as the estimated regression coefficient for the exposure minus the logarithm of the simulated incidence rate ratio.
      (MCSE)
      SE (MCSE)RMSE (95% UI)Coverage, % (MCSE)
      Outcome: Effectiveness (IRR=0.1)
       Poisson0.260.16 (0.15, 0.16)0.94 (0.00)0.02 (0.00)0.94 (0.94, 0.94)0 (0)
       Poisson, adjusted0.100.00 (−0.00, 0.00)0.00 (0.00)0.02 (0.00)0.02 (0.02, 0.02)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      0.380.28 (0.28, 0.28)1.33 (0.00)0.02 (0.00)1.33 (1.32, 1.33)0 (0)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , adjusted
      0.130.03 (0.03, 0.03)0.23 (0.00)0.12 (0.00)0.26 (0.26, 0.27)54 (2)
       Cox, calendar time0.10−0.00 (−0.00, 0.00)−0.00 (0.00)0.02 (0.00)0.02 (0.02, 0.02)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , matching
      0.10−0.00 (−0.00, 0.00)−0.00 (0.00)0.04 (0.00)0.04 (0.04, 0.04)95 (1)
      Outcome: Effectiveness (IRR=0.5)
       Poisson1.260.76 (0.76, 0.76)0.93 (0.00)0.01 (0.00)0.93 (0.93, 0.93)0 (0)
       Poisson, adjusted0.50−0.00 (−0.00, −0.00)−0.00 (0.00)0.01 (0.00)0.01 (0.01, 0.01)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      1.941.44 (1.43, 1.44)1.35 (0.00)0.01 (0.00)1.35 (1.35, 1.35)0 (0)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , adjusted
      0.620.12 (0.11, 0.12)0.21 (0.00)0.07 (0.00)0.22 (0.22, 0.22)10 (1)
       Cox, calendar time0.500.00 (−0.00, 0.00)0.00 (0.00)0.01 (0.00)0.01 (0.01, 0.01)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , matching
      0.500.00 (−0.00, 0.00)0.00 (0.00)0.02 (0.00)0.02 (0.02, 0.02)95 (1)
      Outcome: Effectiveness (IRR=1.0)
       Poisson2.501.50 (1.50, 1.50)0.91 (0.00)0.01 (0.00)0.92 (0.91, 0.92)0 (0)
       Poisson, adjusted0.99−0.01 (−0.01, −0.00)−0.01 (0.00)0.01 (0.00)0.01 (0.01, 0.01)93 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      3.942.94 (2.94, 2.94)1.37 (0.00)0.01 (0.00)1.37 (1.37, 1.37)0 (0)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , adjusted
      1.210.21 (0.20, 0.21)0.19 (0.00)0.05 (0.00)0.19 (0.19, 0.20)4 (1)
       Cox, calendar time1.000.00 (−0.00, 0.00)0.00 (0.00)0.01 (0.00)0.01 (0.01, 0.01)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , matching
      1.000.00 (−0.00, 0.00)0.00 (0.00)0.02 (0.00)0.02 (0.02, 0.02)96 (1)
      Outcome: Safety (IRR=1.0)
       Poisson0.58−0.42 (−0.43, −0.41)−0.55 (0.01)0.23 (0.01)0.59 (0.58, 0.61)31 (1)
       Poisson, adjusted1.050.05 (0.03, 0.08)0.05 (0.01)0.37 (0.01)0.38 (0.36, 0.39)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      0.55−0.45 (−0.46, −0.44)−0.59 (0.01)0.23 (0.01)0.64 (0.62, 0.65)24 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , adjusted
      1.020.02 (−0.06, 0.10)0.01 (0.04)1.27 (0.03)1.27 (1.21, 1.33)94 (1)
       Cox, calendar time1.040.04 (0.01, 0.06)0.04 (0.01)0.39 (0.01)0.39 (0.37, 0.41)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , matching
      0.95−0.05 (−0.15, 0.05)−0.06 (0.05)1.73 (0.04)1.73 (-
      No lower bound reported, due to a negative value being obtained on the mean squared error scale.
      , 2.58)
      97 (1)
      Outcome: Safety (IRR=2.0)
       Poisson1.16−0.84 (−0.85, −0.83)−0.54 (0.01)0.18 (0.00)0.57 (0.56, 0.58)14 (1)
       Poisson, adjusted2.010.01 (−0.03, 0.05)0.01 (0.01)0.31 (0.01)0.31 (0.30, 0.33)96 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      1.12−0.88 (−0.90, −0.87)−0.58 (0.01)0.18 (0.00)0.61 (0.60, 0.62)11 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , adjusted
      2.130.13 (0.01, 0.26)0.06 (0.03)0.95 (0.02)0.95 (0.90, 0.99)95 (1)
       Cox, calendar time2.060.06 (0.01, 0.10)0.03 (0.01)0.34 (0.01)0.34 (0.33, 0.36)96 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , matching
      2.100.10 (0.04, 0.17)0.05 (0.02)0.51 (0.01)0.51 (0.48, 0.54)98 (0)
      Outcome: Safety (IRR=5.0)
       Poisson2.91−2.09 (−2.12, −2.07)−0.54 (0.00)0.15 (0.00)0.56 (0.55, 0.57)6 (1)
       Poisson, adjusted4.72−0.28 (−0.36, −0.20)−0.06 (0.01)0.27 (0.01)0.27 (0.26, 0.29)95 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      2.85−2.15 (−2.18, −2.12)−0.56 (0.00)0.15 (0.00)0.58 (0.57, 0.59)5 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , adjusted
      5.330.33 (0.12, 0.55)0.06 (0.02)0.66 (0.01)0.66 (0.63, 0.69)95 (1)
       Cox, calendar time5.120.12 (0.02, 0.22)0.02 (0.01)0.31 (0.01)0.31 (0.29, 0.32)96 (1)
       Cox
      Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      , matching
      5.360.36 (0.21, 0.52)0.07 (0.01)0.45 (0.01)0.46 (0.43, 0.49)97 (1)
      Emp.SE, empirical standard error; MCSE, monte carlo standard error; RMSE, root mean squared error; UI, Uncertainty interval based on the Monte Carlo standard error of the estimated regression coefficient for bias or the mean squared error.
      a Using a time under observation scale, i.e., time 0 is when coming under observation with the current exposure status.
      b Calculated as the estimated relative risk minus the simulated incidence rate ratio.
      c Calculated as the estimated regression coefficient for the exposure minus the logarithm of the simulated incidence rate ratio.
      d No lower bound reported, due to a negative value being obtained on the mean squared error scale.
      Table 2Accrued person time, number of events and estimates obtained from the empirical data analysis
      Outcome: SARS-CoV-2 infectionOutcome: Myocarditis
      BNT162b2UnvaccinatedmRNA-1273Unvaccinated
      Individuals (N)261,964417,076488,6164,224,150
      Person years42,512196,64172,8611,786,933
      Events1,31121,3671892
      Incidence rate/1000 PY311090.250.05
      Estimators
      Poisson0.28 (0.27-0.30)(ref.)4.80 (2.90-7.95)(ref.)
      Poisson, adjusted0.13 (0.12-0.14)(ref.)3.45 (1.96-6.06)(ref.)
      Cox
      Using a time under observation scale.
      0.40 (0.38-0.43)(ref.)7.84 (3.62-17.0)(ref.)
      Cox
      Using a time under observation scale.
      , adjusted
      0.22 (0.15-0.32)(ref.)NR(ref.)
      Cox, calendar time0.14 (0.13-0.15)(ref.)2.65 (1.42-4.93)(ref.)
      Cox
      Using a time under observation scale.
      , matching
      0.13 (0.11-0.14)(ref.)NR(ref.)
      ref., reference category; PY, person years; NR, not reported due to nonconvergence of the regression; SARS-CoV-2, severe acute respiratory syndrome corona virus 2.
      a Using a time under observation scale.
      All estimators were unbiased when analyzing an event with a constant rate (Supplementary table 1).

      3.2 Empirical data analysis

      In analyses regarding the effectiveness outcome, documented SARS-CoV-2 infection, we included 417,687 individuals living in the Copenhagen municipality of whom 417,076 contributed unvaccinated person time and 261,964 contributed vaccinated person time (Table 2). Overall, 196,641 unvaccinated and 42,512 vaccinated person years were accrued. We observed 22,687 documented SARS-CoV-2 infections, of which 1,311 (incidence rate [IR] 31 events/1,000 person years) were observed among vaccinated individuals and 21,367 events among unvaccinated individuals (IR 109 events/1,000 person years), yielding an unadjusted IRR of 0.28 (95% [CI] 0.27–0.30). Adjusting for calendar time reduced the IRR to 0.13 (0.12–0.14). Likewise, estimates obtained using Cox regression and a time under observation scale were presumably upwards biased (HR 0.40, 0.38 to 0.43; HRadjusted 0.22, 0.15–0.32). Using a calendar time scale the estimated HR was 0.14 (0.13–0.15) and when matching vaccinated and unvaccinated individuals the HR was 0.13 (0.11–0.14). The incidence of SARS-CoV-2 infection in the simulated and empirical data analysis resembled each other (Figure 1, supplementary figure 2).
      For the safety outcome, a first-ever hospital diagnosis of myocarditis, we included 4,224,150 individuals who contributed 72,861 vaccinated person years and 1,786,933 unvaccinated person years. A total 110 myocarditis events were observed, 18 of these occurred during vaccinated person time (IR 25 events/100,000 person years) and 92 events during unvaccinated person time (IR 5 events/100,000 person years). All estimators yielded point estimates of the relative risk between 2.65 and 7.84 and therefore within the range of previously published results (Figure 4). Estimates from adjusted Poisson regression and the CoxCT estimator were the most conservative (IRR 3.45, 1.96–6.06; HR 2.65, 1.42–4.93), while estimates from unadjusted Poisson (IRR 4.80, 2.90–7.95) and time under observation Cox regression (HR 7.84, 3.62–17.0) were markedly higher. No estimates were obtained for Cox regression adjusted for calendar time using restricted cubic splines, as the regression did not converge. No estimates were obtained for the matching estimator, due to having observed too few events among matched unvaccinated individuals to provide a meaningful comparison.

      4. Discussion

      In this combined simulation study and empirical data analysis, we explored the susceptibility of six different estimators toward temporal biases. In the simulation study, we found that the unadjusted Poisson and Cox regression with a time under observation scale yielded biased results. Adjustment for calendar time using restricted cubic splines yielded unbiased results for Poisson regression, but vastly reduced statistical precision in Cox regression compared to all other estimators. Cox regression models using a calendar time scale and matching vaccinated and unvaccinated individuals on calendar time yielded unbiased estimates in simulations of vaccine effectiveness. In analyses of vaccine safety, the Cox regression using a calendar time scale yielded results closest to the true IRR, while the matching estimator exhibited reduced precision. Results from the empirical data analysis were comparable to previously published risk estimates.
      The main strength of our study is that we were able to test findings from the simulation study using empirical data from the Danish nationwide health registries. Although the empirical data are more complex than the simulation at hand, we did observe the same relationship between estimators and estimates as in the simulation study, making it more likely that we did simulate parameters close to the actual data.
      A limitation of our study is that it assumes that the relative effect of vaccination is constant over time. Vaccine effectiveness toward documented SARS-CoV-2 infection decreases over time [
      • Haas E.J.
      • Angulo F.J.
      • McLaughlin J.M.
      • Anis E.
      • Singer S.R.
      • Khan F.
      • et al.
      Impact and effectiveness of mRNA BNT162b2 vaccine against SARS-CoV-2 infections and COVID-19 cases, hospitalisations, and deaths following a nationwide vaccination campaign in Israel: an observational study using national surveillance data.
      ], and adverse events are most likely to occur within days of vaccination [
      • Patone M.
      • Mei X.W.
      • Handunnetthi L.
      • Dixon S.
      • Zaccardi F.
      • Shankar-Hari M.
      • et al.
      Risks of myocarditis, pericarditis, and cardiac arrhythmias associated with COVID-19 vaccination or SARS-CoV-2 infection.
      ]. We chose not to simulate such time-varying effects, as this would complicate the simulation beyond the actual question of interest and complicate the evaluation of different estimators, as a single risk estimate would no longer suffice to quantify bias. Further, a limitation of our empirical data analysis is that the true vaccine effectiveness and risk increase about myocarditis is unknown, although multiple studies on vaccine effectiveness are in agreement with our findings [
      • Haas E.J.
      • Angulo F.J.
      • McLaughlin J.M.
      • Anis E.
      • Singer S.R.
      • Khan F.
      • et al.
      Impact and effectiveness of mRNA BNT162b2 vaccine against SARS-CoV-2 infections and COVID-19 cases, hospitalisations, and deaths following a nationwide vaccination campaign in Israel: an observational study using national surveillance data.
      ,
      • Tartof S.Y.
      • Slezak J.M.
      • Fischer H.
      • Hong V.
      • Ackerson B.K.
      • Ranasinghe O.N.
      • et al.
      Effectiveness of mRNA BNT162b2 COVID-19 vaccine up to 6 months in a large integrated health system in the USA: a retrospective cohort study.
      ,
      • Reis B.Y.
      • Barda N.
      • Leshchinsky M.
      • Kepten E.
      • Hernán M.A.
      • Lipsitch M.
      • et al.
      Effectiveness of BNT162b2 vaccine against delta variant in adolescents.
      ,
      • Dagan N.
      • Barda N.
      • Kepten E.
      • Miron O.
      • Perchik S.
      • Katz M.A.
      • et al.
      BNT162b2 mRNA covid-19 vaccine in a nationwide mass vaccination setting.
      ]. For the event of myocarditis, the true effect size is much more uncertain, making it difficult to determine which estimator may be the least biased in the empirical setting.
      A previous simulation study evaluated the calendar time scale in Cox regression analyses in the presence of a time-varying exposure with a linear trend, but found no benefits compared to an age or time under observation scale [
      • Griffin B.A.
      • Anderson G.L.
      • Shih R.A.
      • Whitsel E.A.
      Use of alternative time scales in Cox proportional hazard models: implications for time-varying environmental exposures.
      ]. Differences in results between our and the previous study could be explained by the nonlinear trends evaluated in our study, demanding more flexible methods of adjusting for calendar time. A popular method to adjust for calendar time in analyses of COVID-19 vaccine effectiveness or safety is matching individuals on calendar time [
      • Dagan N.
      • Barda N.
      • Kepten E.
      • Miron O.
      • Perchik S.
      • Katz M.A.
      • et al.
      BNT162b2 mRNA covid-19 vaccine in a nationwide mass vaccination setting.
      ], which was also tested in this study. Matching ensures that exposed and unexposed person time has been accrued simultaneously, but it may be difficult to find suitable unvaccinated controls in a population with rapid vaccine uptake, effectively reducing the size of the study population and statistical precision, due to discarding unvaccinated person time. This was inconsequential in analyses with a common outcome (vaccine effectiveness) but reduced precision in analyses where the outcome was rare (vaccine safety). Finally, individuals delaying COVID-19 vaccination may differ on important, unmeasured variables from individuals who are immunized as soon as possible.
      Our findings underline the need to use flexible methods to adjust for calendar time in analyses regarding COVID-19 vaccine effectiveness and safety. This can be achieved parametrically, e.g., Poisson regression adjusted for calendar time modeled using restricted cubic splines, or nonparametrically by choosing calendar time as the underlying time scale in a Cox regression model. Currently, not all studies report the time scale used in Cox regression analyses. This may be an indicator of applied researchers not being aware of the implications and importance of the choice of time scale.

      5. Conclusions

      We evaluated the performance of Poisson and Cox regression with and without adjustment for calendar time in analyses regarding COVID-19 vaccine effectiveness and safety. Cox regression using a calendar time scale consistently yielded the least biased results when analyzing simulated data, and provided effect estimates close to previously reported estimates when analyzing empirical data. Matching vaccinated and unvaccinated individuals yielded unbiased results in vaccine effectiveness analyses but reduces precision when analyzing rare outcomes. More analyses are needed regarding the performance of these estimators in the presence of time-varying effects and the implications of nonproportional hazards.

      Supplementary data

      References

        • Haas E.J.
        • Angulo F.J.
        • McLaughlin J.M.
        • Anis E.
        • Singer S.R.
        • Khan F.
        • et al.
        Impact and effectiveness of mRNA BNT162b2 vaccine against SARS-CoV-2 infections and COVID-19 cases, hospitalisations, and deaths following a nationwide vaccination campaign in Israel: an observational study using national surveillance data.
        Lancet. 2021; 397: 1819-1829
        • Husby A.
        • Hansen J.V.
        • Fosbøl E.
        • Thiesson E.M.
        • Madsen M.
        • Thomsen R.W.
        • et al.
        SARS-CoV-2 vaccination and myocarditis or myopericarditis: population based cohort study.
        BMJ. 2021; 375: e068665
        • Rosen B.
        • Waitzberg R.
        • Israeli A.
        Israel’s rapid rollout of vaccinations for COVID-19.
        Isr J Health Policy Res. 2021; 10: 6
      1. Statistics » COVID-19 Vaccinations.
        (Available at:)
        • Reilev M.
        • Olesen M.
        • Kildegaard H.
        • Stovring H.
        • Andersen J.H.
        • Hallas J.
        • et al.
        Changing characteristics over time of individuals receiving COVID-19 vaccines in Denmark: A population-based descriptive study of vaccine uptake.
        Scand J Public Health. 2022; 50: 686-692
        • Griffin B.A.
        • Anderson G.L.
        • Shih R.A.
        • Whitsel E.A.
        Use of alternative time scales in Cox proportional hazard models: implications for time-varying environmental exposures.
        Statist Med. 2012; 31: 3320-3327
        • Rothman K.J.
        • Lash T.L.
        • VanderWeele T.J.
        • Haneuse S.
        Modern epidemiology.
        (Fourth edition) Wolters Kluwer/Lippincott Williams & Wilkins, Philadelphia2021
        • Ritchie H.
        • Mathieu E.
        • Rodés-Guirao L.
        • Appel C.
        • Giattino C.
        • Ortiz-Ospina E.
        • et al.
        Coronavirus pandemic (COVID-19).
        (Available at:)
        https://ourworldindata.org/covid-cases
        Date accessed: February 3, 2022
        • Patone M.
        • Handunnetthi L.
        • Saatci D.
        • Pan J.
        • Katikireddi S.V.
        • Razvi S.
        • et al.
        Neurological complications after first dose of COVID-19 vaccines and SARS-CoV-2 infection.
        Nat Med. 2021; 27: 2144-2153
        • Patone M.
        • Mei X.W.
        • Handunnetthi L.
        • Dixon S.
        • Zaccardi F.
        • Shankar-Hari M.
        • et al.
        Risks of myocarditis, pericarditis, and cardiac arrhythmias associated with COVID-19 vaccination or SARS-CoV-2 infection.
        Nat Med. 2022; 28: 410-422
        • Karlstad Ø.
        • Hovi P.
        • Husby A.
        • Härkänen T.
        • Selmer R.M.
        • Pihlström N.
        • et al.
        SARS-CoV-2 vaccination and myocarditis in a nordic cohort study of 23 million residents.
        JAMA Cardiol. 2022; 7: 600-612
        • Li X.
        • Raventós B.
        • Roel E.
        • Pistillo A.
        • Martinez-Hernandez E.
        • Delmestri A.
        • et al.
        Association between covid-19 vaccination, SARS-CoV-2 infection, and risk of immune mediated neurological events: population based cohort and self-controlled case series analysis.
        BMJ. 2022; 376: e068373
        • Hansen C.H.
        • Schelde A.B.
        • Moustsen-Helm I.R.
        • Emborg H.D.
        • Krause T.G.
        • Mølbak K.
        • et al.
        Vaccine effectiveness against SARS-CoV-2 infection with the Omicron or Delta variants following a two-dose or booster BNT162b2 or mRNA-1273 vaccination series: A Danish cohort study.
        (Available at:)
        • Tartof S.Y.
        • Slezak J.M.
        • Fischer H.
        • Hong V.
        • Ackerson B.K.
        • Ranasinghe O.N.
        • et al.
        Effectiveness of mRNA BNT162b2 COVID-19 vaccine up to 6 months in a large integrated health system in the USA: a retrospective cohort study.
        Lancet. 2021; 398: 1407-1416
        • Harrell FE Jr.,
        Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis.
        Springer International Publishing, Cham2015
        • Lund J.L.
        • Richardson D.B.
        • Stürmer T.
        The active comparator, new user study design in pharmacoepidemiology: historical foundations and contemporary application.
        Curr Epidemiol Rep. 2015; 2: 221-228
        • Suissa S.
        Immortal time bias in pharmaco-epidemiology.
        Am J Epidemiol. 2008; 167: 492-499
        • Morris T.P.
        • White I.R.
        • Crowther M.J.
        Using simulation studies to evaluate statistical methods.
        Stat Med. 2019; 38: 2074-2102
        • Pedersen C.B.
        The Danish civil registration system.
        Scand J Public Health. 2011; 39: 22-25
        • Voldstedlund M.
        • Haarh M.
        • Mølbak K.
        Representatives the MB of. The Danish microbiology database (MiBa) 2010 to 2013.
        Eurosurveillance. 2014; 19: 20667
        • Grove Krause T.
        • Jakobsen S.
        • Haarh M.
        • Mølbak K.
        The Danish vaccination register.
        Euro Surveill. 2012; 17: 20155
        • Polack F.P.
        • Thomas S.J.
        • Kitchin N.
        • Absalon J.
        • Gurtman A.
        • Lockhart S.
        • et al.
        Safety and Efficacy of the BNT162b2 mRNA covid-19 vaccine.
        New Engl J Med. 2020; 383: 2603-2615
        • Baden L.R.
        • El Sahly H.M.
        • Essink B.
        • Kotloff K.
        • Frey S.
        • Novak R.
        • et al.
        Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine.
        New Engl J Med. 2021; 384: 403-416
        • Reis B.Y.
        • Barda N.
        • Leshchinsky M.
        • Kepten E.
        • Hernán M.A.
        • Lipsitch M.
        • et al.
        Effectiveness of BNT162b2 vaccine against delta variant in adolescents.
        N Engl J Med. 2021; 385: 2101-2103
        • Dagan N.
        • Barda N.
        • Kepten E.
        • Miron O.
        • Perchik S.
        • Katz M.A.
        • et al.
        BNT162b2 mRNA covid-19 vaccine in a nationwide mass vaccination setting.
        N Engl J Med. 2021; 384: 1412-1423