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

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.

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.

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 [1], hospitalization for COVID-19 or adverse events related to vaccination [2], 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 [3], the United Kingdom [4] and Scandinavia [5], much of the unvaccinated and vaccinated person time has not been accrued simultaneously. This increases the risk of temporal biases affecting study results [6], 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 [7]. 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.

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.

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

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.

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 [8], 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 e syndrome which have previously been associated with COVID-19 immunization [2,9e12].
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.
Event times for each outcome were obtained by performing a Bernoulli trial for each individual and day based on the above-described probabilities.

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.

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 [13], and vaccine effectiveness against the Delta variant [14]. 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 [2,10,11]. All simulated effects were constant over time.

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 [15]. 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 Cox CT ). 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 [16,17], 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(b exposure ) e IRR, and also on the log-scale, i.e., b exposure elog(IRR). We further calculated Monte Carlo standard errors and 95% Monte Carlo uncertainty intervals for all performance metrics [18].

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 [19], COVID-19 test results [20] and information on COVID-19 immunization from the Danish vaccination register [21]. 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 [22], 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%e95% [1,14].
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 [23], 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 [2,10,11]. Myocarditis and pericarditis are currently the only serious adverse events of the BNT162b2 and mRNA- 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/.

Simulation
In simulated analyses of vaccine effectiveness with a true IRR of 0.1, we obtained a mean point estimate of 0.26 (bias 5 estimated IRRetrue IRR 5 0.16; log bias 5 log(IRR) e log(true IRR) 5 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, Cox CT and matching estimators provided valid 95% CIs (coverage probability of 95%). When evaluating a true IRR of 0.5, adjusted Poisson, Cox CT 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 Cox CT 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 Cox CT 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.
For detailed results of all scenarios, see Table 1. All estimators were unbiased when analyzing an event with a constant rate (Supplementary table 1).

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) HR adjusted 0.22, 0.15e0.32). Using a calendar time scale the estimated HR was 0.14 (0.13e0.15) and when matching vaccinated and unvaccinated individuals the HR was 0.13 (0.11e0.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 Cox CT estimator were the most conservative (IRR 3.45, 1.96e6.06; HR 2.65, 1.42e4.93), while estimates from unadjusted Poisson (IRR 4.80, 2.90e7.95) and time under observation Cox regression (HR 7.84, 3.62e17.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.

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 [1], and adverse events are most likely to occur within days of vaccination [10]. 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 [1,14,24,25]. 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 timevarying exposure with a linear trend, but found no benefits compared to an age or time under observation scale [6]. 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 [25], 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.

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.